triplet.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef TRIPLET_HH
#define TRIPLET_HH
 
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <vector>
#include <memory>
#include <utility>
 
#include "dune/common/fmatrix.hh"
#include "dune/istl/bcrsmatrix.hh"
 
#include "linalg/crsutil.hh"
#include "linalg/umfpack_tools.hh"
 
namespace Kaskade
{
  // forward declaration
  template <class Entry, class Index>
  class NumaBCRSMatrix;
 
  template<class Scalar,class SparseIndexInt=int>
  class MatrixAsTriplet
  {
  public:
    typedef Scalar value_type;
    typedef typename std::vector<Scalar>::iterator iterator;
    typedef typename std::vector<Scalar>::const_iterator const_iterator;
    typedef typename std::vector<SparseIndexInt>::iterator index_iterator;
    typedef typename std::vector<SparseIndexInt>::const_iterator const_index_iterator;
 
    MatrixAsTriplet() {}
 
    explicit MatrixAsTriplet(size_t s) : ridx(s), cidx(s), data(s){}
 
    template <int n, int m, class Allocator>
    explicit MatrixAsTriplet(Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,n,m>,Allocator> const& other)
    {
      fillFromBCRSInterface<n,m>(other);
    }
 
    template <int n, int m, class Index>
    explicit MatrixAsTriplet(NumaBCRSMatrix<Dune::FieldMatrix<Scalar,n,m>,Index> const& other)
    {
      fillFromBCRSInterface<n,m>(other);
    }
 
    MatrixAsTriplet(std::vector<SparseIndexInt>&& row, std::vector<SparseIndexInt>&& col,
                    std::vector<Scalar>&& value)
    {
      assert(row.size()==col.size() && row.size()==value.size());
      assert(*std::min_element(row.begin(),row.end()) >=0);
      assert(*std::min_element(col.begin(),col.end()) >=0);
 
      ridx = std::move(row);
      cidx = std::move(col);
      data = std::move(value);
    }
 
    MatrixAsTriplet(SparseIndexInt nrows, SparseIndexInt ncols, Scalar const& value)
    {
      if(value != 0.0)
      {
        reserve(nrows*ncols);
        resize(0);
        for(size_t i=0; i<nrows; ++i)
          for(size_t k=0; k<ncols; ++k)
          {
            ridx.push_back(i);
            cidx.push_back(k);
            data.push_back(value);
          }
      }
    }
 
    MatrixAsTriplet(SparseIndexInt nrows, SparseIndexInt ncols, SparseIndexInt shift, Scalar const& value)
    {
      if(value != 0.0)
      {
        resize(0);
        for(size_t i=0; i<ncols; ++i)
        {
          size_t rowidx=i+shift;
          if(rowidx >= 0 && rowidx < nrows)
          {
            ridx.push_back(rowidx);
            cidx.push_back(i);
            data.push_back(value);
          }
        }
      }
    }
 
    MatrixAsTriplet(SparseIndexInt nrows, SparseIndexInt ncols, SparseIndexInt shift, std::vector<Scalar> const& value)
    {
      resize(0);
      for(int i=0; i<ncols; ++i)
      {
        int rowidx=i+shift;
        if(rowidx >= 0 && rowidx < nrows)
        {
          ridx.push_back(rowidx);
          cidx.push_back(i);
          data.push_back(value[i]);
        }
      }
    }
 
    MatrixAsTriplet(MatrixAsTriplet const& other) = default;
 
    MatrixAsTriplet(MatrixAsTriplet&& other) = default;
 
    template <class OtherScalar, class OtherIndex>
    MatrixAsTriplet(MatrixAsTriplet<OtherScalar,OtherIndex> const& A)
    : ridx(A.ridx.size())
    , cidx(A.cidx.size())
    , data(A.data.size())
    {
      std::copy(A.ridx.begin(),A.ridx.end(),ridx.begin());
      std::copy(A.cidx.begin(),A.cidx.end(),cidx.begin());
      std::copy(A.data.begin(),A.data.end(),data.begin());
    }
 
 
    void clear()
    {
      cidx.clear();
      ridx.clear();
      data.clear();
    }
 
 
    void setStartToZero()
    {
      SparseIndexInt r0 = *std::min_element(ridx.begin(),ridx.end());
      SparseIndexInt c0 = *std::min_element(cidx.begin(),cidx.end());
      if(r0!=0 || c0!=0)
        shiftIndices(-r0,-c0);
    }
 
    SparseIndexInt nrows() const
    {
      if(ridx.size()==0) return 0;
      return *std::max_element(ridx.begin(),ridx.end())+1;
    }
 
    SparseIndexInt N() const { return nrows(); }
 
    SparseIndexInt ncols() const
    {
      if(cidx.size()==0) return 0;
      return *std::max_element(cidx.begin(),cidx.end())+1;
    }
 
    SparseIndexInt M() const { return ncols(); }
 
    void resize(size_t s) 
    {
      ridx.resize(s); 
      cidx.resize(s); 
      data.resize(s); 
    }
 
    void reserve(size_t s) 
    {
      ridx.reserve(s); 
      cidx.reserve(s); 
      data.reserve(s); 
    }
 
    size_t nnz() const { return ridx.size(); }
 
 
    template <class X, class Y>
    void axpy(Y& out, X const& in, Scalar alpha=1.0,int nvectors=1) const
    {
      assert(ridx.size()==cidx.size() && ridx.size()==data.size());
      assert(in.size() >= ncols()*nvectors);
      assert(out.size() >= nrows()*nvectors);
      SparseIndexInt cols=ncols();
      SparseIndexInt rows=nrows();
      for(int j=0; j<nvectors; ++j)
        for(size_t i=0; i<ridx.size(); ++i)
        {
          assert(ridx[i]+j*rows < out.size());
          assert(cidx[i]+j*cols < in.size());
          out[ridx[i]+j*rows] += alpha*in[cidx[i]+j*cols]*data[i];
        }
    }
 
 
    template <class X, class Y>
    void atxpy(Y& y, X const& x, Scalar alpha=1.0) const
    {
      assert(ridx.size()==cidx.size() && ridx.size()==data.size());
      assert(y.size() == ncols());
      assert(x.size() == nrows());
      for(size_t i=0; i<data.size(); ++i) y[cidx[i]] += alpha*data[i]*x[ridx[i]];
    }
 
 
    template <class X, class Y>
    void usmtv(Scalar const alpha, X const& x, Y& y) const
    {
      atxpy(y,x,alpha);
    }
 
 
    template <class X, class Y>
    void usmv(Scalar const alpha, X const& x, Y& y) const
    {
      axpy(y, x, alpha);
    }
 
 
    template <class X, class Y>
    void ax(Y& out, X const& in) const
    {
      assert(ridx.size()==cidx.size() && ridx.size()==data.size());
      assert(in.size() >= ncols());
      out.resize(0);
      out.resize(nrows(),0.0);
      for(size_t i=0; i<ridx.size(); ++i)
      {
        assert(ridx[i] < out.size());
        assert(cidx[i] < in.size());
        out[ridx[i]] += in[cidx[i]]*data[i];
      }
    }
 
    template <class X, class Y>
    void mv(X const& x, Y& y) const
    {
      axpy(y,x);
    }
 
    void shiftIndices(SparseIndexInt row, SparseIndexInt col)
    {
      assert(ridx.size()==cidx.size() && ridx.size()==data.size());
      for(size_t i=0; i< ridx.size(); ++i)
      {
        ridx[i] += row;
        cidx[i] += col;
        assert(ridx[i] >= 0 && cidx[i] >= 0);
      }
    }
 
    void addEntry(SparseIndexInt row, SparseIndexInt col, Scalar const& value)
    {
      ridx.push_back(row);
      cidx.push_back(col);
      data.push_back(value);
    }
 
    MatrixAsTriplet& operator+=(MatrixAsTriplet const& m)
    {
      if(ridx.size()==0)
      {
        ridx=m.ridx;
        cidx=m.cidx;
        data=m.data;
        return *this;
      }
      ridx.insert(ridx.end(),m.ridx.begin(),m.ridx.end());
      cidx.insert(cidx.end(),m.cidx.begin(),m.cidx.end());
      data.insert(data.end(),m.data.begin(),m.data.end());
      
      // conversion to compresses row and back to triplet to remove multiple entries
      
      std::vector<SparseIndexInt> Ap, Ai;
      std::vector<Scalar> Az;
      umfpack_triplet_to_col(nrows(), ncols(), ridx, cidx, data, Ap, Ai, Az);
      {
        std::vector<SparseIndexInt> nullInt1, nullInt2;
        std::vector<Scalar> nullDouble;
        cidx.swap(nullInt1);
        ridx.swap(nullInt2);
        data.swap(nullDouble);
      }
      umfpack_col_to_triplet(Ap, cidx);
      data.insert(data.begin(),&Az[0],&Az[Ap.back()]);
      ridx.insert(ridx.begin(),&Ai[0],&Ai[Ap.back()]);
      
      return *this;
    }
 
    MatrixAsTriplet& operator-=(MatrixAsTriplet const& m)
                      {
      if(ridx.size()==0)
      {
        ridx=m.ridx;
        cidx=m.cidx;
        data=m.data;
        return *this;
      }
      int s=data.size();
      ridx.insert(ridx.end(),m.ridx.begin(),m.ridx.end());
      cidx.insert(cidx.end(),m.cidx.begin(),m.cidx.end());
      data.insert(data.end(),m.data.begin(),m.data.end());
 
      for(int i=s; i < data.size(); ++i)
        data[i] *= -1.0;
      // conversion to compresses row and back to triplet to remove multiple entries
 
      std::vector<SparseIndexInt> Ap, Ai;
      std::vector<Scalar> Az;
      umfpack_triplet_to_col(nrows(), ncols(), ridx, cidx, data, Ap, Ai, Az);
      {
        std::vector<SparseIndexInt> nullInt1, nullInt2;
        std::vector<Scalar> nullDouble;
        cidx.swap(nullInt1);
        ridx.swap(nullInt2);
        data.swap(nullDouble);
      }
      umfpack_col_to_triplet(Ap, cidx);
      data.insert(data.begin(),&Az[0],&Az[Ap.back()]);
      ridx.insert(ridx.begin(),&Ai[0],&Ai[Ap.back()]);
 
      return *this;
                      }
 
 
    MatrixAsTriplet& operator*=(Scalar const& scalar)
    {
      for(auto& v: data)
        v *= scalar;
 
      return *this;
    }
 
    MatrixAsTriplet& operator/=(Scalar const& scalar)
    {
      return (*this) *= 1/scalar;
    }
 
    MatrixAsTriplet& operator=(MatrixAsTriplet const& m) = default;
 
    MatrixAsTriplet& operator=(MatrixAsTriplet&& other) = default;
 
    void scaleRows(std::vector<Scalar>const& scaling)
    {
      for(size_t i=0; i<ridx.size();++i)
        data[i] *= scaling[ridx[i]];
    }
 
    MatrixAsTriplet& transpose()
    {
      // simple: just exchange row and column indices.
      std::swap(ridx,cidx);
      return *this;
    }
 
    void deleteLowerTriangle()
    {
      size_t count(0);
      for(size_t i=0; i<ridx.size(); ++i)
      {
        if(ridx[i] <= cidx[i]) count++;
      }
      std::vector<SparseIndexInt> r2,c2;
      std::vector<Scalar> d2;
      r2.reserve(count);
      c2.reserve(count);
      d2.reserve(count);
      for(size_t i=0; i<ridx.size(); ++i)
      {
        if(ridx[i] <= cidx[i])
        {
          r2.push_back(ridx[i]);
          c2.push_back(cidx[i]);
          d2.push_back(data[i]);
        }
      }
      std::swap(r2,ridx);
      std::swap(c2,cidx);
      std::swap(d2,data);
    }
 
    std::ostream& print(std::ostream& out = std::cout, double threshold=0.0) const
    {
      out << "[" << std::endl;
      for(size_t i=0; i<ridx.size(); ++i) if(std::fabs(data[i]) > threshold)
        out << ridx[i] << "," << cidx[i] << "," << data[i] << ";" << std::endl;
      out << "]" << std::endl;
      return out;
    }
 
    void toColumns(std::vector<std::vector<Scalar> >& colsvector) const
    {
      SparseIndexInt rows=nrows();
      colsvector.resize(ncols());
      for(size_t i=0; i<colsvector.size();++i)
      {
        colsvector[i].resize(0);
        colsvector[i].resize(rows,0.0);
      }
      for(size_t i=0; i<data.size(); ++i)
        colsvector[cidx[i]][ridx[i]]=data[i];
    }
 
    void toRows(std::vector<std::vector<Scalar> >& rows) const
    {
      SparseIndexInt cols = ncols();
      rows.resize(nrows());
      std::for_each(rows.begin(),rows.end(),
          [&cols](std::vector<Scalar>& row)
          {
            row.resize(0);
            row.resize(cols,0.0);
          });
      for(size_t i=0; i<data.size(); ++i) rows[ridx[i]][cidx[i]] = data[i];
    }
 
    void toVector(std::vector<Scalar>& colsvector)
    {
      SparseIndexInt rows=nrows();
      colsvector.resize(0);
      colsvector.resize(nrows()*ncols(),0.0);
      for(size_t i=0; i<data.size(); ++i)
        colsvector[cidx[i]*rows+ridx[i]]=data[i];
    }
 
    void addColumn(std::vector<Scalar>& colsvector, size_t position)
    {
      for(size_t i=0; i<colsvector.size();++i)
      {
        ridx.push_back(i);
        cidx.push_back(position);
        data.push_back(colsvector[i]);
      }
    }
 
 
    template<class Mat>
    void addToMatrix(Mat& mat)
    {
      for(size_t i=0; i<data.size(); ++i)
        mat[ridx[i]][cidx[i]] += data[i];
    }
 
    void erase(SparseIndexInt i)
    {
      ridx.erase(std::begin(ridx) + i);
      cidx.erase(std::begin(cidx) + i);
      data.erase(std::begin(data) + i);
    }
 
    MatrixAsTriplet lumped()
    {
      SparseIndexInt noRows = nrows(), nnz = data.size();
      std::vector<Scalar> diagonalEntries(noRows,0);
 
      // sum up row entries
      for(SparseIndexInt row=0; row<noRows; ++row)
      {
        SparseIndexInt i=0;
        while(i<data.size())
        {
          if(ridx[i] == row) diagonalEntries[row] += data[i];
          else ++i;
        }
      }
 
      assert(nrows()==0 && ncols()==0); // matrix should be empty now
 
      MatrixAsTriplet<Scalar,SparseIndexInt> result;
      for(SparseIndexInt i=0; i<noRows; ++i) result.addEntry(i,i,diagonalEntries[i]);
      return result;
    }
 
    iterator begin(){ return data.begin(); }
 
    const_iterator begin() const { return data.begin(); }
 
    iterator end() { return data.end(); }
 
    const_iterator end() const { return data.end(); }
 
    size_t size() const { return data.size(); }
 
    void toFile(const char* filename, unsigned int precision=10) const
    {
      std::ofstream myfile;
      myfile.precision(precision);
      myfile.open(filename);
      for(size_t i=0; i<data.size(); ++i)
        myfile << ridx[i] << " " << cidx[i] << " " << data[i] << std::endl;
      myfile.close();
    }  
    
    template <int bcrsN=1>
    std::unique_ptr<Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,bcrsN,bcrsN>>> toBCRS() const
    {
      typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,bcrsN,bcrsN> > BCRSMat;
      
      // get column indices for each row
      SparseIndexInt noRows = nrows()/bcrsN,
                     noCols = ncols()/bcrsN;
      std::vector<std::vector<SparseIndexInt> > bcrsids(noRows);
      getBCRSIndices<bcrsN>(noRows, ridx, cidx, bcrsids);
 
      SparseIndexInt nonz = 0;
      for (auto const& s : bcrsids)
        nonz += s.size();
      std::unique_ptr<BCRSMat> result(new BCRSMat(noRows,noCols,nonz,BCRSMat::row_wise));
 
      // create sparsity pattern
      for (typename BCRSMat::CreateIterator row=result->createbegin(); row!=result->createend(); ++row)
        for(SparseIndexInt col = 0; col < bcrsids[row.index()].size(); ++col) row.insert(bcrsids[row.index()][col]);
 
      // fill matrix
      for (SparseIndexInt i=0; i<data.size(); ++i)
      {
        size_t blockRow = ridx[i]/bcrsN, blockCol = cidx[i]/bcrsN;
        (*result)[blockRow][blockCol][ridx[i]-bcrsN*blockRow][cidx[i]-bcrsN*blockCol] = data[i];
      }
 
      return result;
    }
    
    bool isSymmetric() const
    {
      for(size_t i=0; i<data.size(); ++i)
      {
        // check if symmetric entry exists
        auto tmp = findEntry(cidx[i],ridx[i]);
        if(tmp.second == false || std::fabs(data[i]-data[tmp.first]) > 1e-6)
          return false;
      }
 
      return true;
    }
 
    std::pair<size_t,bool> findEntry(SparseIndexInt row, SparseIndexInt col) const
    {
      for(size_t i=0; i<data.size(); ++i) 
        if(ridx[i]==row && cidx[i]==col) 
          return std::make_pair(i,true);
      
      return std::make_pair(0,false);
    }
 
    
    std::vector<SparseIndexInt>  ridx;
    std::vector<SparseIndexInt> cidx;
    std::vector<Scalar>  data;
 
  private:
    template <int bcrsN, class SparseInt>
    void getBCRSIndices(SparseInt nrows, std::vector<SparseInt> const& ridx, std::vector<SparseInt> const& cidx, std::vector<std::vector<SparseInt> >& rowColumnIndices) const
    {
      rowColumnIndices.clear();
      rowColumnIndices.resize(nrows);
      for(SparseInt i=0; i<ridx.size(); ++i) rowColumnIndices[ridx[i]/bcrsN].push_back(cidx[i]/bcrsN);
 
      for(std::vector<SparseInt>& s : rowColumnIndices)
      {
        std::sort(s.begin(),s.end());
        auto last = std::unique(s.begin(),s.end());
        s.erase(last,s.end());
      }
    }
 
    template <int n, int m, class Matrix>
    void fillFromBCRSInterface(Matrix const& other)
    {
      reserve(other.nonzeroes()*n*m);
 
      auto rend = other.end();
      for (auto r=other.begin(); r!=rend; ++r)  // iterate over rows
      {
        SparseIndexInt ri = r.index() * n;      // block row start index in matrix with scalar entries
        auto cend = r->end();
        for (auto c=r->begin(); c!=cend; ++c)   // iterate over columns
        {
          SparseIndexInt ci = c.index() * m;    // block col start index in matrix with scalar entries
          auto const& entry = *c;
           
          for(size_t k=0; k<n; ++k)             // iterate over block entries
            for(size_t l=0; l<m; ++l)
            {
              ridx.push_back(ri+k);
              cidx.push_back(ci+l);
              data.push_back(entry[k][l]);
            }
        }
      }
    }
    
    template <class Iterator>
    struct Read
    {
      Read(Iterator& in_): in(in_) {}
 
      template <class Element>
      void operator()(Element& e) const { in = vectorFromSequence(e,in); }
 
    private:
      Iterator& in;
    };
 
    template <class Iterator>
    struct Write
    {
      Write(Iterator& out_): out(out_) {}
 
      template <class Element>
      void operator()(Element const& e) const { out = vectorToSequence(e,out); }
 
    private:
      Iterator& out;
    };
  };
 
  template<class Scalar>
  void deleteLowerTriangle(std::vector<int>& ridx, std::vector<int>& cidx, std::vector<Scalar>& data)
  {
    size_t count(0);
    for(size_t i=0; i<ridx.size(); ++i)
      if(ridx[i] <= cidx[i])
        count++;
 
    std::vector<int> r2,c2;
    std::vector<Scalar> d2;
    r2.reserve(count);
    c2.reserve(count);
    d2.reserve(count);
 
    for(size_t i=0; i<ridx.size(); ++i)
      if(ridx[i] <= cidx[i]) {
        r2.push_back(ridx[i]);
        c2.push_back(cidx[i]);
        d2.push_back(data[i]);
      }
 
    std::swap(r2,ridx);
    std::swap(c2,cidx);
    std::swap(d2,data);
  }
 
  template <typename Scalar, typename SparseIndexInt>
  std::ostream& operator<<(std::ostream &s, MatrixAsTriplet<Scalar,SparseIndexInt> const& mat)
  {
    return mat.print(s,0);
  }
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
 
} // end of namespace Kaskade
 
#endif