00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00028
00029 #include <math.h>
00030
00031
00032 #include "toeplitz.h"
00033
00042 void toeplitz_gershgorin(const T *t, unsigned int n, T *low, T *up) {
00043
00044
00045 if(!t || !low || !up) return;
00046
00047
00048 int i;
00049 T sum = 0.0;
00050 for(i=1; i<n; ++i) sum += fabs(t[i]);
00051
00052 T rad = sum;
00053 for(i=1; i<(n+1)/2; ++i) {
00054 sum = sum - fabs(t[n-i]) + fabs(t[i]);
00055 if(sum > rad) rad = sum;
00056 }
00057
00058 *low = t[0] - rad;
00059 *up = t[0] + rad;
00060 }
00061
00070 void toeplitz_mult(int n, const T *t, const T *x, T *y) {
00071
00072
00073 if(n <= 0) return;
00074
00075
00076 memset(y, 0, n * sizeof(T));
00077 _axpy(n, t[0], x, 1, y, 1);
00078
00079
00080 unsigned int i;
00081 for(i=1; i<n; ++i) {
00082 _axpy(n-i, t[i], x + i, 1, y, 1);
00083 _axpy(n-i, t[i], x, 1, y + i, 1);
00084 }
00085 }
00086
00097 void levinson(unsigned int n, T td, const T *t0, const T *b0, T *x) {
00098
00099
00100 T inv_td = (T)1.0 / td;
00101 T *b = (T *) malloc(n * sizeof(T));
00102 memcpy(b, b0, n * sizeof(T));
00103 _scal(n, inv_td, b, 1);
00104
00105
00106 T *t = (T *) malloc((n-1) * sizeof(T));
00107 memcpy(t, t0, (n-1) * sizeof(T));
00108 _scal(n-1, inv_td, t, 1);
00109
00110 T *y = (T *) malloc(n * sizeof(T));
00111 x[0] = b[0];
00112 y[0] = -t[0];
00113
00114 T alpha = -t[0];
00115 T beta = (T)1.0;
00116
00117 T *aux = (T *) malloc(n * sizeof(T));
00118
00119 int k;
00120 for(k=1; k<n; ++k) {
00121 beta = ((T)1.0 - alpha * alpha) * beta;
00122 T mu = (b[k] - _dot(k, t, 1, x, -1)) / beta;
00123 _axpy(k, mu, y, -1, x, 1);
00124 x[k] = mu;
00125
00126 if(k < n-1) {
00127 alpha = -(t[k] + _dot(k, t, 1, y, -1)) / beta;
00128
00129
00130 memcpy(aux, y, k * sizeof(T));
00131 _axpy(k, alpha, aux, -1, y, 1);
00132 y[k] = alpha;
00133 }
00134 }
00135
00136 free(b);
00137 free(t);
00138 free(y);
00139 free(aux);
00140 }
00141