diff --git a/mat.cc b/mat.cc index 10ec8c6..b18c6d8 100644 --- a/mat.cc +++ b/mat.cc @@ -6,6 +6,8 @@ //// forced instantization in the corresponding object file template NRMat; template NRMat< complex >; +template NRMat; +template NRMat; /* @@ -376,6 +378,21 @@ NRMat< complex >::operator*=(const complex &a) return *this; } +//and for general type +template +NRMat & NRMat::operator*=(const T &a) +{ + copyonwrite(); +#ifdef MATPTR + for (int i=0; i< nn*nn; i++) v[0][i] *= a; +#else + for (int i=0; i< nn*nn; i++) v[i] *= a; +#endif + return *this; +} + + + // Mat += Mat NRMat & NRMat::operator+=(const NRMat &rhs) { @@ -399,6 +416,24 @@ NRMat< complex >::operator+=(const NRMat< complex > &rhs) return *this; } +//and for general type +template +NRMat & NRMat::operator+=(const NRMat &rhs) +{ +#ifdef DEBUG + if (nn != rhs.nn || mm!= rhs.mm) + laerror("Mat -= Mat of incompatible matrices"); +#endif + copyonwrite(); +#ifdef MATPTR + for (int i=0; i< nn*nn; i++) v[0][i] += rhs.v[0][i] ; +#else + for (int i=0; i< nn*nn; i++) v[i] += rhs.v[i] ; +#endif + return *this; +} + + // Mat -= Mat NRMat & NRMat::operator-=(const NRMat &rhs) { @@ -422,6 +457,24 @@ NRMat< complex >::operator-=(const NRMat< complex > &rhs) return *this; } +//and for general type +template +NRMat & NRMat::operator-=(const NRMat &rhs) +{ +#ifdef DEBUG + if (nn != rhs.nn || mm!= rhs.mm) + laerror("Mat -= Mat of incompatible matrices"); +#endif + copyonwrite(); +#ifdef MATPTR + for (int i=0; i< nn*nn; i++) v[0][i] -= rhs.v[0][i] ; +#else + for (int i=0; i< nn*nn; i++) v[i] -= rhs.v[i] ; +#endif + return *this; +} + + // Mat += SMat NRMat & NRMat::operator+=(const NRSMat &rhs) { @@ -461,6 +514,28 @@ NRMat< complex >::operator+=(const NRSMat< complex > &rhs) return *this; } +//and for general type +template +NRMat & NRMat::operator+=(const NRSMat &rhs) +{ +#ifdef DEBUG + if (nn!=mm || nn!=rhs.nrows()) laerror("incompatible matrix size in Mat+=SMat"); +#endif + const T *p = rhs; + copyonwrite(); + for (int i=0; i & NRMat::operator-=(const NRSMat &rhs) { @@ -500,6 +575,28 @@ NRMat< complex >::operator-=(const NRSMat< complex > &rhs) return *this; } + +//and for general type +template +NRMat & NRMat::operator-=(const NRSMat &rhs) +{ +#ifdef DEBUG + if (nn!=mm || nn!=rhs.nrows()) laerror("incompatible matrix size in Mat+=SMat"); +#endif + const T *p = rhs; + copyonwrite(); + for (int i=0; i::dot(const NRMat &rhs) const { @@ -801,6 +898,8 @@ template istream & operator>>(istream &s, NRMat< T > &x); \ INSTANTIZE(double) INSTANTIZE(complex) +INSTANTIZE(int) +INSTANTIZE(char) export template diff --git a/nonclass.cc b/nonclass.cc index eae721c..532844a 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -15,6 +15,8 @@ template void lawritemat(FILE *file,const T *a,int r,int c,const char *form0, \ int nodim,int modulo, int issym); INSTANTIZE(double) INSTANTIZE(complex) +INSTANTIZE(int) +INSTANTIZE(char) template void lawritemat(FILE *file,const T *a,int r,int c,const char *form0, diff --git a/smat.cc b/smat.cc index 9880023..9588ae8 100644 --- a/smat.cc +++ b/smat.cc @@ -7,6 +7,8 @@ ////// forced instantization in the corresponding object file template NRSMat; template NRSMat< complex >; +template NRSMat; +template NRSMat; @@ -27,7 +29,7 @@ NRSMat::NRSMat(const NRMat &rhs) v = new T[NN2]; int i, j, k=0; for (i=0; i>(istream &s, NRSMat &x) return s; } +//not implemented yet +const NRVec NRSMat::operator*(NRVec const&rhs) const +{ +laerror("NRSMat::operator*(NRVec const&) not implemented yet"); +return rhs; +} + +const NRVec NRSMat::operator*(NRVec const&rhs) const +{ +laerror("NRSMat::operator*(NRVec const&) not implemented yet"); +return rhs; +} + + ////////////////////////////////////////////////////////////////////////////// //// forced instantization in the corespoding object file @@ -396,4 +412,6 @@ template istream & operator>>(istream &s, NRSMat< T > &x); \ INSTANTIZE(double) INSTANTIZE(complex) +INSTANTIZE(int) +INSTANTIZE(char) diff --git a/smat.h b/smat.h index 7a29511..c097128 100644 --- a/smat.h +++ b/smat.h @@ -120,9 +120,17 @@ inline NRSMat< complex > & NRSMat< complex >::operator*=(const complex & a) { copyonwrite(); - cblas_zscal(nn, (void *)(&a), (void *)v, 1); + cblas_zscal(NN2, (void *)(&a), (void *)v, 1); return *this; } +template +inline NRSMat & NRSMat::operator*=(const T & a) +{ + copyonwrite(); + for(int i=0; i >::operator+=(const NRSMat< complex > & rhs) return *this; } +template +inline NRSMat & NRSMat::operator+=(const NRSMat & rhs) +{ +#ifdef DEBUG + if (nn != rhs.nn) laerror("incompatible SMats in SMat::operator+="); +#endif + copyonwrite(); + for(int i=0; i & NRSMat::operator-=(const NRSMat & rhs) @@ -187,6 +207,18 @@ NRSMat< complex >::operator-=(const NRSMat< complex > & rhs) return *this; } +template +inline NRSMat & NRSMat::operator-=(const NRSMat & rhs) +{ +#ifdef DEBUG + if (nn != rhs.nn) laerror("incompatible SMats in SMat::operator-="); +#endif + copyonwrite(); + for(int i=0; i inline const NRMat NRSMat::operator+(const NRMat &rhs) const diff --git a/vec.cc b/vec.cc index 3187d25..23b5346 100644 --- a/vec.cc +++ b/vec.cc @@ -9,8 +9,12 @@ template istream & operator>>(istream &s, NRVec< T > &x); \ INSTANTIZE(double) INSTANTIZE(complex) +INSTANTIZE(int) +INSTANTIZE(char) template NRVec; template NRVec< complex >; +template NRVec; +template NRVec; /* @@ -281,6 +285,12 @@ NRVec< complex > & NRVec< complex >::normalize() return *this; } +//and for these types it does not make sense to normalize but we have them for linkage +NRVec & NRVec::normalize() {laerror("normalize() impossible for integer types"); return *this;} +NRVec & NRVec::normalize() {laerror("normalize() impossible for integer types"); return *this;} + + + // gemv call void NRVec::gemv(const double beta, const NRMat &A, const char trans, const double alpha, const NRVec &x) diff --git a/vec.h b/vec.h index 330ad4c..33fb761 100644 --- a/vec.h +++ b/vec.h @@ -160,6 +160,7 @@ inline NRVec & NRVec::operator+=(const double &a) cblas_daxpy(nn, 1.0, &a, 0, v, 1); return *this; } + inline NRVec< complex > & NRVec< complex >::operator+=(const complex &a) { @@ -168,6 +169,17 @@ NRVec< complex >::operator+=(const complex &a) return *this; } +//and for general type +template +inline NRVec & NRVec::operator+=(const T &a) +{ + copyonwrite(); + int i; + for(i=0; i & NRVec::operator-=(const double &a) { @@ -183,6 +195,17 @@ NRVec< complex >::operator-=(const complex &a) return *this; } +//and for general type +template +inline NRVec & NRVec::operator-=(const T &a) +{ + copyonwrite(); + int i; + for(i=0; i & NRVec::operator+=(const NRVec &rhs) { @@ -204,6 +227,20 @@ NRVec< complex >::operator+=(const NRVec< complex > &rhs) return *this; } +//and for general type +template +inline NRVec & NRVec::operator+=(const NRVec &rhs) +{ +#ifdef DEBUG + if (nn != rhs.nn) laerror("daxpy of incompatible vectors"); +#endif + copyonwrite(); + int i; + for(i=0; i & NRVec::operator-=(const NRVec &rhs) { @@ -225,6 +262,20 @@ NRVec< complex >::operator-=(const NRVec< complex > &rhs) return *this; } +//and for general type +template +inline NRVec & NRVec::operator-=(const NRVec &rhs) +{ +#ifdef DEBUG + if (nn != rhs.nn) laerror("daxpy of incompatible vectors"); +#endif + copyonwrite(); + int i; + for(i=0; i & NRVec::operator*=(const double &a) { @@ -240,6 +291,17 @@ NRVec< complex >::operator*=(const complex &a) return *this; } +//and for general type +template +inline NRVec & NRVec::operator*=(const T &a) +{ + copyonwrite(); + int i; + for(i=0; i::operator*(const NRVec &rhs) const {