LA_library/sparsesmat.cc
2009-11-12 21:01:19 +00:00

209 lines
5.6 KiB
C++

/*
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"
namespace LA {
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 -- probably the computational bottleneck
for(i=0; i<prod.nrows(); ++i) for(j=0; j<prod.ncols(); ++j)
add(ai[i],bi[j],prod(i,j),false);
}
simplify();
}
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;
}
template <class T>
void SparseSMat<T>::axpy(const T alpha, const SparseSMat &x, const bool transp)
{
if(nn!=x.nn) laerror("incompatible matrix dimensions in SparseSMat::axpy");
if(alpha==(T)0) return;
copyonwrite();
for(SPMatindex i=0; i<nn; ++i) if(x.v[i])
{
if(!v[i]) v[i] = new std::map<SPMatindex,T>;
typename std::map<SPMatindex,T>::iterator p;
for(p=x.v[i]->begin(); p!=x.v[i]->end(); ++p) (*v[i])[p->first] = p->second * alpha;
}
simplify();
}
template <class T>
void SparseSMat<T>::gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const
{
if(nn!=r.size() || nn!= x.size()) laerror("incompatible matrix vector dimensions in SparseSMat::gemv");
r *= beta;
if(alpha == (T)0) return;
r.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) r[i] += x[p->first] * p->second * alpha ;
}
}
template <class T>
SparseSMat<T> & SparseSMat<T>::operator=(const T &a)
{
clear();
for(SPMatindex i=0; i<nn; ++i)
{
if(!v[i]) v[i] = new std::map<SPMatindex,T>;
(*v[i])[i] = a;
}
return *this;
}
template <class T>
SparseSMat<T> & SparseSMat<T>::operator+=(const T &a)
{
copyonwrite();
for(SPMatindex i=0; i<nn; ++i)
{
if(v[i])
{
typename std::map<SPMatindex,T>::iterator p;
p= v[i]->find(i);
if(p!=v[i]->end()) p->second+=a; else (*v[i])[i] = a;
}
else {v[i] = new std::map<SPMatindex,T>; (*v[i])[i] = a;}
}
return *this;
}
template <class T>
SparseSMat<T> & SparseSMat<T>::operator-=(const T &a)
{
copyonwrite();
for(SPMatindex i=0; i<nn; ++i)
{
if(v[i])
{
typename std::map<SPMatindex,T>::iterator p;
p= v[i]->find(i);
if(p!=v[i]->end()) p->second-=a; else (*v[i])[i] = -a;
}
else {v[i] = new std::map<SPMatindex,T>; (*v[i])[i] = -a;}
}
return *this;
}
template <class T>
typename LA_traits<T>::normtype SparseSMat<T>::norm(const T scalar) const
{
typename LA_traits<T>::normtype sum=0;
for(SPMatindex i=0; i<nn; ++i)
if(v[i])
{
typename std::map<SPMatindex,T>::iterator p;
p= v[i]->find(i);
if(p!=v[i]->end()) sum += LA_traits<T>::sqrabs(p->second - scalar);
else sum += LA_traits<T>::sqrabs(scalar);
}
else sum += LA_traits<T>::sqrabs(scalar); //missing diagonal element
return std::sqrt(sum);
}
#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); \
template void SparseSMat<T>::gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const; \
template void SparseSMat<T>::axpy(const T alpha, const SparseSMat &x, const bool transp); \
template SparseSMat<T> & SparseSMat<T>::operator=(const T &a); \
template SparseSMat<T> & SparseSMat<T>::operator+=(const T &a); \
template SparseSMat<T> & SparseSMat<T>::operator-=(const T &a); \
template LA_traits<T>::normtype SparseSMat<T>::norm(const T scalar) const; \
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> >;
}//namespace