diff --git a/bitvector.cc b/bitvector.cc index 6f6ed57..c79a465 100644 --- a/bitvector.cc +++ b/bitvector.cc @@ -17,6 +17,7 @@ */ #include "bitvector.h" +#include namespace LA { @@ -187,5 +188,26 @@ if(modulo) return s+word_popul(a); } +void bitvector::read(int fd, bool dimensions, bool transp) +{ +if(dimensions) + { + int r = ::read(fd,&modulo,sizeof(modulo)); + if(r!=sizeof(modulo)) laerror("cannot read in bitvector"); + } +NRVec::get(fd,dimensions,transp); +} + +void bitvector::write(int fd, bool dimensions, bool transp) +{ +if(dimensions) + { + int r = ::write(fd,&modulo,sizeof(modulo)); + if(r!=sizeof(modulo)) laerror("cannot write in bitvector"); + } +NRVec::put(fd,dimensions,transp); +} + + }//namespace diff --git a/bitvector.h b/bitvector.h index bbc6e54..1bb14d9 100644 --- a/bitvector.h +++ b/bitvector.h @@ -75,6 +75,9 @@ public: //extended, truncated const i.e. not on *this but return new entity, take care of modulo's bits //logical shifts <<= >>= << >> not implemented yet //logical rotations not implemented yet + //unformatted file IO + void read(int fd, bool dimensions=1, bool transp=0); + void write(int fd, bool dimensions=1, bool transp=0); }; extern std::ostream & operator<<(std::ostream &s, const bitvector &x); diff --git a/csrmat.cc b/csrmat.cc index 4f319ce..621f942 100644 --- a/csrmat.cc +++ b/csrmat.cc @@ -48,7 +48,7 @@ template void CSRMat::put(int fd, bool dimen, bool transp) const; \ INSTANTIZE(double) -INSTANTIZE(complex) +INSTANTIZE(std::complex) */ //// forced instantization of functions in the header in the corresponding object file diff --git a/la.h b/la.h index 62d179f..e15c857 100644 --- a/la.h +++ b/la.h @@ -48,9 +48,9 @@ using namespace LA; typedef NRMat NRIMat; typedef NRMat NRDMat; -typedef NRMat > NRCMat; +typedef NRMat > NRCMat; typedef NRVec NRIVec; typedef NRVec NRDVec; -typedef NRVec > NRCVec; +typedef NRVec > NRCVec; #endif /* _LA_H_ */ diff --git a/la_traits.h b/la_traits.h index 6b5c9cb..2687be7 100644 --- a/la_traits.h +++ b/la_traits.h @@ -39,7 +39,6 @@ //using namespace std; -#define complex std::complex #include "laerror.h" @@ -79,6 +78,18 @@ template class SparseMat; template class SparseSMat; template class CSRMat; +//trick to allow real and imag part of complex as l-values +template +T &real(std::complex &c) { + return reinterpret_cast(&c)[0]; +} +template +T &imag(std::complex &c) { + return reinterpret_cast(&c)[1]; +} + + +// typedef class {} Dummy_type; typedef class {} Dummy_type2; @@ -96,7 +107,7 @@ struct LA_traits_complex #define SPECIALIZE_COMPLEX(T) \ template<> \ -struct LA_traits_complex > \ +struct LA_traits_complex > \ { \ typedef T Component_type; \ typedef NRVec NRVec_Noncomplex_type; \ @@ -106,9 +117,9 @@ struct LA_traits_complex > \ SPECIALIZE_COMPLEX(double) -SPECIALIZE_COMPLEX(complex) +SPECIALIZE_COMPLEX(std::complex) SPECIALIZE_COMPLEX(float) -SPECIALIZE_COMPLEX(complex) +SPECIALIZE_COMPLEX(std::complex) SPECIALIZE_COMPLEX(char) SPECIALIZE_COMPLEX(unsigned char) SPECIALIZE_COMPLEX(short) @@ -172,9 +183,9 @@ class isscalar { public: typedef scalar_false scalar_type;}; template<>\ class isscalar {public: typedef scalar_true scalar_type;};\ template<>\ -class isscalar > {public: typedef scalar_true scalar_type;};\ +class isscalar > {public: typedef scalar_true scalar_type;};\ template<>\ -class isscalar > > {public: typedef scalar_true scalar_type;};\ +class isscalar > > {public: typedef scalar_true scalar_type;};\ //declare what is scalar @@ -211,56 +222,56 @@ template struct LA_traits_aux //complex scalars template -struct LA_traits_aux, scalar_true> { -typedef complex elementtype; -typedef complex producttype; +struct LA_traits_aux, scalar_true> { +typedef std::complex elementtype; +typedef std::complex producttype; typedef C normtype; typedef C realtype; -typedef complex complextype; -static inline C sqrabs(const complex x) { return x.real()*x.real()+x.imag()*x.imag();} -static inline bool gencmp(const complex *x, const complex *y, size_t n) {return memcmp(x,y,n*sizeof(complex));} -static bool bigger(const complex &x, const complex &y) {laerror("complex comparison undefined"); return false;} -static bool smaller(const complex &x, const complex &y) {laerror("complex comparison undefined"); return false;} -static inline normtype norm (const complex &x) {return std::abs(x);} -static inline void axpy (complex &s, const complex &x, const complex &c) {s+=x*c;} -static inline void get(int fd, complex &x, bool dimensions=0, bool transp=0) {if(sizeof(complex)!=read(fd,&x,sizeof(complex))) laerror("read error");} -static inline void put(int fd, const complex &x, bool dimensions=0, bool transp=0) {if(sizeof(complex)!=write(fd,&x,sizeof(complex))) laerror("write error");} -static void multiget(size_t n,int fd, complex *x, bool dimensions=0) +typedef std::complex complextype; +static inline C sqrabs(const std::complex x) { return x.real()*x.real()+x.imag()*x.imag();} +static inline bool gencmp(const std::complex *x, const std::complex *y, size_t n) {return memcmp(x,y,n*sizeof(std::complex));} +static bool bigger(const std::complex &x, const std::complex &y) {laerror("std::complex comparison undefined"); return false;} +static bool smaller(const std::complex &x, const std::complex &y) {laerror("std::complex comparison undefined"); return false;} +static inline normtype norm (const std::complex &x) {return std::abs(x);} +static inline void axpy (std::complex &s, const std::complex &x, const std::complex &c) {s+=x*c;} +static inline void get(int fd, std::complex &x, bool dimensions=0, bool transp=0) {if(sizeof(std::complex)!=read(fd,&x,sizeof(std::complex))) laerror("read error");} +static inline void put(int fd, const std::complex &x, bool dimensions=0, bool transp=0) {if(sizeof(std::complex)!=write(fd,&x,sizeof(std::complex))) laerror("write error");} +static void multiget(size_t n,int fd, std::complex *x, bool dimensions=0) { size_t total=0; - size_t system_limit = (1L<<30)/sizeof(complex); //do not expect too much from the system and read at most 1GB at once + size_t system_limit = (1L<<30)/sizeof(std::complex); //do not expect too much from the system and read at most 1GB at once ssize_t r; size_t nn; do{ - r=read(fd,x+total,nn=(n-total > system_limit ? system_limit : n-total)*sizeof(complex)); + r=read(fd,x+total,nn=(n-total > system_limit ? system_limit : n-total)*sizeof(std::complex)); if(r<0 || r==0 && nn!=0 ) {std::cout<<"read returned "<); - if(r%sizeof(complex)) laerror("read error 2"); + else total += r/sizeof(std::complex); + if(r%sizeof(std::complex)) laerror("read error 2"); } while(total < n); } -static void multiput(size_t n, int fd, const complex *x, bool dimensions=0) +static void multiput(size_t n, int fd, const std::complex *x, bool dimensions=0) { size_t total=0; - size_t system_limit = (1L<<30)/sizeof(complex); //do not expect too much from the system and write at most 1GB at once + size_t system_limit = (1L<<30)/sizeof(std::complex); //do not expect too much from the system and write at most 1GB at once ssize_t r; size_t nn; do{ - r=write(fd,x+total,nn=(n-total > system_limit ? system_limit : n-total)*sizeof(complex)); + r=write(fd,x+total,nn=(n-total > system_limit ? system_limit : n-total)*sizeof(std::complex)); if(r<0 || r==0 && nn!=0 ) {std::cout<<"write returned "<); - if(r%sizeof(complex)) laerror("write error 2"); + else total += r/sizeof(std::complex); + if(r%sizeof(std::complex)) laerror("write error 2"); } while(total < n); } -static void copy(complex *dest, complex *src, size_t n) {memcpy(dest,src,n*sizeof(complex));} -static void clear(complex *dest, size_t n) {memset(dest,0,n*sizeof(complex));} -static void copyonwrite(complex &x) {}; -static void clearme(complex &x) {x=0;}; -static void deallocate(complex &x) {}; -static inline complex conjugate(const complex &x) {return complex(x.real(),-x.imag());}; -static inline C realpart(const complex &x) {return x.real();} -static inline C imagpart(const complex &x) {return x.imag();} +static void copy(std::complex *dest, std::complex *src, size_t n) {memcpy(dest,src,n*sizeof(std::complex));} +static void clear(std::complex *dest, size_t n) {memset(dest,0,n*sizeof(std::complex));} +static void copyonwrite(std::complex &x) {}; +static void clearme(std::complex &x) {x=0;}; +static void deallocate(std::complex &x) {}; +static inline std::complex conjugate(const std::complex &x) {return std::complex(x.real(),-x.imag());}; +static inline C realpart(const std::complex &x) {return x.real();} +static inline C imagpart(const std::complex &x) {return x.imag();} }; @@ -271,7 +282,7 @@ typedef C elementtype; typedef C producttype; typedef C normtype; typedef C realtype; -typedef complex complextype; +typedef std::complex complextype; static inline C sqrabs(const C x) { return x*x;} static inline bool gencmp(const C *x, const C *y, size_t n) {return memcmp(x,y,n*sizeof(C));} static inline bool bigger(const C &x, const C &y) {return x>y;} diff --git a/mat.cc b/mat.cc index b7fdfa9..350ef1c 100644 --- a/mat.cc +++ b/mat.cc @@ -276,7 +276,7 @@ NRMat& NRMat::operator=(const double &a){ * @return reference to the modified matrix ******************************************************************************/ template <> -NRMat >& NRMat >::operator=(const complex &a){ +NRMat >& NRMat >::operator=(const std::complex &a){ const int n2 = nn*nn; copyonwrite(); #ifdef DEBUG @@ -286,7 +286,7 @@ NRMat >& NRMat >::operator=(const complex)); + memset(v[0], 0, n2*sizeof(std::complex)); for(register int i=0; i< nn; i++) v[i][i] = a; #else //set all matrix elements equal to zero @@ -340,7 +340,7 @@ NRMat & NRMat::operator+=(const double& a) { * @return reference to the modified matrix ******************************************************************************/ template <> -NRMat > & NRMat >::operator+=(const complex& a) { +NRMat > & NRMat >::operator+=(const std::complex& a) { copyonwrite(); #ifdef DEBUG if(nn != mm) laerror("nonsquare matrix"); @@ -355,7 +355,7 @@ NRMat > & NRMat >::operator+=(const complex* d = gpuputcomplex(a); + std::complex* d = gpuputcomplex(a); cublasZaxpy(nn, CUMONE, (cuDoubleComplex*)d, 0, (cuDoubleComplex*)v, nn+1); TEST_CUBLAS("cublasDaxpy"); gpufree(d); @@ -402,7 +402,7 @@ NRMat& NRMat::operator-=(const double& a) { * @return reference to the modified matrix ******************************************************************************/ template <> -NRMat >& NRMat >::operator-=(const complex& a) { +NRMat >& NRMat >::operator-=(const std::complex& a) { copyonwrite(); #ifdef DEBUG if(nn != mm) laerror("nonsquare matrix"); @@ -417,7 +417,7 @@ NRMat >& NRMat >::operator-=(const complex* d = gpuputcomplex(a); + std::complex* d = gpuputcomplex(a); cublasZaxpy(nn, CUMONE, (cuDoubleComplex*)d, 0, (cuDoubleComplex*)v, nn+1); TEST_CUBLAS("cublasDaxpy"); gpufree(d); @@ -502,16 +502,16 @@ const NRMat NRMat::operator-() const { * @return modified copy of this matrix ******************************************************************************/ template <> -const NRMat > NRMat >::operator-() const { +const NRMat > NRMat >::operator-() const { const size_t nm = (size_t)nn*mm; - NRMat > result(nn, mm, getlocation()); + NRMat > result(nn, mm, getlocation()); #ifdef CUDALA if(location == cpu) { #endif #ifdef MATPTR for(register size_t i=0; i)); + memcpy(result.v, v, nm*sizeof(std::complex)); cblas_zscal(nm, &CMONE, result.v, 1); #endif #ifdef CUDALA @@ -631,8 +631,8 @@ const NRVec NRMat::csum() const { * @return summed columns in a form of a vector ******************************************************************************/ template <> -const NRVec > NRMat >::csum() const { - NRVec > result(nn, getlocation()); +const NRVec > NRMat >::csum() const { + NRVec > result(nn, getlocation()); result = 0.0; #ifdef CUDALA if(location == cpu){ @@ -699,8 +699,8 @@ const NRVec NRMat::rsum() const { * @return summed rows in a form of a vector ******************************************************************************/ template <> -const NRVec > NRMat >::rsum() const { - NRVec > result(mm, getlocation()); +const NRVec > NRMat >::rsum() const { + NRVec > result(mm, getlocation()); result = 0.0; #ifdef CUDALA if(location == cpu){ @@ -873,24 +873,24 @@ laerror("nonsymmetry not implemented on GPU yet"); * compute matrix non-hermiticity ******************************************************************************/ template <> -const double NRMat >::nonhermiticity() const { +const double NRMat >::nonhermiticity() const { #ifdef DEBUG if (nn != mm) laerror("NRMat:nonsymmetry() invalid for non-square matrix"); #endif double sum = 0; -complex tmp; +std::complex tmp; #ifdef CUDALA if(location == cpu){ #endif for(register int i=1; i (v[i][j].real()-v[j][i].real(),v[i][j].imag()+v[j][i].imag()); + tmp = std::complex (v[i][j].real()-v[j][i].real(),v[i][j].imag()+v[j][i].imag()); #else register int a, b; a = i*(size_t)mm + j; b = j*(size_t)mm + i; - tmp = complex (v[a].real() - v[b].real(), v[a].imag()+v[b].imag()); + tmp = std::complex (v[a].real() - v[b].real(), v[a].imag()+v[b].imag()); #endif sum += tmp.real()*tmp.real()+tmp.imag()*tmp.imag(); } @@ -905,12 +905,12 @@ laerror("nonsymmetry not implemented on GPU yet"); } template <> -const double NRMat >::nonsymmetry() const { +const double NRMat >::nonsymmetry() const { #ifdef DEBUG if (nn != mm) laerror("NRMat:nonsymmetry() invalid for non-square matrix"); #endif double sum = 0; -complex tmp; +std::complex tmp; #ifdef CUDALA if(location == cpu){ #endif @@ -945,28 +945,28 @@ laerror("nonsymmetry not implemented on GPU yet"); * or imaginary part of the complex matrix being created ******************************************************************************/ template<> -NRMat >::NRMat(const NRMat &rhs, bool imagpart): nn(rhs.nrows()), mm(rhs.ncols()), count(new int(1)) { +NRMat >::NRMat(const NRMat &rhs, bool imagpart): nn(rhs.nrows()), mm(rhs.ncols()), count(new int(1)) { const size_t nn_mm = (size_t)nn*mm; #ifdef CUDALA if(location == cpu){ #endif #ifdef MATPTR - v = new complex*[n]; - v[0] = new complex[nn_mm]; + v = new std::complex*[n]; + v[0] = new std::complex[nn_mm]; for(register int i=1; i)); + memset(v[0], 0, nn_mm*sizeof(std::complex)); cblas_dcopy(nn_mm, &rhs[0][0], 1, ((double *)v[0]) + (imagpart?1:0), 2); #else - v = new complex[nn_mm]; - memset(v, 0, nn_mm*sizeof(complex)); + v = new std::complex[nn_mm]; + memset(v, 0, nn_mm*sizeof(std::complex)); cblas_dcopy(nn_mm, &rhs[0][0], 1, ((double *)v) + (imagpart?1:0), 2); #endif #ifdef CUDALA }else{ - v = (complex*)gpualloc(sizeof(complex)*nn_mm); - complex *_val = gpuputcomplex(CZERO); + v = (std::complex*)gpualloc(sizeof(std::complex)*nn_mm); + std::complex *_val = gpuputcomplex(CZERO); cublasZcopy(nn_mm, (cuDoubleComplex*)_val, 0, (cuDoubleComplex*)v, 1); TEST_CUBLAS("cublasZcopy"); gpufree(_val); @@ -986,7 +986,7 @@ NRMat >::NRMat(const NRMat &rhs, bool imagpart): nn(rhs. * or imaginary part of the input complex matrix ******************************************************************************/ template<> -NRMat::NRMat(const NRMat > &rhs, bool imagpart): nn(rhs.nrows()), mm(rhs.ncols()), count(new int(1)) { +NRMat::NRMat(const NRMat > &rhs, bool imagpart): nn(rhs.nrows()), mm(rhs.ncols()), count(new int(1)) { const size_t nn_mm = (size_t) nn*mm; #ifdef CUDALA if(location == cpu){ @@ -1109,9 +1109,9 @@ const NRSMat NRMat::transposedtimes() const { * @return complex NRSMat object because of the hermiticity of \f$A^\dagger{}A\f$ ******************************************************************************/ template<> -const NRSMat > NRMat >::transposedtimes() const { +const NRSMat > NRMat >::transposedtimes() const { int i(0), j(0); - NRSMat > r(mm, getlocation()); + NRSMat > r(mm, getlocation()); #ifdef CUDALA if(location == cpu){ #endif @@ -1130,7 +1130,7 @@ const NRSMat > NRMat >::transposedtimes() const for(j=0; j<=i; ++j){ cuDoubleComplex val = cublasZdotc(nn, (const cuDoubleComplex*)(v + i), mm, (const cuDoubleComplex*)(v + j), mm); TEST_CUBLAS("cublasZdotc"); - r(i, j) = *(reinterpret_cast*> (&val)); + r(i, j) = *(reinterpret_cast*> (&val)); } } r.moveto(this->location); @@ -1201,9 +1201,9 @@ const NRSMat NRMat::timestransposed() const { * @return complex NRSMat object because of the hermiticity of \f$AA^\dagger{}\f$ ******************************************************************************/ template<> -const NRSMat > NRMat >::timestransposed() const { +const NRSMat > NRMat >::timestransposed() const { int i(0), j(0); - NRSMat > r(nn, getlocation()); + NRSMat > r(nn, getlocation()); #ifdef CUDALA if(location == cpu){ #endif @@ -1222,7 +1222,7 @@ const NRSMat > NRMat >::timestransposed() const for(j=0; j<=i; ++j){ cuDoubleComplex val = cublasZdotc(nn, (const cuDoubleComplex *)(v + i*(size_t)mm), 1, (const cuDoubleComplex *)(v + j*(size_t)mm), 1); TEST_CUBLAS("cublasZdotc"); - r(i, j) = *(reinterpret_cast*> (&val)); + r(i, j) = *(reinterpret_cast*> (&val)); } } r.moveto(this->location); @@ -1287,7 +1287,7 @@ void NRMat::randomize(const double &x) { * @param[in] x generate random numbers from the interval [0, x] ******************************************************************************/ template<> -void NRMat >::randomize(const double &x) { +void NRMat >::randomize(const double &x) { #ifdef CUDALA if(location == cpu){ #endif @@ -1295,17 +1295,17 @@ void NRMat >::randomize(const double &x) { for(register int j=0; j(re, im); + (*this)(i,j) = std::complex(re, im); } } #ifdef CUDALA }else{ - NRMat > tmp(nn, mm, cpu); - complex *tmp_data = tmp; + NRMat > tmp(nn, mm, cpu); + std::complex *tmp_data = tmp; for(register size_t i=0; i<(size_t)nn*mm; ++i){ const double re = x*(2.*random()/(1. + RAND_MAX) - 1.); const double im = x*(2.*random()/(1. + RAND_MAX) - 1.); - tmp_data[i] = complex(re, im); + tmp_data[i] = std::complex(re, im); } tmp.moveto(this->location); *this |= tmp; @@ -1342,8 +1342,8 @@ NRMat& NRMat::operator*=(const double &a) { * @return reference to the modified matrix ******************************************************************************/ template<> -NRMat > & -NRMat >::operator*=(const complex &a) { +NRMat > & +NRMat >::operator*=(const std::complex &a) { copyonwrite(); #ifdef CUDALA if(location == cpu){ @@ -1409,8 +1409,8 @@ NRMat & NRMat::operator+=(const NRMat &rhs) { * @return reference to the modified matrix ******************************************************************************/ template<> -NRMat > & -NRMat >::operator+=(const NRMat< complex > &rhs) { +NRMat > & +NRMat >::operator+=(const NRMat< std::complex > &rhs) { #ifdef DEBUG if (nn != rhs.nn || mm != rhs.mm) laerror("incompatible matrices"); #endif @@ -1483,8 +1483,8 @@ NRMat & NRMat::operator-=(const NRMat &rhs) { * @return reference to the modified matrix ******************************************************************************/ template<> -NRMat< complex > & -NRMat< complex >::operator-=(const NRMat< complex > &rhs) { +NRMat< std::complex > & +NRMat< std::complex >::operator-=(const NRMat< std::complex > &rhs) { #ifdef DEBUG if (nn != rhs.nn || mm != rhs.mm) laerror("incompatible matrices"); #endif @@ -1584,13 +1584,13 @@ NRMat & NRMat::operator+=(const NRSMat &rhs) { * @see NRSMat ******************************************************************************/ template<> -NRMat< complex > & -NRMat< complex >::operator+=(const NRSMat< complex > &rhs) +NRMat< std::complex > & +NRMat< std::complex >::operator+=(const NRSMat< std::complex > &rhs) { #ifdef DEBUG if (nn != rhs.nn || mm != rhs.nn) laerror("incompatible matrices"); #endif - const complex *p = rhs; + const std::complex *p = rhs; SAME_LOC(*this, rhs); copyonwrite(); @@ -1709,12 +1709,12 @@ NRMat & NRMat::operator-=(const NRSMat &rhs) * @see NRSMat ******************************************************************************/ template<> -NRMat > & -NRMat >::operator-=(const NRSMat > &rhs) { +NRMat > & +NRMat >::operator-=(const NRSMat > &rhs) { #ifdef DEBUG if (nn != rhs.nn || mm != rhs.nn) laerror("incompatible matrices"); #endif - const complex *p = rhs; + const std::complex *p = rhs; SAME_LOC(*this, rhs); copyonwrite(); @@ -1810,13 +1810,13 @@ const double NRMat::dot(const NRMat &rhs) const { * @return computed scalar product ******************************************************************************/ template<> -const complex -NRMat >::dot(const NRMat > &rhs) const { +const std::complex +NRMat >::dot(const NRMat > &rhs) const { #ifdef DEBUG - if(nn != rhs.nn || mm != rhs.mm) laerror("incompatible matrices in NRMat >::dot(const NRMat >&)"); + if(nn != rhs.nn || mm != rhs.mm) laerror("incompatible matrices in NRMat >::dot(const NRMat >&)"); #endif - complex ret(0.0, 0.0); + std::complex ret(0.0, 0.0); #ifdef CUDALA if(location == cpu){ #endif @@ -1824,7 +1824,7 @@ NRMat >::dot(const NRMat > &rhs) const { #ifdef CUDALA }else{ cuDoubleComplex val = cublasZdotc((size_t)nn*mm, (cuDoubleComplex*)v, 1, (cuDoubleComplex*)(rhs.v), 1); - ret = *(reinterpret_cast*> (&val)); + ret = *(reinterpret_cast*> (&val)); } #endif return ret; @@ -1862,14 +1862,14 @@ const NRMat NRMat::operator*(const NRMat &rhs) const { * @return computed product by value ******************************************************************************/ template<> -const NRMat< complex > -NRMat< complex >::operator*(const NRMat< complex > &rhs) const { +const NRMat< std::complex > +NRMat< std::complex >::operator*(const NRMat< std::complex > &rhs) const { #ifdef DEBUG - if(mm != rhs.nn) laerror("incompatible matrices in NRMat >::operator*(const NRMat >&)"); + if(mm != rhs.nn) laerror("incompatible matrices in NRMat >::operator*(const NRMat >&)"); if(rhs.mm <= 0) laerror("illegal matrix dimension in gemm"); #endif SAME_LOC(*this, rhs); - NRMat > result(nn, rhs.mm, getlocation()); + NRMat > result(nn, rhs.mm, getlocation()); #ifdef CUDALA if(location == cpu){ #endif @@ -1914,9 +1914,9 @@ void NRMat::diagmultl(const NRVec &rhs) { * @param[in] rhs complex vector represeting the diagonal of matrix \f$D\f$ ******************************************************************************/ template<> -void NRMat< complex >::diagmultl(const NRVec< complex > &rhs) { +void NRMat< std::complex >::diagmultl(const NRVec< std::complex > &rhs) { #ifdef DEBUG - if (nn != rhs.size()) laerror("incompatible matrices in NRMat >::diagmultl(const NRVec >&)"); + if (nn != rhs.size()) laerror("incompatible matrices in NRMat >::diagmultl(const NRVec >&)"); #endif NOT_GPU(rhs); copyonwrite(); @@ -1965,9 +1965,9 @@ void NRMat::diagmultr(const NRVec &rhs) { * @param[in] rhs complex vector represeting the diagonal of matrix \f$D\f$ ******************************************************************************/ template<> -void NRMat< complex >::diagmultr(const NRVec< complex > &rhs) { +void NRMat< std::complex >::diagmultr(const NRVec< std::complex > &rhs) { #ifdef DEBUG - if(mm != rhs.size()) laerror("incompatible matrices in NRMat >::diagmultr(const NRVec >&)"); + if(mm != rhs.size()) laerror("incompatible matrices in NRMat >::diagmultr(const NRVec >&)"); #endif NOT_GPU(rhs); copyonwrite(); @@ -2027,14 +2027,14 @@ NRMat::operator*(const NRSMat &rhs) const { * @return \f$A\times\S\f$ by value ******************************************************************************/ template<> -const NRMat< complex > -NRMat< complex >::operator*(const NRSMat< complex > &rhs) const { +const NRMat< std::complex > +NRMat< std::complex >::operator*(const NRSMat< std::complex > &rhs) const { #ifdef DEBUG - if(mm != rhs.nrows()) laerror("incompatible matrices int NRMat >::operator*(const NRSMat > &)"); + if(mm != rhs.nrows()) laerror("incompatible matrices int NRMat >::operator*(const NRSMat > &)"); #endif SAME_LOC(*this, rhs); const int rhs_ncols = rhs.ncols(); - NRMat > result(nn, rhs_ncols, getlocation()); + NRMat > result(nn, rhs_ncols, getlocation()); #ifdef CUDALA if(location == cpu){ @@ -2068,7 +2068,7 @@ NRMat& NRMat::conjugateme() { * @return reference to the modified matrix ******************************************************************************/ template<> -NRMat >& NRMat >::conjugateme() { +NRMat >& NRMat >::conjugateme() { copyonwrite(); #ifdef CUDALA if(location == cpu){ @@ -2111,9 +2111,9 @@ const NRMat NRMat::transpose(bool conj) const { * @return transposed (conjugated) matrix by value ******************************************************************************/ template<> -const NRMat > -NRMat >::transpose(bool conj) const { - NRMat > result(mm, nn, getlocation()); +const NRMat > +NRMat >::transpose(bool conj) const { + NRMat > result(mm, nn, getlocation()); #ifdef CUDALA if(location == cpu){ #endif @@ -2177,10 +2177,10 @@ void NRMat::gemm(const double &beta, const NRMat &a, template<> -void NRMat >::gemm(const complex & beta, - const NRMat > & a, const char transa, - const NRMat > & b, const char transb, - const complex & alpha) +void NRMat >::gemm(const std::complex & beta, + const NRMat > & a, const char transa, + const NRMat > & b, const char transb, + const std::complex & alpha) { int k(tolower(transa)=='n'?a.mm:a.nn); @@ -2188,7 +2188,7 @@ void NRMat >::gemm(const complex & beta, int l(tolower(transa)=='n'?a.nn:a.mm); int kk(tolower(transb)=='n'?b.nn:b.mm); int ll(tolower(transb)=='n'?b.mm:b.nn); - if (l!=nn || ll!=mm || k!=kk) laerror("incompatible matrices in NRMat >::gemm(...)"); + if (l!=nn || ll!=mm || k!=kk) laerror("incompatible matrices in NRMat >::gemm(...)"); #endif if (alpha==CZERO && beta==CONE) return; @@ -2246,7 +2246,7 @@ const double NRMat::norm(const double scalar) const { * @return computed norm ******************************************************************************/ template<> -const double NRMat >::norm(const complex scalar) const { +const double NRMat >::norm(const std::complex scalar) const { if(scalar == CZERO){ #ifdef CUDALA if(location == cpu){ @@ -2263,7 +2263,7 @@ const double NRMat >::norm(const complex scalar) const { double sum(0.0); for(register int i=0; i tmp(0.0, 0.0); + register std::complex tmp(0.0, 0.0); #ifdef MATPTR tmp = v[i][j]; #else @@ -2310,10 +2310,10 @@ void NRMat::axpy(const double alpha, const NRMat &mat) { * @param[in] mat complex matrix \f$B\f$ ******************************************************************************/ template<> -void NRMat >::axpy(const complex alpha, - const NRMat > & mat) { +void NRMat >::axpy(const std::complex alpha, + const NRMat > & mat) { #ifdef DEBUG - if (nn != mat.nn || mm != mat.mm) laerror("incompatible matrices in NRMat >::axpy(...)"); + if (nn != mat.nn || mm != mat.mm) laerror("incompatible matrices in NRMat >::axpy(...)"); #endif SAME_LOC(*this, mat); copyonwrite(); @@ -2470,10 +2470,10 @@ int nnmin= nn<=mm?nn:mm; * @return void ******************************************************************************/ template<> -void NRMat >::diagonalset(const NRVec > &r) { +void NRMat >::diagonalset(const NRVec > &r) { int nnmin= nn<=mm?nn:mm; #ifdef DEBUG - if(r.size() != nnmin) laerror("incompatible vectors int NRMat >::diagonalset(...)"); + if(r.size() != nnmin) laerror("incompatible vectors int NRMat >::diagonalset(...)"); #endif SAME_LOC(*this, r); copyonwrite(); @@ -2683,7 +2683,7 @@ NRMat& NRMat::swap_rows(const int i, const int j){ * @return reference to the modified matrix ******************************************************************************/ template<> -NRMat >& NRMat >::swap_rows(){ +NRMat >& NRMat >::swap_rows(){ copyonwrite(); const int n_pul = this->nn >> 1; @@ -2705,7 +2705,7 @@ NRMat >& NRMat >::swap_rows(){ } template<> -NRMat >& NRMat >::swap_rows(const int i, const int j){ +NRMat >& NRMat >::swap_rows(const int i, const int j){ copyonwrite(); #ifdef CUDALA @@ -2802,7 +2802,7 @@ NRMat& NRMat::swap_cols(const int i, const int j){ * @return reference to the modified matrix ******************************************************************************/ template<> -NRMat >& NRMat >::swap_cols(){ +NRMat >& NRMat >::swap_cols(){ copyonwrite(); const int m_pul = mm >> 1; @@ -2824,7 +2824,7 @@ NRMat >& NRMat >::swap_cols(){ } template<> -NRMat >& NRMat >::swap_cols(const int i, const int j){ +NRMat >& NRMat >::swap_cols(const int i, const int j){ copyonwrite(); #ifdef CUDALA @@ -2986,11 +2986,11 @@ NRMat& NRMat::swap_rows_cols(){ * @return reference to the modified matrix ******************************************************************************/ template<> -NRMat >& NRMat >::swap_rows_cols(){ +NRMat >& NRMat >::swap_rows_cols(){ const int n_pul = nn >> 1; const int m_pul = mm >> 1; - complex tmp(0.0, 0.0); + std::complex tmp(0.0, 0.0); copyonwrite(); #ifdef CUDALA @@ -3014,7 +3014,7 @@ NRMat >& NRMat >::swap_rows_cols(){ TEST_CUBLAS("cublasZswap"); } if(nn & 1){ - void *gpu_ptr = gpualloc(sizeof(complex)*mm); + void *gpu_ptr = gpualloc(sizeof(std::complex)*mm); cublasZswap(mm, (cuDoubleComplex*)(v + n_pul*mm + mm - 1), -1, (cuDoubleComplex*)gpu_ptr, 1); cublasZcopy(mm, (cuDoubleComplex*)gpu_ptr, 1, (cuDoubleComplex*)(v + n_pul*mm), 1); gpufree(gpu_ptr); @@ -3079,7 +3079,7 @@ NRMat& NRMat::swap_rows_cols(){ * forced instantization in the corresponding object file ******************************************************************************/ template class NRMat; -template class NRMat >; +template class NRMat >; template class NRMat; template class NRMat; template class NRMat; diff --git a/mat.h b/mat.h index 9fd9fb1..0fd28fe 100644 --- a/mat.h +++ b/mat.h @@ -83,7 +83,7 @@ public: //! complexifying constructor NRMat(const typename LA_traits_complex::NRMat_Noncomplex_type &rhs, bool imagpart = false); //! explicit decomplexifying constructor - explicit NRMat(const NRMat > &rhs, bool imagpart = false); + explicit NRMat(const NRMat > &rhs, bool imagpart = false); //! explicit constructor converting symmetric matrix stored in packed form into a NRMat object explicit NRMat(const NRSMat &rhs); @@ -191,8 +191,8 @@ public: return result; }; //! multiply this matrix of general type T by vector of type complex - const NRVec > operator*(const NRVec > &rhs) const { - NRVec > result(nn, rhs.getlocation()); + const NRVec > operator*(const NRVec > &rhs) const { + NRVec > result(nn, rhs.getlocation()); result.gemv((T)0, *this, 'n', (T)1, rhs); return result; }; @@ -243,7 +243,7 @@ public: //! perform the gemv operation with vector of type T void gemv(const T beta, NRVec &r, const char trans, const T alpha, const NRVec &x) const { r.gemv(beta, *this, trans, alpha, x); }; //! perform the gemv operation with vector of type complex - void gemv(const T beta, NRVec > &r, const char trans, const T alpha, const NRVec > &x) const { r.gemv(beta, *this, trans, alpha, x); }; + void gemv(const T beta, NRVec > &r, const char trans, const T alpha, const NRVec > &x) const { r.gemv(beta, *this, trans, alpha, x); }; //! determine the pointer to the ith row inline T* operator[](const int i); @@ -289,6 +289,7 @@ public: //! get the pointer to the data inline operator T*(); + //! get the const pointer to the data inline operator const T*() const; @@ -858,7 +859,7 @@ inline const double NRMat::amin() const{ * @return \f$A_{l,m}\f$ which maximizes \f$\left\{\left|\Re{}A_{i,j}\right|+\left|\Im{}A_{i,j}\right|\right}\f$ ******************************************************************************/ template<> -inline const complex NRMat >::amax() const{ +inline const std::complex NRMat >::amax() const{ #ifdef CUDALA if(location == cpu){ #endif @@ -869,10 +870,10 @@ inline const complex NRMat >::amax() const{ #endif #ifdef CUDALA }else{ - complex ret(0.0, 0.0); + std::complex ret(0.0, 0.0); const size_t pozice = cublasIzamax((size_t)nn*mm, (cuDoubleComplex*)v, 1) - 1; TEST_CUBLAS("cublasIzamax"); - gpuget(1, sizeof(complex), v + pozice, &ret); + gpuget(1, sizeof(std::complex), v + pozice, &ret); return ret; } #endif @@ -885,8 +886,8 @@ inline const complex NRMat >::amax() const{ * @return \f$A_{l,m}\f$ which minimizes \f$\left\{\left|\Re{}A_{i,j}\right|+\left|\Im{}A_{i,j}\right|\right}\f$ ******************************************************************************/ template<> -inline const complex NRMat >::amin() const{ - complex ret(0.0, 0.0); +inline const std::complex NRMat >::amin() const{ + std::complex ret(0.0, 0.0); #ifdef CUDALA if(location == cpu){ #endif @@ -894,7 +895,7 @@ inline const complex NRMat >::amin() const{ const size_t nm = (size_t)nn*mm; int index(-1); double val(0.0), min_val(0.0); - complex z_val(0.0, 0.0); + std::complex z_val(0.0, 0.0); min_val = std::numeric_limits::max(); for(register int i=0; i < nm; i++){ @@ -915,7 +916,7 @@ inline const complex NRMat >::amin() const{ }else{ const size_t pozice = cublasIzamin((size_t)nn*mm, (cuDoubleComplex*)v, 1) - 1; TEST_CUBLAS("cublasIzamin"); - gpuget(1, sizeof(complex), v + pozice, &ret); + gpuget(1, sizeof(std::complex), v + pozice, &ret); } #endif return ret; @@ -1139,10 +1140,10 @@ void NRMat::resize(int n, int m) { * @return matrix \f$B\f$ where \f$\Re B=A\f$ and \f$\Im B = 0\f$ ******************************************************************************/ template -NRMat > complexify(const NRMat &rhs) { +NRMat > complexify(const NRMat &rhs) { NOT_GPU(rhs); - NRMat > r(rhs.nrows(), rhs.ncols(), rhs.getlocation()); + NRMat > r(rhs.nrows(), rhs.ncols(), rhs.getlocation()); for(register int i=0; i &a, NRMat *u, NRVec &s, NRMat *v, const bool vnotdagger, int m, int n) { diff --git a/quaternion.cc b/quaternion.cc index 134a154..3ed9294 100644 --- a/quaternion.cc +++ b/quaternion.cc @@ -331,10 +331,11 @@ case Euler_case('z','x','y'): } +// template void Quaternion::axis2normquat(const T *axis, const T &angle) { -T a = (T).5*angle; +T a = ((T).5)*angle; q[0]=cos(a); T s=sin(a); q[1]=axis[0]*s; diff --git a/quaternion.h b/quaternion.h index a9db23f..ab0333b 100644 --- a/quaternion.h +++ b/quaternion.h @@ -102,8 +102,8 @@ public: //C-style IO - void fprintf(FILE *f, const char *format) const {::fprintf(f,format,q[0],q[1],q[2],q[3]);}; - void sprintf(char *f, const char *format) const {::sprintf(f,format,q[0],q[1],q[2],q[3]);}; + int fprintf(FILE *f, const char *format) const {return ::fprintf(f,format,q[0],q[1],q[2],q[3]);}; + int sprintf(char *f, const char *format) const {return ::sprintf(f,format,q[0],q[1],q[2],q[3]);}; int fscanf(FILE *f, const char *format) const {return ::fscanf(f,format,q[0],q[1],q[2],q[3]);}; int sscanf(char *f, const char *format) const {return ::sscanf(f,format,q[0],q[1],q[2],q[3]);}; }; diff --git a/smat.cc b/smat.cc index db4dfdc..178613c 100644 --- a/smat.cc +++ b/smat.cc @@ -208,12 +208,12 @@ const NRSMat NRSMat::operator-() const { * @return modified copy of this matrix ******************************************************************************/ template <> -const NRSMat > NRSMat >::operator-() const { - NRSMat > result(nn, getlocation()); +const NRSMat > NRSMat >::operator-() const { + NRSMat > result(nn, getlocation()); #ifdef CUDALA if(location == cpu) { #endif - memcpy(result.v, v, NN2*sizeof(complex)); + memcpy(result.v, v, NN2*sizeof(std::complex)); cblas_zscal(NN2, &CMONE, result.v, 1); #ifdef CUDALA @@ -258,7 +258,7 @@ void NRSMat::randomize(const double &x) { * distribution. The real and imaginary parts are generated independently. ******************************************************************************/ template<> -void NRSMat >::randomize(const double &x) { +void NRSMat >::randomize(const double &x) { for(register size_t i=0; i NRSMat::operator*(const NRMat &rhs) const { * @return matrix produt \f$S\times{}A\f$ ******************************************************************************/ template<> -const NRMat > -NRSMat >::operator*(const NRMat > &rhs) const { +const NRMat > +NRSMat >::operator*(const NRMat > &rhs) const { #ifdef DEBUG - if (nn != rhs.nrows()) laerror("incompatible dimensions in NRSMat >::operator*(const NRMat > &)"); + if (nn != rhs.nrows()) laerror("incompatible dimensions in NRSMat >::operator*(const NRMat > &)"); #endif SAME_LOC(*this, rhs); - NRMat > result(nn, rhs.ncols(), getlocation()); + NRMat > result(nn, rhs.ncols(), getlocation()); #ifdef CUDALA if(location == cpu){ #endif @@ -429,14 +429,14 @@ const NRMat NRSMat::operator*(const NRSMat &rhs) const { * @return matrix produt \f$G\times{}H\f$ (not necessarily symmetric) ******************************************************************************/ template<> -const NRMat > -NRSMat >::operator*(const NRSMat > &rhs) const { +const NRMat > +NRSMat >::operator*(const NRSMat > &rhs) const { #ifdef DEBUG - if (nn != rhs.nn) laerror("incompatible dimensions in NRSMat >::operator*(const NRSMat > &)"); + if (nn != rhs.nn) laerror("incompatible dimensions in NRSMat >::operator*(const NRSMat > &)"); #endif SAME_LOC(*this, rhs); - NRMat > result(nn, nn, getlocation()); - NRMat > rhsmat(rhs); + NRMat > result(nn, nn, getlocation()); + NRMat > rhsmat(rhs); result = *this * rhsmat; return result; } @@ -477,11 +477,11 @@ const double NRSMat::dot(const NRSMat &rhs) const { * @return computed inner product ******************************************************************************/ template<> -const complex NRSMat >::dot(const NRSMat > &rhs) const { +const std::complex NRSMat >::dot(const NRSMat > &rhs) const { #ifdef DEBUG - if (nn != rhs.nn) laerror("incompatible dimensions in complex NRSMat >::dot(const NRSMat > &)"); + if (nn != rhs.nn) laerror("incompatible dimensions in std::complex NRSMat >::dot(const NRSMat > &)"); #endif - complex dot(0., 0.); + std::complex dot(0., 0.); SAME_LOC(*this, rhs); #ifdef CUDALA @@ -491,7 +491,7 @@ const complex NRSMat >::dot(const NRSMat #ifdef CUDALA }else{ const cuDoubleComplex _dot = cublasZdotc(NN2, (cuDoubleComplex*)v, 1, (cuDoubleComplex*)(rhs.v), 1); - dot = complex(cuCreal(_dot), cuCimag(_dot)); + dot = std::complex(cuCreal(_dot), cuCimag(_dot)); TEST_CUBLAS("cublasZdotc"); } #endif @@ -522,6 +522,7 @@ const double NRSMat::dot(const NRVec &rhs) const { TEST_CUBLAS("cublasDdot"); } #endif +return ret; } @@ -532,12 +533,12 @@ const double NRSMat::dot(const NRVec &rhs) const { * @return computed inner product ******************************************************************************/ template<> -const complex -NRSMat >::dot(const NRVec > &rhs) const { +const std::complex +NRSMat >::dot(const NRVec > &rhs) const { #ifdef DEBUG - if(NN2 != rhs.nn) laerror("incompatible dimensions in complex NRSMat >::dot(const NRVec > &)"); + if(NN2 != rhs.nn) laerror("incompatible dimensions in std::complex NRSMat >::dot(const NRVec > &)"); #endif - complex dot(0., 0.); + std::complex dot(0., 0.); SAME_LOC(*this, rhs); #ifdef CUDALA if(location == cpu){ @@ -547,7 +548,7 @@ NRSMat >::dot(const NRVec > &rhs) const { }else{ const cuDoubleComplex _dot = cublasZdotc(NN2, (cuDoubleComplex*)v, 1, (cuDoubleComplex*)rhs.v, 1); TEST_CUBLAS("cublasZdotc"); - dot = complex(cuCreal(_dot), cuCimag(_dot)); + dot = std::complex(cuCreal(_dot), cuCimag(_dot)); } #endif return dot; @@ -593,7 +594,7 @@ const double NRSMat::norm(const double scalar) const { * @param[in] scalar subtract this scalar value from the diagonal elements before the norm computation ******************************************************************************/ template<> -const double NRSMat< complex >::norm(const complex scalar) const { +const double NRSMat< std::complex >::norm(const std::complex scalar) const { if(!(scalar.real()) && !(scalar.imag())){ double ret(0.); #ifdef CUDALA @@ -611,7 +612,7 @@ const double NRSMat< complex >::norm(const complex scalar) const int k(0); double sum(0.); - complex tmp; + std::complex tmp; for(register int i=0; i::axpy(const double alpha, const NRSMat &x) { * \f[H \leftarrow \alpha G + H\f] ******************************************************************************/ template<> -void NRSMat >::axpy(const complex alpha, const NRSMat > & x) { +void NRSMat >::axpy(const std::complex alpha, const NRSMat > & x) { #ifdef DEBUG - if(nn != x.nn) laerror("incompatible dimensions in void NRSMat >::axpy(const complex , const NRSMat >&)"); + if(nn != x.nn) laerror("incompatible dimensions in void NRSMat >::axpy(const std::complex , const NRSMat >&)"); #endif SAME_LOC(*this, x); copyonwrite(); @@ -682,21 +683,21 @@ void NRSMat >::axpy(const complex alpha, const NRSMat -NRSMat >::NRSMat(const NRSMat &rhs, bool imagpart): nn(rhs.nrows()), count(new int(1)) { +NRSMat >::NRSMat(const NRSMat &rhs, bool imagpart): nn(rhs.nrows()), count(new int(1)) { //inconsistent in general case? const int nnp1 = nn*(nn + 1)/2; #ifdef CUDALA location = rhs.getlocation(); if(location == cpu){ #endif - v = new complex[nnp1]; - memset(v, 0, nnp1*sizeof(complex)); + v = new std::complex[nnp1]; + memset(v, 0, nnp1*sizeof(std::complex)); cblas_dcopy(nnp1, &rhs(0, 0), 1, ((double *)v) + (imagpart?1:0), 2); #ifdef CUDALA }else{ - v = (complex*) gpualloc(nnp1*sizeof(complex)); + v = (std::complex*) gpualloc(nnp1*sizeof(std::complex)); - complex *_val = gpuputcomplex(CZERO); + std::complex *_val = gpuputcomplex(CZERO); cublasZcopy(nnp1, (cuDoubleComplex*)_val, 0, (cuDoubleComplex*)v, 1); TEST_CUBLAS("cublasZcopy"); gpufree(_val); @@ -711,7 +712,7 @@ NRSMat >::NRSMat(const NRSMat &rhs, bool imagpart): nn(r * forced instantization in the corresponding object file ******************************************************************************/ template class NRSMat; -template class NRSMat >; +template class NRSMat >; template class NRSMat; template class NRSMat; diff --git a/smat.h b/smat.h index 38afed9..aba27a1 100644 --- a/smat.h +++ b/smat.h @@ -127,12 +127,12 @@ public: const T dot(const NRVec &rhs) const; const NRVec operator*(const NRVec &rhs) const {NRVec result(nn,rhs.getlocation()); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; - const NRVec > operator*(const NRVec > &rhs) const {NRVec > result(nn,rhs.getlocation()); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; + const NRVec > operator*(const NRVec > &rhs) const {NRVec > result(nn,rhs.getlocation()); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; const T* diagonalof(NRVec &, const bool divide = 0, bool cache = false) const; void gemv(const T beta, NRVec &r, const char trans, const T alpha, const NRVec &x) const {r.gemv(beta,*this,trans,alpha,x);}; - void gemv(const T beta, NRVec > &r, const char trans, const T alpha, const NRVec > &x) const {r.gemv(beta,*this,trans,alpha,x);}; + void gemv(const T beta, NRVec > &r, const char trans, const T alpha, const NRVec > &x) const {r.gemv(beta,*this,trans,alpha,x);}; inline const T& operator[](const size_t ij) const; inline T& operator[](const size_t ij); @@ -309,8 +309,8 @@ inline NRSMat & NRSMat::operator*=(const double &a) { * @return reference to the modified matrix ******************************************************************************/ template<> -inline NRSMat > & -NRSMat >::operator*=(const complex &a) { +inline NRSMat > & +NRSMat >::operator*=(const std::complex &a) { copyonwrite(); #ifdef CUDALA if(location == cpu){ @@ -320,7 +320,7 @@ NRSMat >::operator*=(const complex &a) { }else{ const cuDoubleComplex _a = make_cuDoubleComplex(a.real(), a.imag()); cublasZscal(NN2, _a, (cuDoubleComplex*)v, 1); - TEST_CUBLAS("cublasZscal");//"NRSMat >& NRSMat >::operator*=(const complex &)" + TEST_CUBLAS("cublasZscal");//"NRSMat >& NRSMat >::operator*=(const std::complex &)" } #endif return *this; @@ -404,9 +404,9 @@ inline NRSMat& NRSMat::operator+=(const NRSMat & rhs) { * @return reference to the modified matrix ******************************************************************************/ template<> -inline NRSMat >& NRSMat >::operator+=(const NRSMat > & rhs) { +inline NRSMat >& NRSMat >::operator+=(const NRSMat > & rhs) { #ifdef DEBUG - if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat >& NRSMat >::operator+=(const NRSMat > &)"); + if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat >& NRSMat >::operator+=(const NRSMat > &)"); #endif SAME_LOC(*this, rhs); copyonwrite(); @@ -418,7 +418,7 @@ inline NRSMat >& NRSMat >::operator+=(const NRSM #ifdef CUDALA }else{ cublasZaxpy(NN2, CUONE, (cuDoubleComplex*)(rhs.v), 1, (cuDoubleComplex*)v, 1); - TEST_CUBLAS("cublasZaxpy");//"NRSMat >& NRSMat >::operator+=(const NRSMat > &)" + TEST_CUBLAS("cublasZaxpy");//"NRSMat >& NRSMat >::operator+=(const NRSMat > &)" } #endif return *this; @@ -474,9 +474,9 @@ inline NRSMat& NRSMat::operator-=(const NRSMat& rhs) { * @return reference to the modified matrix ******************************************************************************/ template<> -inline NRSMat >& NRSMat >::operator-=(const NRSMat >& rhs) { +inline NRSMat >& NRSMat >::operator-=(const NRSMat >& rhs) { #ifdef DEBUG - if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat >& NRSMat >::operator-=(const NRSMat > &)"); + if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat >& NRSMat >::operator-=(const NRSMat > &)"); #endif SAME_LOC(*this, rhs); copyonwrite(); @@ -488,7 +488,7 @@ inline NRSMat >& NRSMat >::operator-=(const NRSM #ifdef CUDALA }else{ cublasZaxpy(NN2, CUMONE, (cuDoubleComplex*)(rhs.v), 1, (cuDoubleComplex*)v, 1); - TEST_CUBLAS("cublasZaxpy");//"NRSMat >& NRSMat >::operator-=(const NRSMat > &)" + TEST_CUBLAS("cublasZaxpy");//"NRSMat >& NRSMat >::operator-=(const NRSMat > &)" } #endif return *this; @@ -779,8 +779,8 @@ inline const double NRSMat::amin() const { * @return \f$A_{l,m}\f$ which maximizes \f$\left|\Re{}A_{i,j}\right| + \left|\Im{}A_{i,j}\right|\f$ ******************************************************************************/ template<> -inline const complex NRSMat< complex >::amax() const{ - complex ret(0., 0.); +inline const std::complex NRSMat< std::complex >::amax() const{ + std::complex ret(0., 0.); #ifdef CUDALA if(location == cpu){ #endif @@ -788,8 +788,8 @@ inline const complex NRSMat< complex >::amax() const{ #ifdef CUDALA }else{ const int pozice = cublasIzamax(NN2, (cuDoubleComplex*)v, 1) - 1; - TEST_CUBLAS("cublasIzamax");//"complex NRSMat >::amax()" - gpuget(1, sizeof(complex), v + pozice, &ret); + TEST_CUBLAS("cublasIzamax");//"std::complex NRSMat >::amax()" + gpuget(1, sizeof(std::complex), v + pozice, &ret); } #endif return ret; @@ -801,15 +801,15 @@ inline const complex NRSMat< complex >::amax() const{ * @return \f$A_{l,m}\f$ which minimizes \f$\left|\Re{}A_{i,j}\right| + \left|\Im{}A_{i,j}\right|\f$ ******************************************************************************/ template<> -inline const complex NRSMat >::amin() const{ - complex ret(0., 0.); +inline const std::complex NRSMat >::amin() const{ + std::complex ret(0., 0.); #ifdef CUDALA if(location == cpu){ #endif // izamin seems not to be supported int index(0); double val(0.0), min_val(0.0); - complex z_val(0.0, 0.0); + std::complex z_val(0.0, 0.0); min_val = std::numeric_limits::max(); for(register size_t i = 0; i < NN2; i++){ @@ -821,8 +821,8 @@ inline const complex NRSMat >::amin() const{ #ifdef CUDALA }else{ const int pozice = cublasIzamin(nn, (cuDoubleComplex*)v, 1) - 1; - TEST_CUBLAS("cublasIzamin");//"complex NRSMat >::amin()" - gpuget(1, sizeof(complex), v + pozice, &ret); + TEST_CUBLAS("cublasIzamin");//"std::complex NRSMat >::amin()" + gpuget(1, sizeof(std::complex), v + pozice, &ret); } #endif return ret; @@ -1056,10 +1056,10 @@ void NRSMat::moveto(const GPUID dest) { * @return matrix \f$B\f$ where \f$\Re B=A\f$ and \f$\Im B = 0\f$ ******************************************************************************/ template -NRSMat > complexify(const NRSMat &rhs) { +NRSMat > complexify(const NRSMat &rhs) { NOT_GPU(rhs); - NRSMat > r(rhs.nrows()); + NRSMat > r(rhs.nrows()); for(register int i = 0; i > complexify(const NRSMat &rhs) { ******************************************************************************/ /* template<> -NRSMat > complexify(const NRSMat &rhs) { - NRSMat > r(rhs.nrows(), rhs.getlocation()); +NRSMat > complexify(const NRSMat &rhs) { + NRSMat > r(rhs.nrows(), rhs.getlocation()); #ifdef CUDALA if(rhs.getlocation() == cpu){ #endif @@ -1082,7 +1082,7 @@ NRSMat > complexify(const NRSMat &rhs) { #ifdef CUDALA }else{ cublasDcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2); - TEST_CUBLAS("cublasDcopy");//"NRSMat > complexify(const NRSMat &)" + TEST_CUBLAS("cublasDcopy");//"NRSMat > complexify(const NRSMat &)" } #endif return r; diff --git a/sparsemat.cc b/sparsemat.cc index 41758e3..d5c4227 100644 --- a/sparsemat.cc +++ b/sparsemat.cc @@ -764,7 +764,7 @@ return r; } template<> -void NRMat >::gemm(const complex &beta, const SparseMat > &a, const char transa, const NRMat > &b, const char transb, const complex &alpha) +void NRMat >::gemm(const std::complex &beta, const SparseMat > &a, const char transa, const NRMat > &b, const char transb, const std::complex &alpha) { laerror("not implemented yet"); } @@ -1291,13 +1291,13 @@ template void SparseMat::permuteindices(const NRVec &p);\ INSTANTIZE(double) -INSTANTIZE(complex) //some functions are not OK for hermitean matrices, needs a revision!!! +INSTANTIZE(std::complex) //some functions are not OK for hermitean matrices, needs a revision!!! */ ////////////////////////////////////////////////////////////////////////////// //// forced instantization in the corresponding object file template class SparseMat; -template class SparseMat >; +template class SparseMat >; #define INSTANTIZE(T) \ template NRMat::NRMat(const SparseMat &rhs); \ @@ -1305,7 +1305,7 @@ template NRSMat::NRSMat(const SparseMat &rhs); \ template NRVec::NRVec(const SparseMat &rhs); INSTANTIZE(double) -INSTANTIZE(complex) +INSTANTIZE(std::complex) }//namespace diff --git a/sparsesmat.cc b/sparsesmat.cc index 586e3ed..cac53ba 100644 --- a/sparsesmat.cc +++ b/sparsesmat.cc @@ -357,7 +357,7 @@ INSTANTIZE(complex) //// forced instantization of functions in the header in the corresponding object file template class SparseSMat; -template class SparseSMat >; +template class SparseSMat >; /*activate this if needed template void SparseSMat >::put(int fd, bool dimen, bool transp) const; diff --git a/t.cc b/t.cc index 80e5c0c..984412c 100644 --- a/t.cc +++ b/t.cc @@ -922,6 +922,24 @@ cout < a,b,c; +cin >>a>>b; +c=a+b; +cout< NRVec::operator-() const { * @return the modified vector by value ******************************************************************************/ template<> -const NRVec > NRVec >::operator-() const { - NRVec > result(*this); +const NRVec > NRVec >::operator-() const { + NRVec > result(*this); result.copyonwrite(); #ifdef CUDALA if(location == cpu){ @@ -262,11 +262,11 @@ void NRVec::randomize(const double &x){ * \li \c true current vector is smaller than vector \c rhs ******************************************************************************/ template<> -void NRVec >::randomize(const double &x) { +void NRVec >::randomize(const double &x) { NOT_GPU(*this); for(register int i=0; i(x*(2.*random()/(1. + RAND_MAX) - 1.), x*(2.*random()/(1. + RAND_MAX) - 1.)); + v[i] = std::complex(x*(2.*random()/(1. + RAND_MAX) - 1.), x*(2.*random()/(1. + RAND_MAX) - 1.)); } } @@ -281,7 +281,7 @@ void NRVec >::randomize(const double &x) { * \li \c true current vector is smaller than vector \c rhs ******************************************************************************/ template<> -NRVec >::NRVec(const NRVec &rhs, bool imagpart): nn(rhs.size()){ +NRVec >::NRVec(const NRVec &rhs, bool imagpart): nn(rhs.size()){ count = new int; *count = 1; @@ -289,12 +289,12 @@ NRVec >::NRVec(const NRVec &rhs, bool imagpart): nn(rhs. location = rhs.getlocation(); if(location == cpu){ #endif - v = (complex*)new complex[nn]; - memset(v, 0, nn*sizeof(complex)); + v = (std::complex*)new std::complex[nn]; + memset(v, 0, nn*sizeof(std::complex)); cblas_dcopy(nn, &rhs[0], 1, ((double *)v) + (imagpart?1:0), 2); #ifdef CUDALA }else{ - v = (complex*) gpualloc(nn*sizeof(complex)); + v = (std::complex*) gpualloc(nn*sizeof(std::complex)); cublasZscal(nn, CUZERO, (cuDoubleComplex*)v, 1); TEST_CUBLAS("cublasZscal"); @@ -338,7 +338,7 @@ void NRVec::axpy(const double alpha, const NRVec &x) { * @param[in] x complex vector \f$\vec{x}\f$ ******************************************************************************/ template<> -void NRVec >::axpy(const complex alpha, const NRVec > &x){ +void NRVec >::axpy(const std::complex alpha, const NRVec > &x){ #ifdef DEBUG if (nn != x.nn) laerror("incompatible vectors"); #endif @@ -382,7 +382,7 @@ void NRVec::axpy(const double alpha, const double *x, const int stride){ * @param[in] stride sets the stride ******************************************************************************/ template<> -void NRVec >::axpy(const complex alpha, const complex *x, const int stride){ +void NRVec >::axpy(const std::complex alpha, const std::complex *x, const int stride){ NOT_GPU(*this); copyonwrite(); @@ -416,7 +416,7 @@ copyonwrite(); * @return reference to the modified vector ******************************************************************************/ template<> -NRVec >& NRVec >::operator=(const complex &a){ +NRVec >& NRVec >::operator=(const std::complex &a){ copyonwrite(); #ifdef CUDALA @@ -425,7 +425,7 @@ copyonwrite(); cblas_zcopy(nn, &a, 0, v, 1); #ifdef CUDALA }else{ - smart_gpu_set(nn, (complex)0, v); + smart_gpu_set(nn, (std::complex)0, v); } #endif return *this; @@ -492,7 +492,7 @@ NRVec& NRVec::normalize(double *norm){ * @return reference to the modified vector ******************************************************************************/ template<> -NRVec > & NRVec >::normalize(double *norm){ +NRVec > & NRVec >::normalize(double *norm){ double tmp(0.0); #ifdef CUDALA if(location == cpu){ @@ -566,8 +566,8 @@ void NRVec::gemv(const double beta, const NRMat &A, * @see gemm ******************************************************************************/ template<> -void NRVec >::gemv(const double beta, const NRMat &A, - const char trans, const double alpha, const NRVec > &x) { +void NRVec >::gemv(const double beta, const NRMat &A, + const char trans, const double alpha, const NRVec > &x) { #ifdef DEBUG if ((tolower(trans) == 'n'?A.ncols():A.nrows()) != x.size()){ laerror("incompatible vectors"); } #endif @@ -604,9 +604,9 @@ void NRVec >::gemv(const double beta, const NRMat &A, * @see gemm ******************************************************************************/ template<> -void NRVec >::gemv(const complex beta, - const NRMat > &A, const char trans, - const complex alpha, const NRVec > &x) { +void NRVec >::gemv(const std::complex beta, + const NRMat > &A, const char trans, + const std::complex alpha, const NRVec > &x) { #ifdef DEBUG if ((tolower(trans) == 'n'?A.ncols():A.nrows()) != x.size()){ laerror("incompatible vectors"); } #endif @@ -673,8 +673,8 @@ void NRVec::gemv(const double beta, const NRSMat &A, * @see gemm, NRSMat ******************************************************************************/ template<> -void NRVec >::gemv(const double beta, const NRSMat &A, - const char trans, const double alpha, const NRVec > &x) { +void NRVec >::gemv(const double beta, const NRSMat &A, + const char trans, const double alpha, const NRVec > &x) { #ifdef DEBUG if(A.ncols() != x.size()){ laerror("incompatible dimensions"); } #endif @@ -708,9 +708,9 @@ void NRVec >::gemv(const double beta, const NRSMat &A, * @see gemm, NRSMat ******************************************************************************/ template<> -void NRVec >::gemv(const complex beta, - const NRSMat > &A, const char trans, - const complex alpha, const NRVec > &x){ +void NRVec >::gemv(const std::complex beta, + const NRSMat > &A, const char trans, + const std::complex alpha, const NRVec > &x){ #ifdef DEBUG if(A.ncols() != x.size()) laerror("incompatible dimensions"); #endif @@ -726,7 +726,7 @@ void NRVec >::gemv(const complex beta, const cuDoubleComplex _alpha = make_cuDoubleComplex(alpha.real(), alpha.imag()); const cuDoubleComplex _beta = make_cuDoubleComplex(beta.real(), beta.imag()); - cublasZhpmv('U', A.ncols(), _alpha, (cuDoubleComplex*)((const complex*)A), (cuDoubleComplex*)(x.v), 1, _beta, (cuDoubleComplex*)(this->v), 1); + cublasZhpmv('U', A.ncols(), _alpha, (cuDoubleComplex*)((const std::complex*)A), (cuDoubleComplex*)(x.v), 1, _beta, (cuDoubleComplex*)(this->v), 1); TEST_CUBLAS("cublasZhpmv"); } #endif @@ -771,11 +771,11 @@ const NRMat NRVec::otimes(const NRVec &b,const bool conj * @param[in] scale complex scaling factor \f$\alpha\f$ ******************************************************************************/ template<> -const NRMat > -NRVec >::otimes(const NRVec > &b, const bool conj, const complex &scale) const { +const NRMat > +NRVec >::otimes(const NRVec > &b, const bool conj, const std::complex &scale) const { SAME_LOC(*this, b); - NRMat > result(0., nn, b.nn, this->getlocation()); + NRMat > result(0., nn, b.nn, this->getlocation()); #ifdef CUDALA if(location == cpu){ @@ -816,8 +816,8 @@ int NRVec::sort(int direction, int from, int to, int *perm) { } template<> -NRVec > complexify(const NRVec &rhs) { - NRVec > r(rhs.size(), rhs.getlocation()); +NRVec > complexify(const NRVec &rhs) { + NRVec > r(rhs.size(), rhs.getlocation()); #ifdef CUDALA if(rhs.getlocation() == cpu){ @@ -844,7 +844,7 @@ template void NRVec::put(int fd, bool dim, bool transp) const; \ template void NRVec::get(int fd, bool dim, bool transp); \ INSTANTIZE(double) -INSTANTIZE(complex) +INSTANTIZE(std::complex) INSTANTIZE(char) INSTANTIZE(short) INSTANTIZE(int) @@ -869,7 +869,7 @@ template<> const NRMat NRVec::otimes(const NRVec &b,const bool conj, co // Roman // following gemv are not implemented template<> void NRVec::gemv(const double beta, const SparseMat &a, const char trans, const double alpha, const NRVec &x, bool s) { laerror("gemv on unsupported types"); } -template<> void NRVec< complex >::gemv(const complex beta, const SparseMat< complex > &a, const char trans, const complex alpha, const NRVec< complex > &x, bool s) { laerror("gemv on unsupported types"); } +template<> void NRVec< std::complex >::gemv(const std::complex beta, const SparseMat< std::complex > &a, const char trans, const std::complex alpha, const NRVec< std::complex > &x, bool s) { laerror("gemv on unsupported types"); } INSTANTIZE_DUMMY(char) @@ -882,21 +882,21 @@ INSTANTIZE_DUMMY(unsigned short) INSTANTIZE_DUMMY(unsigned int) INSTANTIZE_DUMMY(unsigned long) INSTANTIZE_DUMMY(unsigned long long) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex) -INSTANTIZE_DUMMY(complex >) -INSTANTIZE_DUMMY(complex >) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex) +INSTANTIZE_DUMMY(std::complex >) +INSTANTIZE_DUMMY(std::complex >) template class NRVec; -template class NRVec >; +template class NRVec >; template class NRVec; template class NRVec; template class NRVec; diff --git a/vec.h b/vec.h index 6aa5546..3bc34aa 100644 --- a/vec.h +++ b/vec.h @@ -33,7 +33,7 @@ template void lawritemat(FILE *file, const T *a, int r, int c, /***************************************************************************//** * static constants used in several cblas-routines ******************************************************************************/ -const static complex CONE = 1.0, CMONE = -1.0, CZERO = 0.0; +const static std::complex CONE = 1.0, CMONE = -1.0, CZERO = 0.0; #ifdef CUDALA const static cuDoubleComplex CUONE = {1.,0.}, CUMONE = {-1.,0.}, CUZERO = {0.,0.}; #endif @@ -73,7 +73,7 @@ protected: public: friend class NRSMat; friend class NRMat; - template friend NRVec > complexify(const NRVec&); + template friend NRVec > complexify(const NRVec&); typedef T ROWTYPE; @@ -272,6 +272,7 @@ public: //! get the pointer to the underlying data structure inline operator T*(); + //! get the constant pointer to the underlying data structure inline operator const T*() const; @@ -963,14 +964,14 @@ NRVec & NRVec::operator|=(const NRVec &rhs) { * @see NRVec::copyonwrite() ******************************************************************************/ template -NRVec > complexify(const NRVec &rhs) { +NRVec > complexify(const NRVec &rhs) { NOT_GPU(rhs); - NRVec > r(rhs.size(), rhs.getlocation()); + NRVec > r(rhs.size(), rhs.getlocation()); for(register int i=0; i NRVec > complexify(const NRVec &rhs); +template<> NRVec > complexify(const NRVec &rhs); /***************************************************************************//** * routine for moving vector data between CPU and GPU memory @@ -1039,7 +1040,7 @@ inline NRVec& NRVec::operator+=(const double &a) { * @return reference to the modified vector ******************************************************************************/ template<> -inline NRVec >& NRVec >::operator+=(const complex &a) { +inline NRVec >& NRVec >::operator+=(const std::complex &a) { copyonwrite(); #ifdef CUDALA if(location == cpu){ @@ -1047,7 +1048,7 @@ inline NRVec >& NRVec >::operator+=(const comple cblas_zaxpy(nn, &CONE, &a, 0, v, 1); #ifdef CUDALA }else{ - complex *d = gpuputcomplex(a); + std::complex *d = gpuputcomplex(a); cublasZaxpy(nn, CUONE, (cuDoubleComplex *)d, 0, (cuDoubleComplex *)v, 1); TEST_CUBLAS("cublasZaxpy"); gpufree(d); @@ -1087,7 +1088,7 @@ inline NRVec& NRVec::operator-=(const double &a) { * @return reference to the modified vector ******************************************************************************/ template<> -inline NRVec >& NRVec >::operator-=(const complex &a) { +inline NRVec >& NRVec >::operator-=(const std::complex &a) { copyonwrite(); #ifdef CUDALA if(location == cpu){ @@ -1095,7 +1096,7 @@ inline NRVec >& NRVec >::operator-=(const comple cblas_zaxpy(nn, &CMONE, &a, 0, v, 1); #ifdef CUDALA }else{ - complex *d = gpuputcomplex(a); + std::complex *d = gpuputcomplex(a); cublasZaxpy(nn, CUMONE, (cuDoubleComplex *)d, 0, (cuDoubleComplex *)v, 1); TEST_CUBLAS("cublasZaxpy"); gpufree(d); @@ -1137,7 +1138,7 @@ inline NRVec& NRVec::operator+=(const NRVec &rhs) { * @return reference to the modified vector ******************************************************************************/ template<> -inline NRVec >& NRVec >::operator+=(const NRVec > &rhs) { +inline NRVec >& NRVec >::operator+=(const NRVec > &rhs) { #ifdef DEBUG if (nn != rhs.nn) laerror("incompatible dimensions"); #endif @@ -1189,7 +1190,7 @@ inline NRVec & NRVec::operator-=(const NRVec &rhs) { * @return reference to the modified vector ******************************************************************************/ template<> -inline NRVec >& NRVec >::operator-=(const NRVec > &rhs) { +inline NRVec >& NRVec >::operator-=(const NRVec > &rhs) { #ifdef DEBUG if (nn != rhs.nn) laerror("incompatible dimensions"); #endif @@ -1237,7 +1238,7 @@ inline NRVec& NRVec::operator*=(const double &a) { * @return reference to the modified vector ******************************************************************************/ template<> -inline NRVec >& NRVec >::operator*=(const complex &a) { +inline NRVec >& NRVec >::operator*=(const std::complex &a) { copyonwrite(); #ifdef CUDALA if(location == cpu){ @@ -1285,11 +1286,11 @@ inline const double NRVec::operator*(const NRVec &rhs) const { * @return \f$\sum_{i=1}^N\overbar{\vec{x}_i}\cdot\vec{y}_i\f$ ******************************************************************************/ template<> -inline const complex NRVec >::operator*(const NRVec< complex > &rhs) const { +inline const std::complex NRVec >::operator*(const NRVec< std::complex > &rhs) const { #ifdef DEBUG if(nn != rhs.nn) laerror("incompatible dimensions"); #endif - complex dot; + std::complex dot; SAME_LOC(*this, rhs); #ifdef CUDALA if(location == cpu){ @@ -1299,7 +1300,7 @@ inline const complex NRVec >::operator*(const NRVec< com }else{ const cuDoubleComplex val = cublasZdotc(nn, (cuDoubleComplex*)v, 1, (cuDoubleComplex*)rhs.v, 1); TEST_CUBLAS("cublasZdotc"); - dot = complex(cuCreal(val), cuCimag(val)); + dot = std::complex(cuCreal(val), cuCimag(val)); } #endif return dot; @@ -1324,8 +1325,8 @@ inline const double NRVec::dot(const double *y, const int stride) const * @return \f$\sum_{i=1}^N\vec{x}_{i}\cdot \overbar{y_{\mathrm{stride}\cdot(i-1) + 1}}\f$ ******************************************************************************/ template<> -inline const complex NRVec >::dot(const complex *y, const int stride) const { - complex dot; +inline const std::complex NRVec >::dot(const std::complex *y, const int stride) const { + std::complex dot; NOT_GPU(*this); cblas_zdotc_sub(nn, y, stride, v, 1, &dot); return dot; @@ -1358,7 +1359,7 @@ inline const double NRVec::asum() const { * @return the value of this sum ******************************************************************************/ template<> -inline const double NRVec >::asum() const { +inline const double NRVec >::asum() const { double ret(0.0); #ifdef CUDALA if(location == cpu){ @@ -1398,7 +1399,7 @@ inline const double NRVec::norm() const { * @return \f$\sum_{i=1}^N\left|\vec{x}_i\right|^2\f$ ******************************************************************************/ template<> -inline const double NRVec< complex >::norm() const { +inline const double NRVec< std::complex >::norm() const { double ret(0.); #ifdef CUDALA if(location == cpu){ @@ -1469,8 +1470,8 @@ inline const double NRVec::amin() const { * @return \f$\vec{v}_{j}\f$ which maximizes \f$\left\{\left|\Re{}\vec{v}_{i}\right|+\left|\Im{}\vec{v}_{i}\right|\right}\f$ ******************************************************************************/ template<> -inline const complex NRVec >::amax() const { - complex ret(0., 0.); +inline const std::complex NRVec >::amax() const { + std::complex ret(0., 0.); #ifdef CUDALA if(location == cpu){ #endif @@ -1479,7 +1480,7 @@ inline const complex NRVec >::amax() const { }else{ const int pozice = cublasIzamax(nn, (cuDoubleComplex*)v, 1) - 1; TEST_CUBLAS("cublasIzamax"); - gpuget(1, sizeof(complex), v + pozice, &ret); + gpuget(1, sizeof(std::complex), v + pozice, &ret); } #endif return ret; @@ -1491,15 +1492,15 @@ inline const complex NRVec >::amax() const { * @return \f$\vec{v}_{j}\f$ which minimizes \f$\left\{\left|\Re{}\vec{v}_{i}\right|+\left|\Im{}\vec{v}_{i}\right|\right}\f$ ******************************************************************************/ template<> -inline const complex NRVec >::amin() const { - complex ret(0., 0.); +inline const std::complex NRVec >::amin() const { + std::complex ret(0., 0.); #ifdef CUDALA if(location == cpu){ #endif // izamin seems not to be supported int index(0); double val(0.0), min_val(std::numeric_limits::max()); - complex z_val(0.0, 0.0); + std::complex z_val(0.0, 0.0); for(register int i=0; i < nn; i++){ z_val = v[i]; @@ -1511,7 +1512,7 @@ inline const complex NRVec >::amin() const { }else{ const int pozice = cublasIzamin(nn, (cuDoubleComplex*)v, 1) - 1; TEST_CUBLAS("cublasIzamin"); - gpuget(1, sizeof(complex), v + pozice, &ret); + gpuget(1, sizeof(std::complex), v + pozice, &ret); } #endif return ret; diff --git a/vecmat3.h b/vecmat3.h index 5891989..81e016b 100644 --- a/vecmat3.h +++ b/vecmat3.h @@ -82,8 +82,8 @@ public: Vec3& fast_normalize(void); const Vec3 operator*(const Mat3 &rhs) const; //C-style IO - void fprintf(FILE *f, const char *format) const {::fprintf(f,format,q[0],q[1],q[2]);}; - void sprintf(char *f, const char *format) const {::sprintf(f,format,q[0],q[1],q[2]);}; + int fprintf(FILE *f, const char *format) const {return ::fprintf(f,format,q[0],q[1],q[2]);}; + int sprintf(char *f, const char *format) const {return ::sprintf(f,format,q[0],q[1],q[2]);}; int fscanf(FILE *f, const char *format) const {return ::fscanf(f,format,q[0],q[1],q[2]);}; int sscanf(char *f, const char *format) const {return ::sscanf(f,format,q[0],q[1],q[2]);}; }; @@ -141,7 +141,7 @@ public: const Mat3 transpose() const {Mat3 r(*this); r.transposeme(); return r;}; const Mat3 inverse() const; //C-style IO - void fprintf(FILE *f, const char *format) const {::fprintf(f,format,q[0][0],q[0][1],q[0][2]); ::fprintf(f,format,q[1][0],q[1][1],q[1][2]); ::fprintf(f,format,q[2][0],q[2][1],q[2][2]);}; + int fprintf(FILE *f, const char *format) const {int n= ::fprintf(f,format,q[0][0],q[0][1],q[0][2]); n+=::fprintf(f,format,q[1][0],q[1][1],q[1][2]); n+=::fprintf(f,format,q[2][0],q[2][1],q[2][2]); return n;}; int fscanf(FILE *f, const char *format) const {return ::fscanf(f,format,q[0][0],q[0][1],q[0][2]) + ::fscanf(f,format,q[1][0],q[1][1],q[1][2]) + ::fscanf(f,format,q[2][0],q[2][1],q[2][2]);}; };