jiri/LA_library
1
0 Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

float support for smat

Browse Source
This commit is contained in:
jiri committed 2026-11-05 15:42:29 +01:00
1 parent ee8c6279e1
commit 5266767a4c
2 files changed
+657 -25
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

smat.cc
+448 -3
View File
@@ -167,7 +167,7 @@ const T* NRSMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const
if(divide){ if(divide){
for(register int i=0; i<nn; i++){ for(register int i=0; i<nn; i++){
const T a = v[(size_t)i*(i+1)/2+i]; const T a = v[(size_t)i*(i+1)/2+i];
if(a != 0.) r[i] /= a; if(a != (T)0) r[i] /= a;
} }
}else{ }else{
for(register int i=0; i<nn; i++) r[i] = v[(size_t)i*(i+1)/2+i]; for(register int i=0; i<nn; i++) r[i] = v[(size_t)i*(i+1)/2+i];
@@ -215,6 +215,26 @@ const NRSMat<double> NRSMat<double>::operator-() const {
return result; return result;
} }
template <>
const NRSMat<float> NRSMat<float>::operator-() const {
NRSMat<float> result(nn, getlocation());
#ifdef CUDALA
if(location == cpu){
#endif
memcpy(result.v, v, NN2*sizeof(float));
cblas_sscal(NN2, -1., result.v, 1);
#ifdef CUDALA
}else{
cublasScopy(NN2, v, 1, result.v, 1);
TEST_CUBLAS("cublasScopy");
cublasSscal(NN2, -1., result.v, 1);
TEST_CUBLAS("cublasSscal");
}
#endif
return result;
}
/***************************************************************************//** /***************************************************************************//**
* implements unary minus operator for this hermitian matrix * implements unary minus operator for this hermitian matrix
* @return modified copy of this matrix * @return modified copy of this matrix
@@ -240,6 +260,27 @@ const NRSMat<std::complex<double> > NRSMat<std::complex<double> >::operator-() c
return result; return result;
} }
template <>
const NRSMat<std::complex<float> > NRSMat<std::complex<float> >::operator-() const {
NRSMat<std::complex<float> > result(nn, getlocation());
#ifdef CUDALA
if(location == cpu) {
#endif
memcpy(result.v, v, NN2*sizeof(std::complex<float>));
cblas_cscal(NN2, &CMONEF, result.v, 1);
#ifdef CUDALA
}else{
cublasCcopy(NN2, (cuComplex*)v, 1, (cuComplex*)result.v, 1);
TEST_CUBLAS("cublasCcopy");
cublasCscal(NN2, CUMONE, (cuComplex*)result.v, 1);
TEST_CUBLAS("cublasCscal");
}
#endif
return result;
}
/***************************************************************************//** /***************************************************************************//**
* @return the sum of the diagonal elements * @return the sum of the diagonal elements
******************************************************************************/ ******************************************************************************/
@@ -266,6 +307,16 @@ void NRSMat<double>::randomize(const double &x) {
} }
} }
template<>
void NRSMat<float>::randomize(const float &x) {
NOT_GPU(*this);
copyonwrite();
for(size_t i=0; i<NN2; ++i){
v[i] = x*RANDDOUBLESIGNED();
}
}
/***************************************************************************//** /***************************************************************************//**
* Fill this hermitian matrix with pseudorandom numbers generated from uniform * Fill this hermitian matrix with pseudorandom numbers generated from uniform
* distribution. The real and imaginary parts are generated independently. * distribution. The real and imaginary parts are generated independently.
@@ -282,6 +333,18 @@ void NRSMat<std::complex<double> >::randomize(const double &x) {
} }
} }
template<>
void NRSMat<std::complex<float> >::randomize(const float &x) {
copyonwrite();
for(register size_t i=0; i<NN2; ++i) v[i].real(x*RANDDOUBLESIGNED());
for(register size_t i=0; i<NN2; ++i) v[i].imag(x*RANDDOUBLESIGNED());
for(register int i=0; i<nn; ++i){
for(register int j=0; j<=i; ++j){
if(i == j) v[i*(size_t)(i+1)/2+j].imag(0.); //hermitean
}
}
}
/***************************************************************************//** /***************************************************************************//**
* routine for formatted output via lawritemat * routine for formatted output via lawritemat
* @param[in] file pointer to <tt>FILE</tt> structure representing the output file * @param[in] file pointer to <tt>FILE</tt> structure representing the output file
@@ -390,6 +453,30 @@ const NRMat<double> NRSMat<double>::operator*(const NRMat<double> &rhs) const {
return result; return result;
} }
template<>
const NRMat<float> NRSMat<float>::operator*(const NRMat<float> &rhs) const {
#ifdef DEBUG
if(nn != rhs.nrows()) laerror("incompatible dimensions in NRMat<float> NRSMat<float>::operator*(const NRMat<float> &)");
#endif
SAME_LOC(*this, rhs);
NRMat<float> result(nn, rhs.ncols(), getlocation());
#ifdef CUDALA
if(location == cpu){
#endif
for(register int k = 0; k<rhs.ncols(); k++){
cblas_sspmv(CblasRowMajor, CblasLower, nn, 1.0, v, rhs[0] + k, rhs.ncols(), 0.0, result[0] + k, rhs.ncols());
}
#ifdef CUDALA
}else{
for(register int k = 0; k<rhs.ncols(); k++){
cublasSspmv('U', nn, 1.0, v, rhs[0] + k, rhs.ncols(), 0.0, result[0] + k, rhs.ncols());
TEST_CUBLAS("cublasSspmv");
}
}
#endif
return result;
}
/***************************************************************************//** /***************************************************************************//**
* multiply this real double-precision symmetric matrix \f$S\f$ stored in packed form * multiply this real double-precision symmetric matrix \f$S\f$ stored in packed form
@@ -425,6 +512,34 @@ NRSMat<std::complex<double> >::operator*(const NRMat<std::complex<double> > &rhs
return result; return result;
} }
template<>
const NRMat<std::complex<float> >
NRSMat<std::complex<float> >::operator*(const NRMat<std::complex<float> > &rhs) const {
#ifdef DEBUG
if (nn != rhs.nrows()) laerror("incompatible dimensions in NRSMat<std::complex<float> >::operator*(const NRMat<std::complex<float> > &)");
#endif
SAME_LOC(*this, rhs);
NRMat<std::complex<float> > result(nn, rhs.ncols(), getlocation());
#ifdef CUDALA
if(location == cpu){
#endif
for(register int k=0; k<rhs.ncols(); k++){
cblas_chpmv(CblasRowMajor, CblasLower, nn, &CONEF, v, rhs[0]+k, rhs.ncols(), &CZEROF, result[0]+k, rhs.ncols());
}
#ifdef CUDALA
}else{
for(register int k = 0; k<rhs.ncols(); k++){
cublasChpmv('U', nn,
CUONE, (cuComplex*)v, (cuComplex*)(rhs[0] + k), rhs.ncols(),
CUZERO, (cuComplex*)(result[0] + k), rhs.ncols());
TEST_CUBLAS("cublasSspmv");
}
}
#endif
return result;
}
/***************************************************************************//** /***************************************************************************//**
* multiply this real double-precision symmetric matrix \f$S\f$ stored in packed form * multiply this real double-precision symmetric matrix \f$S\f$ stored in packed form
* with real double-precision symmetric matrix \f$T\f$ * with real double-precision symmetric matrix \f$T\f$
@@ -478,6 +593,54 @@ const NRMat<double> NRSMat<double>::operator*(const NRSMat<double> &rhs) const {
return result; return result;
} }
template<>
const NRMat<float> NRSMat<float>::operator*(const NRSMat<float> &rhs) const {
#ifdef DEBUG
if (nn != rhs.nn) laerror("incompatible dimensions in NRMat<float> NRSMat<float>::operator*(const NRSMat<float> &)");
#endif
NRMat<float> result(0.0, nn, nn);
float *p, *q;
p = v;
for (int i=0; i<nn;i++) {
q = rhs.v;
for (int k=0; k<=i; k++) {
cblas_saxpy(k+1, *p++, q, 1, result[i], 1);
q += k+1;
}
}
p = v;
for (int i=0; i<nn;i++) {
q = rhs.v+1;
for (int j=1; j<nn; j++) {
result[i][j] += cblas_sdot(i+1<j ? i+1 : j, p, 1, q, 1);
q += j+1;
}
p += i+1;
}
p = v;
q = rhs.v;
for (int i=0; i<nn; i++) {
cblas_sger(CblasRowMajor, i, i+1, 1., p, 1, q, 1, result, nn);
p += i+1;
q += i+1;
}
q = rhs.v+3;
for (int j=2; j<nn; j++) {
p = v+1;
for (int i=1; i<j; i++) {
cblas_saxpy(i, *++q, p, 1, result[0]+j, nn);
p += i+1;
}
q += 2;
}
return result;
}
/***************************************************************************//** /***************************************************************************//**
* multiply this complex double-precision symmetric matrix \f$G\f$ stored in packed form * multiply this complex double-precision symmetric matrix \f$G\f$ stored in packed form
@@ -497,6 +660,19 @@ NRSMat<std::complex<double> >::operator*(const NRSMat<std::complex<double> > &rh
return result; return result;
} }
template<>
const NRMat<std::complex<float> >
NRSMat<std::complex<float> >::operator*(const NRSMat<std::complex<float> > &rhs) const {
#ifdef DEBUG
if (nn != rhs.nn) laerror("incompatible dimensions in NRSMat<std::complex<float> >::operator*(const NRSMat<std::complex<float> > &)");
#endif
SAME_LOC(*this, rhs);
NRMat<std::complex<float> > result(nn, nn, getlocation());
NRMat<std::complex<float> > rhsmat(rhs);
result = *this * rhsmat;
return result;
}
/***************************************************************************//** /***************************************************************************//**
* compute inner product of this real symmetric matrix \f$A\f$ with given real symmetric matrix \f$B\f$ * compute inner product of this real symmetric matrix \f$A\f$ with given real symmetric matrix \f$B\f$
@@ -524,6 +700,25 @@ const double NRSMat<double>::dot(const NRSMat<double> &rhs) const {
return ret; return ret;
} }
template<>
const float NRSMat<float>::dot(const NRSMat<float> &rhs) const {
float ret(0.);
#ifdef DEBUG
if (nn != rhs.nn) laerror("incompatible dimensions in float NRSMat<float>::dot(const NRSMat<float> &)");
#endif
SAME_LOC(*this, rhs);
#ifdef CUDALA
if(location == cpu){
#endif
ret = cblas_sdot(NN2, v, 1, rhs.v, 1);
#ifdef CUDALA
}else{
ret = cublasSdot(NN2, v, 1, rhs.v, 1);
}
#endif
return ret;
}
/***************************************************************************//** /***************************************************************************//**
* compute inner product of this complex symmetric matrix \f$A\f$ with given complex symmetric matrix \f$B\f$ * compute inner product of this complex symmetric matrix \f$A\f$ with given complex symmetric matrix \f$B\f$
@@ -554,6 +749,28 @@ const std::complex<double> NRSMat<std::complex<double> >::dot(const NRSMat<std::
return dot; return dot;
} }
template<>
const std::complex<float> NRSMat<std::complex<float> >::dot(const NRSMat<std::complex<float> > &rhs) const {
#ifdef DEBUG
if (nn != rhs.nn) laerror("incompatible dimensions in std::complex<float> NRSMat<std::complex<float> >::dot(const NRSMat<std::complex<float> > &)");
#endif
std::complex<float> dot(0., 0.);
SAME_LOC(*this, rhs);
#ifdef CUDALA
if(location == cpu){
#endif
cblas_cdotc_sub(NN2, v, 1, rhs.v, 1, &dot);
#ifdef CUDALA
}else{
const cuComplex _dot = cublasCdotc(NN2, (cuComplex*)v, 1, (cuComplex*)(rhs.v), 1);
dot = std::complex<float>(cuCreal(_dot), cuCimag(_dot));
TEST_CUBLAS("cublasCdotc");
}
#endif
return dot;
}
/***************************************************************************//** /***************************************************************************//**
* compute inner product of this real double-precision symmetric matrix \f$S\f$ of order \f$n\f$ * compute inner product of this real double-precision symmetric matrix \f$S\f$ of order \f$n\f$
@@ -581,6 +798,26 @@ const double NRSMat<double>::dot(const NRVec<double> &rhs) const {
return ret; return ret;
} }
template<>
const float NRSMat<float>::dot(const NRVec<float> &rhs) const {
float ret(0.0);
#ifdef DEBUG
if(NN2 != rhs.nn) laerror("incompatible dimensions in float NRSMat<float>::dot(const NRVec<float> &)");
#endif
SAME_LOC(*this, rhs);
#ifdef CUDALA
if(location == cpu){
#endif
ret = cblas_sdot(NN2, v, 1, rhs.v, 1);
#ifdef CUDALA
}else{
ret = cublasSdot(NN2, v, 1, rhs.v, 1);
TEST_CUBLAS("cublasSdot");
}
#endif
return ret;
}
/***************************************************************************//** /***************************************************************************//**
* compute inner product of this complex double-precision hermitian matrix \f$H\f$ of order \f$n\f$ * compute inner product of this complex double-precision hermitian matrix \f$H\f$ of order \f$n\f$
@@ -610,6 +847,28 @@ NRSMat<std::complex<double> >::dot(const NRVec<std::complex<double> > &rhs) cons
return dot; return dot;
} }
template<>
const std::complex<float>
NRSMat<std::complex<float> >::dot(const NRVec<std::complex<float> > &rhs) const {
#ifdef DEBUG
if(NN2 != rhs.nn) laerror("incompatible dimensions in std::complex<float> NRSMat<std::complex<float> >::dot(const NRVec<std::complex<float> > &)");
#endif
std::complex<float> dot(0., 0.);
SAME_LOC(*this, rhs);
#ifdef CUDALA
if(location == cpu){
#endif
cblas_cdotc_sub(NN2, v, 1, rhs.v, 1, &dot);
#ifdef CUDALA
}else{
const cuComplex _dot = cublasCdotc(NN2, (cuComplex*)v, 1, (cuComplex*)rhs.v, 1);
TEST_CUBLAS("cublasCdotc");
dot = std::complex<float>(cuCreal(_dot), cuCimag(_dot));
}
#endif
return dot;
}
/***************************************************************************//** /***************************************************************************//**
* compute the Frobenius norm of this real double-precision symmetric matrix * compute the Frobenius norm of this real double-precision symmetric matrix
* @param[in] scalar subtract this scalar value from the diagonal elements before the norm computation * @param[in] scalar subtract this scalar value from the diagonal elements before the norm computation
@@ -645,6 +904,37 @@ const double NRSMat<double>::norm(const double scalar) const {
return std::sqrt(sum); return std::sqrt(sum);
} }
template<>
const float NRSMat<float>::norm(const float scalar) const {
if(!scalar){
float ret(0.);
#ifdef CUDALA
if(location == cpu){
#endif
ret = cblas_snrm2(NN2, v, 1);
#ifdef CUDALA
}else{
ret = cublasSnrm2(NN2, v, 1);
TEST_CUBLAS("cublasSnrm2");
}
#endif
return ret;
}
NOT_GPU(*this);
float sum(0.);
int k(0);
for(register int i=0; i<nn; ++i){
for(register int j=0; j<=i; ++j) {
register float tmp = v[k++];
if(i == j) tmp -= scalar;
sum += tmp*tmp;
}
}
return std::sqrt(sum);
}
/***************************************************************************//** /***************************************************************************//**
* compute the Frobenius norm of this complex double-precision hermitian matrix * compute the Frobenius norm of this complex double-precision hermitian matrix
* @param[in] scalar subtract this scalar value from the diagonal elements before the norm computation * @param[in] scalar subtract this scalar value from the diagonal elements before the norm computation
@@ -680,6 +970,37 @@ const double NRSMat< std::complex<double> >::norm(const std::complex<double> sca
return std::sqrt(sum); return std::sqrt(sum);
} }
template<>
const float NRSMat< std::complex<float> >::norm(const std::complex<float> scalar) const {
if(!(scalar.real()) && !(scalar.imag())){
float ret(0.);
#ifdef CUDALA
if(location == cpu){
#endif
ret = cblas_scnrm2(NN2, v, 1);
#ifdef CUDALA
}else{
ret = cublasScnrm2(NN2, (cuComplex*)v, 1);
TEST_CUBLAS("cublasScnrm2");
}
#endif
return ret;
}
int k(0);
float sum(0.);
std::complex<float> tmp;
for(register int i=0; i<nn; ++i){
for(register int j=0; j<=i; ++j){
tmp = v[k++];
if (i == j) tmp -= scalar;
sum += tmp.real()*tmp.real() + tmp.imag()*tmp.imag();
}
}
return std::sqrt(sum);
}
/***************************************************************************//** /***************************************************************************//**
* symmetric (hermitian) rank-2 update of this matrix stored in packed form * symmetric (hermitian) rank-2 update of this matrix stored in packed form
@@ -709,6 +1030,16 @@ void NRSMat<double>::rank2update(const double alpha, const NRVec<double> &x, con
cblas_dspr2(CblasRowMajor, CblasLower, nn, alpha, x, 1, y, 1, v); cblas_dspr2(CblasRowMajor, CblasLower, nn, alpha, x, 1, y, 1, v);
} }
template<>
void NRSMat<float>::rank2update(const float alpha, const NRVec<float> &x, const NRVec<float> &y) {
#ifdef DEBUG
if(nn != x.size() || nn != y.size()) laerror("incompatible dimensions in NRSMat<float>::rank2update");
#endif
NOT_GPU(*this); NOT_GPU(x); NOT_GPU(y);
copyonwrite();
cblas_sspr2(CblasRowMajor, CblasLower, nn, alpha, x, 1, y, 1, v);
}
template<> template<>
void NRSMat<std::complex<double> >::rank2update(const std::complex<double> alpha, const NRVec<std::complex<double> > &x, const NRVec<std::complex<double> > &y) { void NRSMat<std::complex<double> >::rank2update(const std::complex<double> alpha, const NRVec<std::complex<double> > &x, const NRVec<std::complex<double> > &y) {
#ifdef DEBUG #ifdef DEBUG
@@ -719,6 +1050,16 @@ void NRSMat<std::complex<double> >::rank2update(const std::complex<double> alpha
cblas_zhpr2(CblasRowMajor, CblasLower, nn, &alpha, x, 1, y, 1, v); cblas_zhpr2(CblasRowMajor, CblasLower, nn, &alpha, x, 1, y, 1, v);
} }
template<>
void NRSMat<std::complex<float> >::rank2update(const std::complex<float> alpha, const NRVec<std::complex<float> > &x, const NRVec<std::complex<float> > &y) {
#ifdef DEBUG
if(nn != x.size() || nn != y.size()) laerror("incompatible dimensions in NRSMat<std::complex<float> >::rank2update");
#endif
NOT_GPU(*this); NOT_GPU(x); NOT_GPU(y);
copyonwrite();
cblas_chpr2(CblasRowMajor, CblasLower, nn, &alpha, x, 1, y, 1, v);
}
/***************************************************************************//** /***************************************************************************//**
* for this real double-precision symmetric matrix \f$S\f$ stored in packed form, * for this real double-precision symmetric matrix \f$S\f$ stored in packed form,
@@ -744,6 +1085,25 @@ void NRSMat<double>::axpy(const double alpha, const NRSMat<double> &x) {
#endif #endif
} }
template<>
void NRSMat<float>::axpy(const float alpha, const NRSMat<float> &x) {
#ifdef DEBUG
if(nn != x.nn) laerror("incompatible dimensions in void NRSMat<float>::axpy(const float, const NRSMat<float>&)");
#endif
SAME_LOC(*this, x);
copyonwrite();
#ifdef CUDALA
if(location == cpu){
#endif
cblas_saxpy(NN2, alpha, x.v, 1, v, 1);
#ifdef CUDALA
}else{
cublasSaxpy(NN2, alpha, x.v, 1, v, 1);
TEST_CUBLAS("cublasSaxpy");
}
#endif
}
/***************************************************************************//** /***************************************************************************//**
* for this complex double-precision hermitian matrix \f$H\f$ stored in packed form, * for this complex double-precision hermitian matrix \f$H\f$ stored in packed form,
@@ -771,6 +1131,27 @@ void NRSMat<std::complex<double> >::axpy(const std::complex<double> alpha, const
} }
template<>
void NRSMat<std::complex<float> >::axpy(const std::complex<float> alpha, const NRSMat<std::complex<float> > & x) {
#ifdef DEBUG
if(nn != x.nn) laerror("incompatible dimensions in void NRSMat<std::complex<float> >::axpy(const std::complex<float> , const NRSMat<std::complex<float> >&)");
#endif
SAME_LOC(*this, x);
copyonwrite();
#ifdef CUDALA
if(location == cpu){
#endif
cblas_caxpy(NN2, &alpha, x.v, 1, v, 1);
#ifdef CUDALA
}else{
const cuComplex _alpha = make_cuComplex(alpha.real(), alpha.imag());
cublasCaxpy(NN2, _alpha, (cuComplex*)x.v, 1, (cuComplex*)v, 1);
TEST_CUBLAS("cublasCaxpy");
}
#endif
}
/***************************************************************************//** /***************************************************************************//**
* create hermitian matrix \f$H\f$ from given real double-precision symmetric * create hermitian matrix \f$H\f$ from given real double-precision symmetric
* matrix \f$S\f$ * matrix \f$S\f$
@@ -803,6 +1184,32 @@ NRSMat<std::complex<double> >::NRSMat(const NRSMat<double> &rhs, bool imagpart):
#endif #endif
} }
template<>
NRSMat<std::complex<float> >::NRSMat(const NRSMat<float> &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 std::complex<float>[nnp1];
memset(v, 0, nnp1*sizeof(std::complex<float>));
cblas_scopy(nnp1, &rhs(0, 0), 1, ((float *)v) + (imagpart?1:0), 2);
#ifdef CUDALA
}else{
v = (std::complex<float>*) gpualloc(nnp1*sizeof(std::complex<float>));
std::complex<float> *_val = gpuputcomplex(CZEROF);
cublasCcopy(nnp1, (cuComplex*)_val, 0, (cuComplex*)v, 1);
TEST_CUBLAS("cublasCcopy");
gpufree(_val);
cublasScopy(nnp1, (float*)(&rhs(0,0)), 1, ((float*)v) + (imagpart?1:0), 2);
TEST_CUBLAS("cublasScopy");
}
#endif
}
/***************************************************************************//** /***************************************************************************//**
@@ -978,13 +1385,51 @@ NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::conjugateme() {
/***************************************************************************//**
* complexify a given real matrix \f$A\f$ (double and float, by BLAS; the complex default constructor zeroes the imaginary parts)
* @param[in] rhs matrix \f$A\f$ intended for this operation
* @return matrix \f$B\f$ where \f$\Re B=A\f$ and \f$\Im B = 0\f$
******************************************************************************/
template<>
NRSMat<std::complex<double> > complexify<double>(const NRSMat<double> &rhs) {
NRSMat<std::complex<double> > r(rhs.nrows(), rhs.getlocation());
#ifdef CUDALA
if(rhs.getlocation() == cpu){
#endif
cblas_dcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2);
#ifdef CUDALA
}else{
cublasDcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2);
TEST_CUBLAS("cublasDcopy");//"NRSMat<std::complex<double> > complexify(const NRSMat<double> &)"
}
#endif
return r;
}
template<>
NRSMat<std::complex<float> > complexify<float>(const NRSMat<float> &rhs) {
NRSMat<std::complex<float> > r(rhs.nrows(), rhs.getlocation());
#ifdef CUDALA
if(rhs.getlocation() == cpu){
#endif
cblas_scopy(rhs.size(), &(rhs[0]), 1, (float*)(&(r[0])), 2);
#ifdef CUDALA
}else{
cublasScopy(rhs.size(), &(rhs[0]), 1, (float*)(&(r[0])), 2);
TEST_CUBLAS("cublasScopy");//"NRSMat<std::complex<float> > complexify(const NRSMat<float> &)"
}
#endif
return r;
}
/***************************************************************************//** /***************************************************************************//**
* forced instantization in the corresponding object file * forced instantization in the corresponding object file
******************************************************************************/ ******************************************************************************/
template class NRSMat<double>; template class NRSMat<double>;
template class NRSMat<std::complex<double> >; template class NRSMat<std::complex<double> >;
//template class NRSMat<float>; template class NRSMat<float>;
//template class NRSMat<std::complex<float> >; template class NRSMat<std::complex<float> >;
template class NRSMat<long long>; template class NRSMat<long long>;
template class NRSMat<long>; template class NRSMat<long>;
template class NRSMat<int>; template class NRSMat<int>;
smat.h
+209 -22
View File
@@ -387,6 +387,22 @@ inline NRSMat<double> & NRSMat<double>::operator*=(const double &a) {
return *this; return *this;
} }
template<>
inline NRSMat<float> & NRSMat<float>::operator*=(const float &a) {
copyonwrite();
#ifdef CUDALA
if(location == cpu){
#endif
cblas_sscal(NN2, a, v, 1);
#ifdef CUDALA
}else{
cublasSscal(NN2, a, v, 1);
TEST_CUBLAS("cublasSscal");//"NRSMat<float>& NRSMat<float>::operator*=(const float &)"
}
#endif
return *this;
}
/***************************************************************************//** /***************************************************************************//**
* multiply this complex symmetric matrix with complex scalar value * multiply this complex symmetric matrix with complex scalar value
* @param[in] a complex multiplicative factor * @param[in] a complex multiplicative factor
@@ -410,6 +426,24 @@ NRSMat<std::complex<double> >::operator*=(const std::complex<double> &a) {
return *this; return *this;
} }
template<>
inline NRSMat<std::complex<float> > &
NRSMat<std::complex<float> >::operator*=(const std::complex<float> &a) {
copyonwrite();
#ifdef CUDALA
if(location == cpu){
#endif
cblas_cscal(NN2, &a, v, 1);
#ifdef CUDALA
}else{
const cuComplex _a = make_cuComplex(a.real(), a.imag());
cublasCscal(NN2, _a, (cuComplex*)v, 1);
TEST_CUBLAS("cublasCscal");//"NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator*=(const std::complex<float> &)"
}
#endif
return *this;
}
/***************************************************************************//** /***************************************************************************//**
* multiply this symmetric matrix of general type <code>T</code> stored in packed form * multiply this symmetric matrix of general type <code>T</code> stored in packed form
@@ -508,6 +542,27 @@ inline NRSMat<double>& NRSMat<double>::operator+=(const NRSMat<double> & rhs) {
return *this; return *this;
} }
template<>
inline NRSMat<float>& NRSMat<float>::operator+=(const NRSMat<float> & rhs) {
#ifdef DEBUG
if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat<float>& NRSMat<float>::operator+=(const NRSMat<float> &)");
#endif
SAME_LOC(*this, rhs);
copyonwrite();
#ifdef CUDALA
if(location == cpu){
#endif
cblas_saxpy(NN2, 1.0, rhs.v, 1, v, 1);
#ifdef CUDALA
}else{
cublasSaxpy(NN2, 1.0, rhs.v, 1, v, 1);
TEST_CUBLAS("cublasSaxpy");//" NRSMat<float>& NRSMat<float>::operator+=(const NRSMat<float> &)"
}
#endif
return *this;
}
/***************************************************************************//** /***************************************************************************//**
* add up this complex symmetric matrix with given symmetric matrix * add up this complex symmetric matrix with given symmetric matrix
* @param[in] rhs complex symmetric matrix to be added * @param[in] rhs complex symmetric matrix to be added
@@ -534,6 +589,27 @@ inline NRSMat<std::complex<double> >& NRSMat<std::complex<double> >::operator+=(
return *this; return *this;
} }
template<>
inline NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator+=(const NRSMat<std::complex<float> > & rhs) {
#ifdef DEBUG
if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator+=(const NRSMat<std::complex<float> > &)");
#endif
SAME_LOC(*this, rhs);
copyonwrite();
#ifdef CUDALA
if(location == cpu){
#endif
cblas_caxpy(NN2, &CONEF, rhs.v, 1, v, 1);
#ifdef CUDALA
}else{
cublasCaxpy(NN2, CUONE, (cuComplex*)(rhs.v), 1, (cuComplex*)v, 1);
TEST_CUBLAS("cublasCaxpy");//"NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator+=(const NRSMat<std::complex<float> > &)"
}
#endif
return *this;
}
/***************************************************************************//** /***************************************************************************//**
* add up this symmetric matrix of general type <code>T</code> with given symmetric matrix * add up this symmetric matrix of general type <code>T</code> with given symmetric matrix
* @param[in] rhs complex matrix of general type <code>T</code> to be added * @param[in] rhs complex matrix of general type <code>T</code> to be added
@@ -578,6 +654,27 @@ inline NRSMat<double>& NRSMat<double>::operator-=(const NRSMat<double>& rhs) {
return *this; return *this;
} }
template<>
inline NRSMat<float>& NRSMat<float>::operator-=(const NRSMat<float>& rhs) {
#ifdef DEBUG
if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat<float>& NRSMat<float>::operator-=(const NRSMat<float> &)");
#endif
SAME_LOC(*this, rhs);
copyonwrite();
#ifdef CUDALA
if(location == cpu){
#endif
cblas_saxpy(NN2, -1.0, rhs.v, 1, v, 1);
#ifdef CUDALA
}else{
cublasSaxpy(NN2, -1.0, rhs.v, 1, v, 1);
TEST_CUBLAS("cublasSaxpy");//" NRSMat<float>& NRSMat<float>::operator-=(const NRSMat<float> &)"
}
#endif
return *this;
}
/***************************************************************************//** /***************************************************************************//**
* subtracts given complex symmetric matrix from this complex symmetric matrix * subtracts given complex symmetric matrix from this complex symmetric matrix
* @param[in] rhs complex symmetric matrix to be subtracted * @param[in] rhs complex symmetric matrix to be subtracted
@@ -604,6 +701,27 @@ inline NRSMat<std::complex<double> >& NRSMat<std::complex<double> >::operator-=(
return *this; return *this;
} }
template<>
inline NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator-=(const NRSMat<std::complex<float> >& rhs) {
#ifdef DEBUG
if(nn != rhs.nn) laerror("incompatible dimensions in NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator-=(const NRSMat<std::complex<float> > &)");
#endif
SAME_LOC(*this, rhs);
copyonwrite();
#ifdef CUDALA
if(location == cpu){
#endif
cblas_caxpy(NN2, &CMONEF, rhs.v, 1, v, 1);
#ifdef CUDALA
}else{
cublasCaxpy(NN2, CUMONE, (cuComplex*)(rhs.v), 1, (cuComplex*)v, 1);
TEST_CUBLAS("cublasCaxpy");//"NRSMat<std::complex<float> >& NRSMat<std::complex<float> >::operator-=(const NRSMat<std::complex<float> > &)"
}
#endif
return *this;
}
/***************************************************************************//** /***************************************************************************//**
* subtracts given symmetric matrix of general type <code>T</code> from this symmetric matrix of type <code>T</code> * subtracts given symmetric matrix of general type <code>T</code> from this symmetric matrix of type <code>T</code>
* @param[in] rhs symmetric matrix of general type <code>T</code> to be subtracted * @param[in] rhs symmetric matrix of general type <code>T</code> to be subtracted
@@ -871,6 +989,24 @@ inline const double NRSMat<double>::amax() const {
return ret; return ret;
} }
template<>
inline const float NRSMat<float>::amax() const {
float ret(0.0);
#ifdef CUDALA
if(location == cpu){
#endif
ret = v[cblas_isamax(NN2, v, 1)];
#ifdef CUDALA
}else{
const int pozice = cublasIsamax(NN2, v, 1) ;
TEST_CUBLAS("cublasIsamax");//"float NRSMat<float>::amax()"
gpuget(1, sizeof(float), v + pozice, &ret);
}
#endif
return ret;
}
/***************************************************************************//** /***************************************************************************//**
* for this real symmetric matrix \f$A\f$, determine the * for this real symmetric matrix \f$A\f$, determine the
* first element with smallest absolute value * first element with smallest absolute value
@@ -901,6 +1037,31 @@ inline const double NRSMat<double>::amin() const {
return ret; return ret;
} }
template<>
inline const float NRSMat<float>::amin() const {
float ret(0.0);
#ifdef CUDALA
if(location == cpu){
#endif
// idamin seems not to be supported
float val(0.0);
int index(-1);
ret = std::numeric_limits<float>::max();
for(register size_t i = 0; i < NN2; i++){
val = std::abs(v[i]);
if(val < ret){ index = i; ret = val; }
}
ret = v[index];
#ifdef CUDALA
}else{
const int pozice = cublasIsamin(nn, v, 1);
TEST_CUBLAS("cublasIsamin");//"float NRSMat<float>::amin()"
gpuget(1, sizeof(float), v + pozice, &ret);
}
#endif
return ret;
}
/***************************************************************************//** /***************************************************************************//**
* for this complex symmetric matrix \f$A\f$, determine the * for this complex symmetric matrix \f$A\f$, determine the
* first element with largest "absolute value" * first element with largest "absolute value"
@@ -923,6 +1084,23 @@ inline const std::complex<double> NRSMat< std::complex<double> >::amax() const{
return ret; return ret;
} }
template<>
inline const std::complex<float> NRSMat< std::complex<float> >::amax() const{
std::complex<float> ret(0., 0.);
#ifdef CUDALA
if(location == cpu){
#endif
ret = v[cblas_icamax(NN2, v, 1)];
#ifdef CUDALA
}else{
const int pozice = cublasIcamax(NN2, (cuComplex*)v, 1);
TEST_CUBLAS("cublasIcamax");//"std::complex<float> NRSMat<std::complex<float> >::amax()"
gpuget(1, sizeof(std::complex<float>), v + pozice, &ret);
}
#endif
return ret;
}
/***************************************************************************//** /***************************************************************************//**
* for this complex symmetric matrix \f$A\f$, determine the * for this complex symmetric matrix \f$A\f$, determine the
* first element with smallest "absolute value" * first element with smallest "absolute value"
@@ -956,6 +1134,34 @@ inline const std::complex<double> NRSMat<std::complex<double> >::amin() const{
return ret; return ret;
} }
template<>
inline const std::complex<float> NRSMat<std::complex<float> >::amin() const{
std::complex<float> ret(0., 0.);
#ifdef CUDALA
if(location == cpu){
#endif
// izamin seems not to be supported
int index(0);
float val(0.0), min_val(0.0);
std::complex<float> z_val(0.0, 0.0);
min_val = std::numeric_limits<float>::max();
for(register size_t i = 0; i < NN2; i++){
z_val = v[i];
val = std::abs(z_val.real()) + std::abs(z_val.imag());
if(val < min_val){ index = i; min_val = val; }
}
ret = v[index];
#ifdef CUDALA
}else{
const int pozice = cublasIcamin(nn, (cuComplex*)v, 1);
TEST_CUBLAS("cublasIcamin");//"std::complex<float> NRSMat<std::complex<float> >::amin()"
gpuget(1, sizeof(std::complex<float>), v + pozice, &ret);
}
#endif
return ret;
}
/***************************************************************************//** /***************************************************************************//**
* @return pointer of general type T to the underlying data structure * @return pointer of general type T to the underlying data structure
******************************************************************************/ ******************************************************************************/
@@ -1215,28 +1421,9 @@ NRSMat<std::complex<T> > complexify(const NRSMat<T> &rhs) {
} }
/***************************************************************************//** //BLAS versions for double and float are in smat.cc (a full specialization defined in the header would be multiply defined)
* complexify a given double-precision real matrix \f$A\f$ template<> NRSMat<std::complex<double> > complexify<double>(const NRSMat<double> &rhs);
* @param[in] rhs matrix \f$A\f$ intended for this operation template<> NRSMat<std::complex<float> > complexify<float>(const NRSMat<float> &rhs);
* @return matrix \f$B\f$ where \f$\Re B=A\f$ and \f$\Im B = 0\f$
******************************************************************************/
/*
template<>
NRSMat<std::complex<double> > complexify(const NRSMat<double> &rhs) {
NRSMat<std::complex<double> > r(rhs.nrows(), rhs.getlocation());
#ifdef CUDALA
if(rhs.getlocation() == cpu){
#endif
cblas_dcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2);
#ifdef CUDALA
}else{
cublasDcopy(rhs.size(), &(rhs[0]), 1, (double*)(&(r[0])), 2);
TEST_CUBLAS("cublasDcopy");//"NRSMat<std::complex<double> > complexify(const NRSMat<double> &)"
}
#endif
return r;
}
*/
/***************************************************************************//** /***************************************************************************//**
* output operator * output operator
* @param[in,out] s output stream * @param[in,out] s output stream