00001
00002
00003 #include "ctypes_def.h"
00004 #include "cblas.h"
00005
00006 #if TYPE == FLOAT
00007 #define LEVINSON slevinson
00008 #elif TYPE == DOUBLE
00009 #define LEVINSON dlevinson
00010 #elif TYPE == COMPLEXFLOAT
00011 #define LEVINSON clevinson
00012 #elif TYPE == COMPLEXDOUBLE
00013 #define LEVINSON zlevinson
00014 #endif
00015
00016 #ifdef NOBLAS
00017 void LEVINSON( int n, CTYPE *r, CTYPE *b, CTYPE *x)
00018 {
00019 int k,i;
00020 CTYPE alpha, mu;
00021 CTYPE *b_tmp;
00022 RTYPE beta;
00023
00024 if (cimag(r[0]) != RZERO || creal(r[0]) <= RZERO )
00025 {
00026 fprintf(stderr, "%s: Toeplitz matrix must be Hermitian.\n",__func__);
00027 fprintf(stderr, "imag = %f, real = %f\n",cimag(r[0]),creal(r[0]));
00028 abort();
00029 }
00030 b_tmp = (CTYPE *)malloc(n*sizeof(CTYPE));
00031 beta = creal(r[0]);
00032 x[0] = b[0] / r[0];
00033 b[0] = -r[1]/beta;
00034 alpha = b[0];
00035 for (k=1; k<=n-1; k++)
00036 {
00037 beta = (RONE - CONJ(alpha)*alpha)*beta;
00038 mu = b[k];
00039 for (i=2; i<=(k+1); i++)
00040 mu -= conj(r[i-1])*x[k+1-i];
00041 mu /= beta;
00042 for (i=1; i<=k; i++)
00043 x[i-1] += mu*b[k-i];
00044 x[k] = mu;
00045
00046 if (k<n-1)
00047 {
00048 alpha = -r[k+1];
00049 for (i=2; i<=(k+1); i++)
00050 alpha -= r[i-1] * b[k+1-i];
00051 alpha /= beta;
00052 memcpy(b_tmp, b, k*sizeof(CTYPE));
00053 for (i=1; i<=k; i++)
00054 b[i-1] += alpha * CONJ(b_tmp[k-i]);
00055 b[k] = alpha;
00056 }
00057 }
00058 free(b_tmp);
00059 }
00060
00061 #else
00062
00063 void LEVINSON( int n, CTYPE *r, CTYPE *b, CTYPE *x)
00064 {
00065 int k,i;
00066 CTYPE alpha, mu;
00067 CTYPE *b_tmp;
00068 RTYPE beta;
00069
00070 if (cimag(r[0]) != RZERO || creal(r[0]) <= RZERO )
00071 {
00072 fprintf(stderr, "%s: Toeplitz matrix must be Hermitian.\n",__func__);
00073 fprintf(stderr, "imag = %f, real = %f\n",cimag(r[0]),creal(r[0]));
00074 abort();
00075 }
00076 b_tmp = (CTYPE *)malloc(n*sizeof(CTYPE));
00077 beta = creal(r[0]);
00078 x[0] = b[0] / r[0];
00079 b[0] = -r[1]/beta;
00080 alpha = b[0];
00081 for (k=1; k<=n-1; k++)
00082 {
00083 beta = (RONE - CONJ(alpha)*alpha)*beta;
00084 mu = b[k] - DOTC(k,&r[1],1,x,-1);
00085 mu /= beta;
00086
00087 AXPY(k,mu,b,-1,x,1);
00088 x[k] = mu;
00089
00090 if (k<n-1)
00091 {
00092 alpha = -r[k+1] - DOTU(k,&r[1],1,b,-1);
00093 alpha /= beta;
00094 for (i=0; i<k; i++)
00095 b_tmp[i] = CONJ(b[i]);
00096 AXPY(k,alpha,b_tmp,-1,b,1);
00097 b[k] = alpha;
00098 }
00099 }
00100 free(b_tmp);
00101 }
00102
00103 #endif
00104 #undef LEVINSON
00105 #undef AXPY
00106 #undef DOTU
00107 #undef DOTC
00108 #include "ctypes_undef.h"
