00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef _LA_NONCLASS_H_
00019 #define _LA_NONCLASS_H_
00020 #include "vec.h"
00021 #include "smat.h"
00022 #include "mat.h"
00023 #include "la_traits.h"
00024
00025 namespace LA {
00026
00027
00028 template <class T>
00029 const NRSMat<T> twoside_transform(const NRSMat<T> &S, const NRMat<T> &C, bool transp=0)
00030 {
00031 if(transp)
00032 {
00033 NRMat<T> tmp = C * S;
00034 NRMat<T> result(C.nrows(),C.nrows());
00035 result.gemm((T)0,tmp,'n',C,'t',(T)1);
00036 return NRSMat<T>(result);
00037 }
00038 NRMat<T> tmp = S * C;
00039 NRMat<T> result(C.ncols(),C.ncols());
00040 result.gemm((T)0,C,'t',tmp,'n',(T)1);
00041 return NRSMat<T>(result);
00042 }
00043
00044
00045
00046
00047 template <class T>
00048 const NRMat<T> diagonalmatrix(const NRVec<T> &x)
00049 {
00050 int n=x.size();
00051 NRMat<T> result((T)0,n,n);
00052 T *p = result[0];
00053 for(int j=0; j<n; j++) {*p = x[j]; p+=(n+1);}
00054 return result;
00055 }
00056
00057
00058
00059 template<class T>
00060 inline const NRMat<T> commutator ( const NRMat<T> &x, const NRMat<T> &y, const bool trx=0, const bool tryy=0)
00061 {
00062 NRMat<T> r(trx?x.ncols():x.nrows(), tryy?y.nrows():y.ncols());
00063 r.gemm((T)0,x,trx?'t':'n',y,tryy?'t':'n',(T)1);
00064 r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)-1);
00065 return r;
00066 }
00067
00068
00069 template<class T>
00070 inline const NRMat<T> anticommutator ( const NRMat<T> &x, const NRMat<T> &y, const bool trx=0, const bool tryy=0)
00071 {
00072 NRMat<T> r(trx?x.ncols():x.nrows(), tryy?y.nrows():y.ncols());
00073 r.gemm((T)0,x,trx?'t':'n',y,tryy?'t':'n',(T)1);
00074 r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)1);
00075 return r;
00076 }
00077
00078
00079
00080
00082
00084
00085 #define declare_la(T) \
00086 extern const NRVec<T> diagofproduct(const NRMat<T> &a, const NRMat<T> &b,\
00087 bool trb=0, bool conjb=0); \
00088 extern T trace2(const NRMat<T> &a, const NRMat<T> &b, bool trb=0); \
00089 extern T trace2(const NRSMat<T> &a, const NRSMat<T> &b, const bool diagscaled=0);\
00090 extern T trace2(const NRSMat<T> &a, const NRMat<T> &b, const bool diagscaled=0);\
00091 extern void linear_solve(NRMat<T> &a, NRMat<T> *b, double *det=0,int n=0); \
00092 extern void linear_solve(NRSMat<T> &a, NRMat<T> *b, double *det=0, int n=0); \
00093 extern void linear_solve(NRMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
00094 extern void linear_solve(NRSMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
00095 extern void diagonalize(NRMat<T> &a, NRVec<LA_traits<T>::normtype> &w, const bool eivec=1, const bool corder=1, int n=0, NRMat<T> *b=NULL, const int itype=1); \
00096 extern void diagonalize(NRSMat<T> &a, NRVec<LA_traits<T>::normtype> &w, NRMat<T> *v, const bool corder=1, int n=0, NRSMat<T> *b=NULL, const int itype=1);\
00097 extern void singular_decomposition(NRMat<T> &a, NRMat<T> *u, NRVec<T> &s,\
00098 NRMat<T> *v, const bool corder=1, int m=0, int n=0);
00099
00100
00101
00102 declare_la(double)
00103 declare_la(complex<double>)
00104
00105
00106
00107 extern void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
00108 NRMat<double> *vl, NRMat<double> *vr, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
00109 NRMat<double> *b=NULL, NRVec<double> *beta=NULL);
00110 extern void gdiagonalize(NRMat<double> &a, NRVec< complex<double> > &w,
00111 NRMat< complex<double> >*vl, NRMat< complex<double> > *vr,
00112 const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
00113 NRMat<double> *b=NULL, NRVec<double> *beta=NULL);
00114 extern NRMat<double> matrixfunction(NRSMat<double> a, double (*f) (double));
00115 extern NRMat<double> realmatrixfunction(NRMat<double> a, double (*f) (double));
00116 extern NRMat<complex<double> > complexmatrixfunction(NRMat<double> a, double (*fre) (double), double (*fim) (double));
00117 extern NRMat<double> matrixfunction(NRMat<double> a, complex<double> (*f)(const complex<double> &),const bool adjust=0);
00118
00119 extern complex<double> sqrtinv(const complex<double> &);
00120 extern double sqrtinv(const double);
00121
00122
00123 inline NRMat<double> sqrt(const NRSMat<double> &a) { return matrixfunction(a,&std::sqrt); }
00124 inline NRMat<double> sqrtinv(const NRSMat<double> &a) { return matrixfunction(a,&sqrtinv); }
00125 inline NRMat<double> realsqrt(const NRMat<double> &a) { return realmatrixfunction(a,&std::sqrt); }
00126 inline NRMat<double> realsqrtinv(const NRMat<double> &a) { return realmatrixfunction(a,&sqrtinv); }
00127 inline NRMat<double> log(const NRSMat<double> &a) { return matrixfunction(a,&std::log); }
00128 extern NRMat<double> log(const NRMat<double> &a);
00129 extern NRMat<double> exp0(const NRMat<double> &a);
00130
00131
00132 extern const NRMat<double> realpart(const NRMat< complex<double> >&);
00133 extern const NRMat<double> imagpart(const NRMat< complex<double> >&);
00134 extern const NRMat< complex<double> > realmatrix (const NRMat<double>&);
00135 extern const NRMat< complex<double> > imagmatrix (const NRMat<double>&);
00136 extern const NRMat< complex<double> > complexmatrix (const NRMat<double>&, const NRMat<double>&);
00137
00138
00139 extern void cholesky(NRMat<double> &a, bool upper=1);
00140 extern void cholesky(NRMat<complex<double> > &a, bool upper=1);
00141
00142
00143 template<typename T>
00144 const NRMat<T> inverse(NRMat<T> a, T *det=0)
00145 {
00146 #ifdef DEBUG
00147 if(a.nrows()!=a.ncols()) laerror("inverse() for non-square matrix");
00148 #endif
00149 NRMat<T> result(a.nrows(),a.nrows());
00150 result = (T)1.;
00151 a.copyonwrite();
00152 linear_solve(a, &result, det);
00153 result.transposeme();
00154 return result;
00155 }
00156
00157
00158 template<class MAT>
00159 typename LA_traits<MAT>::normtype MatrixNorm(const MAT &A, const char norm);
00160
00161
00162 template<class MAT>
00163 typename LA_traits<MAT>::normtype CondNumber(const MAT &A, const char norm);
00164
00165
00166
00167 template<class MAT>
00168 const typename LA_traits<MAT>::elementtype determinant(MAT a)
00169 {
00170 typename LA_traits<MAT>::elementtype det;
00171 if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
00172 linear_solve(a,NULL,&det);
00173 return det;
00174 }
00175
00176
00177 template<class MAT>
00178 const typename LA_traits<MAT>::elementtype determinant_destroy(MAT &a)
00179 {
00180 typename LA_traits<MAT>::elementtype det;
00181 if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
00182 linear_solve(a,NULL,&det);
00183 return det;
00184 }
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 template<class T>
00202 int linear_solve_x(NRMat<T> &A, T *B, const int rhsCount, const int eqCount, const bool eq, const bool saveA, double *rcond);
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 template<class T>
00218 int multiply_by_inverse(NRMat<T> &A, NRMat<T> &B, bool useEq, double *rcond);
00219
00220
00221
00222
00223
00224
00225
00226 template<class MAT, class INDEX>
00227 const NRMat<typename LA_traits<MAT>::elementtype> submatrix(const MAT a, const int nrows, const INDEX rows, const int ncols, const INDEX cols, int indexshift=0, bool ignoresign=false, bool equalsigns=false)
00228 {
00229 NRMat<typename LA_traits<MAT>::elementtype> r(nrows,ncols);
00230
00231 if(equalsigns)
00232 {
00233 if(ignoresign)
00234 {
00235 for(int i=0; i<nrows; ++i)
00236 for(int j=0; j<ncols; ++j)
00237 r(i,j) = rows[i]*cols[j]<0?0.:a(std::abs(rows[i])+indexshift,std::abs(cols[j])+indexshift);
00238 }
00239 else
00240 {
00241 for(int i=0; i<nrows; ++i)
00242 for(int j=0; j<ncols; ++j)
00243 r(i,j) = rows[i]*cols[j]<0?0.:a(rows[i]+indexshift,cols[j]+indexshift);
00244 }
00245 }
00246 else
00247 {
00248 if(ignoresign)
00249 {
00250 for(int i=0; i<nrows; ++i)
00251 for(int j=0; j<ncols; ++j)
00252 r(i,j) = a(std::abs(rows[i])+indexshift,std::abs(cols[j])+indexshift);
00253 }
00254 else
00255 {
00256 for(int i=0; i<nrows; ++i)
00257 for(int j=0; j<ncols; ++j)
00258 r(i,j) = a(rows[i]+indexshift,cols[j]+indexshift);
00259 }
00260 }
00261
00262 return r;
00263 }
00264
00265
00266
00267 extern void adjustphases(NRMat<double> &v);
00268
00269
00270
00271
00272
00273 template<class T> inline void xcopy(int n, const T *x, int incx, T *y, int incy);
00274 template<class T> inline void xaxpy(int n, const T &a, const T *x, int incx, T *y, int incy);
00275 template<class T> inline T xdot(int n, const T *x, int incx, const T *y, int incy);
00276
00277
00278
00279 template<>
00280 inline void xcopy<double> (int n, const double *x, int incx, double *y, int incy)
00281 {
00282 cblas_dcopy(n, x, incx, y, incy);
00283 }
00284
00285 template<>
00286 inline void xaxpy<double>(int n, const double &a, const double *x, int incx, double *y, int incy)
00287 {
00288 cblas_daxpy(n, a, x, incx, y, incy);
00289 }
00290
00291 template<>
00292 inline double xdot<double>(int n, const double *x, int incx, const double *y, int incy)
00293 {
00294 return cblas_ddot(n,x,incx,y,incy);
00295 }
00296
00297
00298
00299
00300
00301 template<typename M, typename T>
00302 NRMat<T> reconstructmatrix(const M &implicitmat)
00303 {
00304 NRMat<T> r(implicitmat.nrows(),implicitmat.ncols());
00305 NRVec<T> rhs(0.,implicitmat.ncols());
00306 NRVec<T> tmp(implicitmat.nrows());
00307 for(int i=0; i<implicitmat.ncols(); ++i)
00308 {
00309 rhs[i]=1.;
00310 implicitmat.gemv(0.,tmp,'n',1.,rhs);
00311 for(int j=0; j<implicitmat.nrows(); ++j) r(j,i)=tmp[j];
00312 rhs[i]=0.;
00313 }
00314 return r;
00315 }
00316
00317
00318 }
00319 #endif