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 2009-11-11 22:07:25 +00:00
parent eda9171eee
commit 91d130e0f7
2 changed files with 132 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
sparsesmat.cc Normal file
View File

@ -0,0 +1,101 @@
/*
LA: linear algebra C++ interface library
Copyright (C) 2008 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <string>
#include <cmath>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <errno.h>
#include "sparsesmat.h"
template <typename T>
void SparseSMat<T>::gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha)
{
(*this) *= beta;
if(alpha==(T)0) return;
if(a.nn!=b.nn || a.nn!=nn) laerror("incompatible sizes in SparseSMat::gemm");
copyonwrite();
for(SPMatindex k=0; k<nn; ++k) //summation loop
if(a.v[k] && b.v[k]) //nonempty in both
{
NRVec<T> av(a.v[k]->size());
NRVec<T> bv(b.v[k]->size());
NRVec<SPMatindex> ai(a.v[k]->size());
NRVec<SPMatindex> bi(b.v[k]->size());
//gather the data
typename std::map<SPMatindex,T>::iterator p;
int i,j;
for(p=a.v[k]->begin(), i=0; p!=a.v[k]->end(); ++p,++i)
{
ai[i] = p->first;
av[i] = p->second;
}
for(p=b.v[k]->begin(), i=0; p!=b.v[k]->end(); ++p,++i)
{
bi[i] = p->first;
bv[i] = p->second;
}
//make multiply via blas
NRMat<T> prod=av.otimes(bv,false,alpha);
//scatter the results
for(i=0; i<prod.nrows(); ++i) for(j=0; j<prod.ncols(); ++j)
add(ai[i],bi[j],prod(i,j),false);
}
simplify(); //erase elements below threshold
}
template <class T>
SparseSMat<T> & SparseSMat<T>::operator*=(const T &a)
{
if(!count) laerror("operator*= on undefined lhs");
if(a==(T)1) return *this;
if(a==(T)0) {clear(); return *this;}
copyonwrite();
for(SPMatindex i=0; i<nn; ++i) if(v[i])
{
typename std::map<SPMatindex,T>::iterator p;
for(p=v[i]->begin(); p!=v[i]->end(); ++p) p->second *= a;
}
return *this;
}
#define INSTANTIZE(T) \
template void SparseSMat<T>::gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha); \
template SparseSMat<T> & SparseSMat<T>::operator*=(const T &a); \
INSTANTIZE(double)
INSTANTIZE(complex<double>)
//// forced instantization of functions in the header in the corresponding object file
template class SparseSMat<double>;
template class SparseSMat<complex<double> >;

61
sparsesmat.h
View File

@ -62,24 +62,25 @@ public:
void simplify();
~SparseSMat();
inline int getcount() const {return count?*count:0;}
//
SparseSMat & operator*=(const T &a); //multiply by a scalar
inline const SparseSMat operator*(const T &rhs) const {return SparseSMat(*this) *= rhs;}
/*@@@to be done
inline const SparseSMat operator+(const T &rhs) const {return SparseSMat(*this) += rhs;}
inline const SparseSMat operator-(const T &rhs) const {return SparseSMat(*this) -= rhs;}
inline const SparseSMat operator*(const T &rhs) const {return SparseSMat(*this) *= rhs;}
inline const SparseSMat operator+(const SparseSMat &rhs) const {return SparseSMat(*this) += rhs;}
inline const SparseSMat operator-(const SparseSMat &rhs) const {return SparseSMat(*this) -= rhs;}
inline const SparseSMat operator*(const SparseSMat &rhs) const; //!!!NOT A GENERAL ROUTINE, JUST FOR THE CASES WHEN THE RESULT STAYS SYMMETRIC
SparseSMat & operator=(const T &a); //assign a to diagonal
SparseSMat & operator+=(const T &a); //assign a to diagonal
SparseSMat & operator-=(const T &a); //assign a to diagonal
SparseSMat & operator*=(const T &a); //multiply by a scalar
SparseSMat & operator+=(const SparseSMat &rhs);
SparseSMat & operator-=(const SparseSMat &rhs);
void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const;
void axpy(const T alpha, const SparseSMat &x, const bool transp=0); // this+= a*x
const typename LA_traits<T>::normtype norm(const T scalar=(T)0) const;
void add(const SPMatindex n, const SPMatindex m, const T elem, const bool both=true);
*/
inline const SparseSMat operator*(const SparseSMat &rhs) const {SparseSMat<T> r(nn); r.gemm(0,*this,'n',rhs,'n',1); return r;}; //!!!NOT A GENERAL ROUTINE, JUST FOR THE CASES WHEN THE RESULT STAYS SYMMETRIC
void gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha); //this := alpha*op( A )*op( B ) + beta*this !!!NOT A GENERAL ROUTINE, JUST FOR THE CASES WHEN THE RESULT STAYS SYMMETRIC
inline void add(const SPMatindex n, const SPMatindex m, const T elem, const bool both=true);
unsigned int length() const;
void transposeme() const {};
int nrows() const {return nn;}
@ -293,7 +294,6 @@ void SparseSMat<T>::add(const SPMatindex n, const SPMatindex m, const T elem, co
if(n>=nn || m>=nn) laerror("illegal index in SparseSMat::add()");
#endif
if(!v[n]) v[n] = new std::map<SPMatindex,T>;
if(!v[m]) v[m] = new std::map<SPMatindex,T>;
typename std::map<SPMatindex,T>::iterator p;
@ -301,6 +301,7 @@ p= v[n]->find(m);
if(p!=v[n]->end()) p->second+=elem; else (*v[n])[m] = elem;
if(n!=m && both) //add also transposed
{
if(!v[m]) v[m] = new std::map<SPMatindex,T>;
p= v[m]->find(n);
if(p!=v[m]->end()) p->second+=elem; else (*v[m])[n] = elem;
}
@ -317,42 +318,42 @@ for(SPMatindex i=0; i<nn; ++i) if(v[i])
//build a list since we are not sure whether erase from within the traversal loop is safe
std::list<SPMatindex> l;
typename std::map<SPMatindex,T>::iterator p;
for(p=v[i]->begin(); p!=v[i]->end(); ++p)
if(std::abs(p->second) < SPARSEEPSILON) l.push_front(p->first);
typename std::list<SPMatindex>::iterator q;
for(q=l.begin(); q!=l.end(); ++q) v[i]->erase(*q);
if(v[i]->size() == 0) delete v[i];
}
for(p=v[i]->begin(); p!=v[i]->end(); ++p)
if(std::abs(p->second) < SPARSEEPSILON) l.push_front(p->first);
typename std::list<SPMatindex>::iterator q;
for(q=l.begin(); q!=l.end(); ++q) v[i]->erase(*q);
if(v[i]->size() == 0) delete v[i];
}
}
template <typename T>
std::ostream & operator<<(std::ostream &s, const SparseSMat<T> &x)
{
SPMatindex n;
SPMatindex n;
n = x.nrows();
s << n << " "<< n<< std::endl;
n = x.nrows();
s << n << " "<< n<< std::endl;
typename SparseSMat<T>::iterator p(x);
for(; p.notend(); ++p) s << (int)p->row << ' ' << (int)p->col << ' ' << (typename LA_traits_io<T>::IOtype) p->elem << '\n';
s << "-1 -1\n";
return s;
s << "-1 -1\n";
return s;
}
template <class T>
std::istream& operator>>(std::istream &s, SparseSMat<T> &x)
{
SPMatindex n,m;
long i,j;
s >> n >> m;
if(n!=m) laerror("SparseSMat must be square");
x.resize(n);
s >> i >> j;
typename LA_traits_io<T>::IOtype tmp;
while(i>=0 && j>=0)
{
s>>tmp;
x.add(i,j,tmp,false);
{
SPMatindex n,m;
long i,j;
s >> n >> m;
if(n!=m) laerror("SparseSMat must be square");
x.resize(n);
s >> i >> j;
typename LA_traits_io<T>::IOtype tmp;
while(i>=0 && j>=0)
{
s>>tmp;
x.add(i,j,tmp,false);
s >> i >> j;
}
return s;