simpleLAPmatrix.hh

Go to the documentation of this file.

#ifndef SIMPLELAPMATRIX_HH
#define SIMPLELAPMATRIX_HH
 
#include <vector>
#include <dune/grid/config.h>
#include <linalg/linearsystem.hh>
 
namespace Kaskade
{
 
class UNCMINVector
{
public:
  UNCMINVector(const std::vector<double>& v) : data(v) {};
  UNCMINVector(int n, double v=0.0) : data(n,v) {};
  UNCMINVector() {};
  UNCMINVector(double* start, double* stop) { for(double* s=start; s<stop; ++s) data.push_back(*s); };
  UNCMINVector& operator=(const UNCMINVector &A) { data=A.data; return *this; }
  UNCMINVector& operator=(const double &s) { for(int i=0; i<data.size(); ++i) data[i]=s; return *this; }
  UNCMINVector& operator*=(const double &s) { for(int i=0; i<data.size(); ++i) data[i]*=s; return *this; }
  double& operator()(int i) { return data[i-1]; }
  int size() { return data.size(); }
  void newsize(int n) {data.resize(n); }
 
  std::vector<double> data;
};
 
template<class Num, int offset=0>
class SLAPMatrix
{
public:
template<class Mat>
  SLAPMatrix(Mat const& mat)
  {
    col = mat.M();
    row = mat.N();
    lda = mat.N();
    data.reserve(col*row);
    // Insert elements in COLUMN-MAJOR ordering
    for(int i=0; i<col; ++i)
      for(int j=0; j<row; ++j)
        data.push_back(mat[j][i]);
  }
 
  SLAPMatrix(SLAPMatrix<Num,offset> const& mat)
  {
    col=mat.cols();
    row=mat.rows();
    lda=mat.LDA();
    data=mat.data;
  }
 
  SLAPMatrix(int row_, int col_, int lda_ = -1):
    data(std::max(lda_,row_)*col_,0), col(col_), row(row_), lda(std::max(lda_,row_)) 
  {
    assert(row>=0 && col>=0);
  }
    
 
  SLAPMatrix(MatrixAsTriplet<Num> const& mat)
  {
    col = mat.ncols();
    row = mat.nrows();
    lda = mat.nrows();
    data.resize(col*row);
    // This should be COLUMN-MAJOR ordering
    for(int i=0; i<mat.data.size(); ++i)
      (*this)(mat.ridx[i],mat.cidx[i]) = mat.data[i];
  }
 
  void print() const
  {
    std::cout << "[" << std::endl;
    for(int i=0; i<rows(); ++i)
    {
      for(int j=0; j<cols();++j)
        std::cout << (*this)(i,j) << ",";
      std::cout << std::endl;
    }
    std::cout << "]" << std::endl;
  }
 
 
  int cols() const {return col;}
  int rows() const {return row;}
  
  int size() {return col*row;}
  
  int LDA() const {return lda;}
 
template<class Mat>
  void toMatrix(Mat& M)
  {
    for(int i=0; i<rows(); ++i)
      for(int k=0; k<cols(); ++k)
        M[i][k] = (*this)(i,k);
  }
 
template<class RT>
  void toTriplet(MatrixAsTriplet<RT>& M)
  {
    for(int i=0; i<rows(); ++i)
      for(int k=0; k<cols(); ++k)
      {
        M.cidx.push_back(i);
        M.ridx.push_back(k);
        M.data.push_back((*this)(k,i));
      }
  }
 
  Num* ptrToData() {return &data[0];}
 
  Num const& operator()(int r, int c) const { return data[lda*(c-offset)+(r-offset)]; }
  Num&       operator()(int r, int c)       { return data[lda*(c-offset)+(r-offset)]; }
  
  int num_columns() const {return col;}
  int num_rows() const {return row;}
 
  SLAPMatrix<Num,offset> operator*=(Num alpha)
  {
    for(int i=0; i<rows(); ++i)
      for(int j=0; j<cols(); ++j)
        (*this)(i,j) *= alpha;
    return *this;
  }
 
  SLAPMatrix<Num,offset> scaleCols(std::vector<Num> scaling)
  {
    for(int i=0; i<rows(); ++i)
      for(int j=0; j<cols(); ++j)
        (*this)(i,j) *= scaling[j];
    return *this;
  }
 
  SLAPMatrix<Num,offset> scaleRows(std::vector<Num> scaling)
  {
    for(int i=0; i<rows(); ++i)
      for(int j=0; j<cols(); ++j)
        (*this)(i,j) *= scaling[i];
    return *this;
  }
 
private:
  std::vector<Num> data;
  int col,row,lda;
};
 
 
void LeastSquares(SLAPMatrix<double> a, std::vector<double> const& b, std::vector<double>& x);
 
void SymmetricEigenvalues(SLAPMatrix<double> a, std::vector<double>& eig);
 
template <class Scalar>
void pseudoinverse(SLAPMatrix<Scalar> A, SLAPMatrix<Scalar>& Ainv);
 
 
void MatrixMultiply(SLAPMatrix<double> A, std::vector<double>& in, std::vector<double>& out);
 
void TransposedMatrixMultiply(SLAPMatrix<double> A, std::vector<double>& in, std::vector<double>& out);
 
void MatrixKT_A_K(SLAPMatrix<double> A, SLAPMatrix<double>& K, SLAPMatrix<double>& out);
 
void MatrixMultiplication(SLAPMatrix<double> A, SLAPMatrix<double>& B, SLAPMatrix<double>& AB);
 
}  // namespace Kaskade
#endif