dynamicMatrixOps.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) 2013-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef DYNAMICMATRIXOPS_HH
#define DYNAMICMATRIXOPS_HH
 
#include "linalg/dynamicMatrix.hh"
 
 
namespace Kaskade 
{
  template <class MEntry, class VEntry>
  auto operator*(DynamicMatrix<MEntry> const& A, Dune::DynamicVector<VEntry> const& x)
  {
    using YEntry = decltype(A[0][0]*x[0]);
    Dune::DynamicVector<YEntry> y(A.N());
    A.mv(x,y);
    return y;
  }
  
  template <class MEntry, class VEntry>
  auto operator*(DynamicMatrix<MEntry> const& A, Dune::BlockVector<VEntry> const& x)
  {
    using YEntry = decltype(A[0][0]*x[0]);
    Dune::BlockVector<YEntry> y(A.N());
    A.mv(x,y);
    return y;
  }
  
  template <class MEntry>
  auto operator*(typename Dune::FieldTraits<DynamicMatrix<MEntry>>::field_type x, DynamicMatrix<MEntry> A)
  {
    // Dune 2.5.1 appears to have a bug in densematrix.hh, where DenseMatrix<M>::field_type is defined 
    // as DenseMatVecTraits<M>::value_type (which is the entry type, and maybe a matrix). Thus, A *= x
    // fails, as x is interpreted to be the entry type instead of a scalar.
    for (int j=0; j<A.M(); ++j)
      for (int i=0; i<A.N(); ++i)
        A[i][j] *= x;
    return A;
  }
  
  template <class MEntry>
  auto operator-(DynamicMatrix<MEntry> const& A)
  {
    return -1.0*A;
  }
  
  template <class MEntry>
  auto operator-(DynamicMatrix<MEntry> A, DynamicMatrix<MEntry> const& B)
  {
    assert(A.N()==B.N() && A.M()==B.M());
//     A -= B; <--- strange, this runs into a loop (2016-09-29).
    for (int i=0; i<A.N(); ++i)
      for (int j=0; j<A.N(); ++j)
        A[i][j] -= B[i][j];
    return A;
  }
 
  template <class MEntry>
  auto operator+(DynamicMatrix<MEntry> A, DynamicMatrix<MEntry> const& B)
  {
    assert(A.N()==B.N() && A.M()==B.M());
    A += B;
    return A;
  }
 
  template <class EntryA, class EntryB>
  void matrixMatrixProduct(DynamicMatrix<EntryA> const& A, DynamicMatrix<EntryB> const& B,
                           DynamicMatrix<typename DynamicMatrixDetail::ProductTraits<EntryA,EntryB>::type>& C)
  {
    using namespace DynamicMatrixDetail;
    using EntryC = typename ProductTraits<EntryA,EntryB>::type;
    assert(A.M()==B.N() && C.N()==A.N() && C.M()==B.M());
 
    if (EntryTraits<EntryA>::lapackLayout && EntryTraits<EntryB>::lapackLayout
        && EntryTraits<EntryC>::lapackLayout && A.N()>0 && A.M()>0 && B.M()>0)
      gemm(false,false,A.N(),B.M(),A.M(),1.0,A.data(),A.lda(),B.data(),B.lda(),0.0,C.data(),C.lda());
    else
      for (int j=0; j<B.M(); ++j)
        for (int i=0; i<A.N(); ++i)
        {
          C[i][j] = 0;
          for (int l=0; l<A.M(); ++l)
            C[i][j] += A[i][l]*B[l][j];
        }
  }
  
  template <class EntryA, class EntryB>
  DynamicMatrix<typename DynamicMatrixDetail::ProductTraits<EntryA,EntryB>::type>
  operator*(DynamicMatrix<EntryA> const& A, DynamicMatrix<EntryB> const& B)
  {
    using namespace DynamicMatrixDetail;
    using EntryC = typename ProductTraits<EntryA,EntryB>::type;
 
    DynamicMatrix<EntryC> C(A.N(),B.M());
    matrixMatrixProduct(A,B,C);
 
    return C;
  }
 
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
  
  template <class MEntry>
  DynamicMatrix<MEntry> submatrix(DynamicMatrix<MEntry> const& A, 
                                  int rowstart, int rowend, int colstart, int colend)
  {
    assert(rowend>=rowstart && colend>=colstart);
    assert(0<=rowstart && 0<=colstart);
    assert(rowend<=A.N() && colend<=A.M());
    DynamicMatrix<MEntry> S(rowend-rowstart,colend-colstart);
    for (int j=0; j<colend-colstart; ++j)
      for (int i=0; i<rowend-rowstart; ++i)
        S[i][j] = A[i+rowstart][j+colstart];
    return S;
  }
  
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
  
  template <class EntryA, class EntryB>
  auto kron(DynamicMatrix<EntryA> const& A, DynamicMatrix<EntryB> const& B)
  {
    using EntryC = decltype(A[0][0]*B[0][0]);
    DynamicMatrix<EntryC> C(A.N()*B.N(),A.M()*B.M());
    
    for (int i=0; i<A.N(); ++i)
      for (int j=0; j<A.M(); ++j)
        for (int k=0; k<B.N(); ++k)
          for (int l=0; l<B.M(); ++l)
            C[i*B.N()+k][j*B.M()+l] = A[i][j]*B[k][l];
    return C;
  }
  
  // --------------------------------------------------------------------------------------------------
  // --------------------------------------------------------------------------------------------------
 
  namespace DynamicMatrixDetail
  {
    // TODO: move this into global horzcat/vertcat functions using constexpr if once we switch to C++17
    template <class EntryA, class EntryB>
    struct cat 
    {
      static auto horz(DynamicMatrix<EntryA> const& A, DynamicMatrix<EntryB> const& B)
      {
        using field_type = typename EntryA::field_type;
        static_assert(std::is_same<field_type,typename EntryB::field_type>::value,"A and B need to have the same underlying field type!");
        
        int const N = A.N()*EntryA::rows;
        assert(N==B.N()*EntryB::rows);
        
        using EntryAB = Dune::FieldMatrix<field_type,1,1>;
        
        DynamicMatrix<EntryAB> C(N,A.M()*EntryA::cols+B.M()*EntryB::cols);
        unflatten(C,horzcat(flatMatrix(A),flatMatrix(B)));  // TODO: a lot of copying around -> improve
        return C;
      }
      
      static auto vert(Dune::DynamicVector<EntryA> const& x, Dune::DynamicVector<EntryB> const& y)
      {
      }
    };
    
    template <class Entry>
    struct cat<Entry,Entry>
    {
      static auto horz(DynamicMatrix<Entry> const& A, DynamicMatrix<Entry> const& B) 
      {
        assert(A.N()==B.N());
        DynamicMatrix<Entry> C(A.N(),A.M()+B.M());
        for (int i=0; i<A.N(); ++i)
        {
          for (int j=0; j<A.M(); ++j)
            C[i][j] = A[i][j];
          for (int j=0; j<B.M(); ++j)
            C[i][j+A.M()] = B[i][j];
        }
        return C;
      }
      
      static auto vert(Dune::DynamicVector<Entry> const& x, Dune::DynamicVector<Entry> const& y)
      {
        Dune::DynamicVector<Entry> z(x.N()+y.N());
        for (int i=0; i<x.N(); ++i)
          z[i] = x[i];
        for (int i=0; i<y.N(); ++i)
          z[i+x.N()] = y[i];
        return z;
      }
    };
  }
 
  template <class EntryA, class EntryB>
  auto horzcat(DynamicMatrix<EntryA> const& A, DynamicMatrix<EntryB> const& B)
  {
    return DynamicMatrixDetail::cat<EntryA,EntryB>::horz(A,B);
  }
 
  template <class EntryX, class EntryY>
  auto vertcat(Dune::DynamicVector<EntryX> const& x, Dune::DynamicVector<EntryY> const& y)
  {
    return DynamicMatrixDetail::cat<EntryX,EntryY>::horz(x,y);
  }
 
  template <class Entry>
  DynamicMatrix<Entry> vertcat(DynamicMatrix<Entry> const& A, DynamicMatrix<Entry> const& B)
  {
    assert(A.M()==B.M());
    DynamicMatrix<Entry> C(A.N()+B.N(),A.M());
    for (int j=0; j<A.M(); ++j)
    {
      for (int i=0; i<A.N(); ++i)
        C[i][j] = A[i][j];
      for (int i=0; i<B.N(); ++i)
        C[i+A.N()][j] = B[i][j];
    }
    return C;
  }
  
  template <class Entry>
  auto dynamicUnitMatrix(Entry const& d, int n)
  {
    Entry zero = 0.0;
    DynamicMatrix<Entry> I(n,n,zero);  // make sure it is zero-initialized
    for (int i=0; i<n; ++i)
      I[i][i] = d;
    return I;
  }
  
  // --------------------------------------------------------------------------------------------------
  // --------------------------------------------------------------------------------------------------
  
  template <class Entry, int n>
  Dune::FieldMatrix<Entry,n,n> diag(Dune::FieldVector<Entry,n> const& d)
  {
    Dune::FieldMatrix<Entry,n,n> A;     // Create matrix
    A = Entry();                        // and zero-initialize it
    for (int i=0; i<n; ++i)
      A[i][i] = d[i];
    return A;
  }
}
 
#endif