00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "matrix.cpp"
00027
00028 namespace PLib {
00029
00030 double Matrix<Complex>::norm(void){
00031 int i,j ;
00032 double sumR, sumI, maxsum;
00033 int init=0 ;
00034 Complex *ptr ;
00035 ptr = m-1 ;
00036 maxsum = -1 ;
00037 for(i=0;i<rows();++i){
00038 sumR = 0.0 ;
00039 sumI = 0.0 ;
00040 for ( j = 0; j < cols(); ++j) {
00041 sumR += real(*ptr)*real(*ptr) ;
00042 sumI += imag(*ptr)*imag(*ptr) ;
00043 }
00044 if(init)
00045 maxsum = (maxsum>(sumR+sumI)) ? maxsum : (sumR+sumI);
00046 else{
00047 maxsum = (sumR+sumI) ;
00048 init = 1;
00049 }
00050 ++ptr ;
00051 }
00052 return sqrt(maxsum);
00053 }
00054
00055 template <>
00056 Matrix<Complex> operator*(const double d,
00057 const Matrix<Complex> &a)
00058 {
00059 int i, size=a.rows()*a.cols() ;
00060 Matrix<Complex> b(a.rows(),a.cols());
00061
00062 Complex *bptr,*aptr ;
00063 bptr = b.m - 1 ;
00064 aptr = a.m - 1 ;
00065 for (i = size; i > 0 ; --i)
00066 *(++bptr) = (Complex)(d * (*(++aptr))) ;
00067 return b;
00068
00069 }
00070
00071 template <>
00072 Matrix<Complex> operator*(const Complex &d,
00073 const Matrix<Complex> &a)
00074 {
00075 int i, size=a.rows()*a.cols() ;
00076 Matrix<Complex> b(a.rows(),a.cols());
00077
00078 Complex *bptr,*aptr ;
00079 bptr = b.m - 1 ;
00080 aptr = a.m - 1 ;
00081 for (i = size; i > 0 ; --i)
00082 *(++bptr) = d * *(++aptr) ;
00083 return b;
00084 }
00085
00086
00087 template <>
00088 Matrix<Complex> Matrix<Complex>::herm() const
00089 {
00090 int i, j, r = cols(), c = rows();
00091 Matrix<Complex> adj(r,c);
00092
00093 for (i = 0; i < r; ++i)
00094 for (j = 0; j < c; ++j)
00095 adj.elem(i,j) = conj(elem(j,i)) ;
00096
00097 return adj;
00098 }
00099
00100
00101 #ifdef NO_IMPLICIT_TEMPLATES
00102
00103
00104
00105 template class Matrix<Complex> ;
00106
00107 template Matrix<Complex> operator+(const Matrix<Complex>&,const Matrix<Complex>&);
00108 template Matrix<Complex> operator-(const Matrix<Complex>&,const Matrix<Complex>&);
00109 template Matrix<Complex> operator*(const Matrix<Complex>&,const Matrix<Complex>&);
00110 template Vector<Complex> operator*(const Matrix<Complex>&,const Vector<Complex>&);
00111 template int operator==(const Matrix<Complex>&,const Matrix<Complex>&);
00112
00113 template Matrix<Complex> comm(const Matrix<Complex>&,const Matrix<Complex>&);
00114
00115
00116 #endif
00117
00118 }