A Toeplitz matrix is an n Times n matrix A = (a_ij) such that a_ij = a_i - 1, j

Question

A Toeplitz matrix is an n Times n matrix A = (a_ij) such that a_ij = a_i - 1, j - 1 for i = 2, 3, , n and j = 2, 3, , n (in other words, entries remain constant along each diagonal of the matrix). Give an algorithm that multiplies an n Times n Toeplitz matrix by a vector of size n in O(n log n) time. To achieve this, you need to think of an efficient way to represent the matrix, and you will have to use the FFT algorithm. Write a program that implements the "generic" O(n^2) matrix-vector multiplication, and a program that implements your algorithm above. Run your code on random samples of size n = 10^k for k = 2, 3, 4, 5, 6. Write a report describing the results of your experiments. You are strongly encouraged to code the entire assigment from scratch. However, you are allowed to use "off-the-shelf" implementations of the FFT algorithm, provided you disclose any source you use for the solution of this assignment.

Explanation / Answer

#include<math.h>

#include<"mes.h">

function h=topsolver(r,q);

r*h=q

void mfunction(

int n,array *p[],int nr,

const array *p[])

}

double *A ,h* ;

int j,k;

double eps = mxGetEps();
   int n = (mxGetN(prhs[0])>=mxGetM(prhs[0])) ? mxGetN(prhs[0]) :
mxGetM(prhs[0]) ;
   double *r = mxGetPr(prhs[0]);
   double *q = mxGetPr(prhs[1]);
   double alpha = r[0];

   n = n - 1;

   plhs[0] = mxCreateDoubleMatrix(n+1,1,0);
   h = mxGetPr(plhs[0]);

   h[0] = q[0]/r[0];

   a = mxCalloc((n+1)*(n+1),sizeof(double));
   if (a == NULL) {
      (".");
   }

   a[(0*(n+1))+0] = 1.0;

   for (k = 1; k <= 6; k++) {
       a[(k*(n+1))+k-1] = 0;
       a[(0*(n+1))+k] = 1.0;
       beta = 0.0;
       for (j = 0; j <= k-1; j++) {
           beta += r[k-j]*a[(j*(n+1))+k-1];
       }
       beta /= alpha;
       for (j = 1; j <= k; j++) {
           a[(j*(n+1))+k] = a[(j*(n+1))+k-1] - beta*a[((k-j)*(n+1))+k-1];
       }
       alpha *= (1 - beta*beta);
       h[k] = q[k];
       for (j = 0; j <= k-1; j++) {
           h[k] -= r[k-j]*h[j];
       }
       h[k] /= alpha;
       for (j = 0; j <= k-1; j++) {
           h[j] += a[((k-j)*(n+1))+k]*h[k];
       }

   mxFree(a);

   return;
}

Related Questions

Navigate

• Browse (All)
• Browse A
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.