/******************************************************************** * This program solves a Toeplitz linear system: * * N * * S R x = y (i=1,...,N) * * j=1 N+i-j j i * * ----------------------------------------------------------------- * * SAMPLE RUN: * * N=10 * * R = (-10,-12,3,4,5,6,7,8,9,5,11,12,13,14,15,16,17,18,19) * * Y = (1,1,1,1,1,1,1,1,30,1) * * * * The solution vector is: * * 1 4.145258 * * 2 -7.202165 * * 3 0.136163 * * 4 0.473885 * * 5 0.811607 * * 6 1.149329 * * 7 1.487051 * * 8 1.824773 * * 9 -3.637505 * * 10 2.500217 * * * * ----------------------------------------------------------------- * * Reference: "Numerical Recipes by W.H. Press, B.P. Flannery, S.A. * * Teukolsky, W.T. Vetterling, Cambridge University * * Press, 1987" * * * * C++ Version By J-P Moreau, Paris * * (www.jpmoreau.fr) * ********************************************************************/ #include void TOEPLZ(double *R, double *X, double *Y, int N) { /*-------------------------------------------------------------------- * This program solves a Toeplitz linear system: * * N * * S R x = y (i=1,...,N) * * j=1 N+i-j j i * * * * Input consists of vectors R of size 2N-1 and Y of size N. * * Output is solution vector X of size N. * * N is the size of the linear system. * --------------------------------------------------------------------*/ const int NMAX=100; double ZERO=0.0; // NMAX is the maximum expected value of N. double G[NMAX], H[NMAX]; double PP, PT1, PT2, QQ, QT1, QT2, SD, SGD, SGN, SHN, SXN; int J, K, M, M1, M2; if (R[N] == ZERO) goto e99; X[1]=Y[1]/R[N]; // initialize for the recursion. if(N == 1) return; G[1]=R[N-1]/R[N]; H[1]=R[N+1]/R[N]; for (M=1; M<=N; M++) { // main loop over the recursion. M1=M+1; SXN=-Y[M1]; SD=-R[N]; for (J=1; J<=M; J++) { SXN=SXN+R[N+M1-J]*X[J]; SD=SD + R[N+M1-J]*G[M-J+1]; } if (SD == ZERO) goto e99; X[M1]=SXN/SD; // whence x. for (J=1; J<=M; J++) X[J] -= X[M1]*G[M-J+1]; if (M1 == N) return; // normal exit; SGN=-R[N-M1]; // compute numerator and denominator for G and H. SHN=-R[N+M1]; SGD=-R[N]; for (J=1; J<=M; J++) { SGN += R[N+J-M1]*G[J]; SHN += R[N-J+M1]*H[J]; SGD += R[N+J-M1]*H[M-J+1]; } if(SD == ZERO || SGD == ZERO) goto e99; G[M1]=SGN/SGD; // whence G and H H[M1]=SHN/SD; K=M; M2=(M+1)/2; PP=G[M1]; QQ=H[M1]; for (J=1; J<=M2; J++) { PT1=G[J]; PT2=G[K]; QT1=H[J]; QT2=H[K]; G[J]=PT1-PP*QT2; G[K]=PT2-PP*QT1; H[J]=QT1-QQ*PT2; H[K]=QT2-QQ*PT1; K--; } } // M loop e99:printf(" Levinson method fails: singular principal minor.\n"); } // TOEPLZ // main program void main() { int i, N; double R[20] = { 0.0,-10.0, -12.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 5.0, 11.0, 12.0, 13.0,14.0,15.0,16.0,17.0,18.0,19.0}; double Y[11] = { 0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,30.0, 1.0}; double X[11]; N = 10; TOEPLZ(R,X,Y,N); printf("\n The solution vector is:\n"); for (i=1; i<=N; i++) printf(" %2d %f\n", i, X[i]); } // end of file toeplitz.cpp