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