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

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2008-04-16 12:56:02 +00:00
parent 9b60f9ddc0
commit b1f415fe04
2 changed files with 11 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
mat.cc
View File

@ -913,24 +913,6 @@ NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
// sum of rows
template<>
const NRVec<double> NRMat<double>::rsum() const
{
NRVec<double> result(mm);
for (int i=0; i<mm; i++) result[i] = cblas_dasum(nn,(*this)[0]+i,mm);
return result;
}
// sum of columns
template<>
const NRVec<double> NRMat<double>::csum() const
{
NRVec<double> result(nn);
for (int i=0; i<nn; i++) result[i] = cblas_dasum(mm, (*this)[i], 1);
return result;
}
// complex conjugate of Mat // complex conjugate of Mat
template<> template<>
NRMat<double> &NRMat<double>::conjugateme() {return *this;} NRMat<double> &NRMat<double>::conjugateme() {return *this;}
@ -1082,37 +1064,24 @@ void NRMat< complex<double> >::axpy(const complex<double> alpha,
// trace of Mat // trace of Mat
template<> template <typename T>
const double NRMat<double>::trace() const const T NRMat<T>::trace() const
{ {
#ifdef DEBUG #ifdef DEBUG
if (nn != mm) laerror("no-square matrix in Mat::trace()"); if (nn != mm) laerror("no-square matrix in Mat::trace()");
#endif #endif
return cblas_dasum(nn, (*this)[0], nn+1); T sum=0;
}
template<>
const complex<double> NRMat< complex<double> >::trace() const
{
#ifdef DEBUG
if (nn != mm) laerror("no-square matrix in Mat::trace()");
#endif
register complex<double> sum = CZERO;
for (int i=0; i<nn*nn; i+=(nn+1))
#ifdef MATPTR #ifdef MATPTR
sum += v[0][i]; for (int i=0; i<nn; ++i) sum += v[i][i];
#else #else
sum += v[i]; for (int i=0; i<nn*nn; i+=(nn+1)) sum += v[i];
#endif #endif
return sum; return sum;
} }
//get diagonal; for compatibility with large matrices do not return newly created object //get diagonal; for compatibility with large matrices do not return newly created object
//for non-square get diagonal of A^TA, will be used as preconditioner //for non-square get diagonal of A^TA, will be used as preconditioner
template<> template<>

14
vec.h
View File

@ -103,7 +103,8 @@ public:
const NRVec operator*(const NRSMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;}; const NRVec operator*(const NRSMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
const NRVec operator*(const SparseMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;}; const NRVec operator*(const SparseMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
const NRMat<T> operator|(const NRVec<T> &rhs) const; const NRMat<T> operator|(const NRVec<T> &rhs) const;
inline const T sum() const; //sum of its elements inline const T sum() const {T sum=0; for(int i=0; i<nn; i++) sum += v[i]; return sum;}; //sum of its elements
inline const T asum() const; //sum of its elements absolute values
inline const T dot(const T *a, const int stride=1) const; // ddot with a stride-vector inline const T dot(const T *a, const int stride=1) const; // ddot with a stride-vector
inline T & operator[](const int i); inline T & operator[](const int i);
inline const T & operator[](const int i) const; inline const T & operator[](const int i) const;
@ -447,18 +448,11 @@ inline const T NRVec<T>::operator*(const NRVec<T> &rhs) const
// Sum of elements // Sum of elements
template<> template<>
inline const double NRVec<double>::sum() const inline const double NRVec<double>::asum() const
{ {
return cblas_dasum(nn, v, 1); return cblas_dasum(nn, v, 1);
} }
template<>
inline const complex<double>
NRVec< complex<double> >::sum() const
{
complex<double> sum = CZERO;
for (int i=0; i<nn; i++) sum += v[i];
return sum;
}
// Dot product: x * y // Dot product: x * y
template<> template<>