matlab.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-2018 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef MATLAB_HH
#define MATLAB_HH
 
#include <fstream>
#include <string>
 
#include "linalg/triplet.hh"
 
namespace Kaskade
{
  
  struct vec {
      size_t N;
      std::vector<double> data;
 
      vec(size_t N) : N(N), data(N) {};
 
      vec(std::ifstream &in) {
          in.read(reinterpret_cast<char*>(&N), sizeof(size_t));
          data.resize(N);
          in.read(reinterpret_cast<char*>(data.data()), N * sizeof(double));
      }
 
      void write(std::ofstream &out) {
          out.write(reinterpret_cast<char*>(&N), sizeof(size_t));
          out.write(reinterpret_cast<char*>(data.data()), N * sizeof(double));
      }
  };
 
  struct csr_matrix {
      size_t N, nnz;
      std::vector<size_t> row_ptrs;
      std::vector<size_t> col_idxs;
      std::vector<double> values;
 
      csr_matrix(size_t N, size_t nnz) : N(N), nnz(nnz), row_ptrs(N + 1, 0), col_idxs(nnz), values(nnz)
      {}
 
      csr_matrix(std::ifstream &in) {
          in.read(reinterpret_cast<char*>(&N), sizeof(size_t));
          in.read(reinterpret_cast<char*>(&nnz), sizeof(size_t));
          row_ptrs.resize(N + 1);
          col_idxs.resize(nnz);
          values.resize(nnz);
          in.read(reinterpret_cast<char*>(row_ptrs.data()), (N + 1) * sizeof(size_t));
          in.read(reinterpret_cast<char*>(col_idxs.data()), nnz * sizeof(size_t));
          in.read(reinterpret_cast<char*>(values.data()), nnz * sizeof(double));
      }
 
      void write(std::ofstream &out) {
          out.write(reinterpret_cast<char*>(&N), sizeof(size_t));
          out.write(reinterpret_cast<char*>(&nnz), sizeof(size_t));
          out.write(reinterpret_cast<char*>(row_ptrs.data()), (N + 1) * sizeof(size_t));
          out.write(reinterpret_cast<char*>(col_idxs.data()), nnz * sizeof(size_t));
          out.write(reinterpret_cast<char*>(values.data()), nnz * sizeof(double));
      }
  };
 
  template <class Entry, class Index, class VEntry>
  void writeToMatlabBinary(NumaBCRSMatrix<Entry,Index> const& A, Dune::BlockVector<VEntry> const& b, std::string const& basename, std::string const& path, int precision=16)
  {
    vec b_binary(b.N());
    for (size_t i = 0; i < b.N(); i++) {
        b_binary.data[i] = b[i];
    }
 
    int n = Entry::rows;
    int m = Entry::cols;
 
    csr_matrix A_binary(b.N(), 0);
    int count = 0;
    for (auto row = A.begin(); row!=A.end(); ++row)
    {
      for (int i=0; i<n; ++i)
      {
        auto ridx = row.index()*n + i;
        for (auto col=row.begin(); col!=row.end(); ++col)
        {
          for (int j=0; j<m; ++j)
          {
            auto val = (*col)[i][j];
            if (val != 0.0) {
              auto cidx = col.index()*m + j;
              A_binary.col_idxs.emplace_back(cidx);
              A_binary.values.emplace_back(val);
              count++;
            }
          }
        }
        A_binary.row_ptrs[ridx+1] = count;
      }
    }
    A_binary.nnz = count;
 
    // write the vector and matrix to disk
    std::string fname_vec = path+"/"+basename+"_vec.bin";
    std::string fname_A = path+"/"+basename+"_A.bin";
    std::ofstream out(fname_vec.c_str(), std::ios::binary);
    b_binary.write(out);
    out.close();
 
    out.open(fname_A.c_str(), std::ios::binary);
    A_binary.write(out);
    out.close();
  }
 
  template <class Entry, class Index>
  void writeToMatlab(NumaBCRSMatrix<Entry,Index> const& A, std::string const& basename, std::string const& outPath,int precision=16)
  {
    std::string fname = outPath + basename + ".m";
    std::ofstream f(fname.c_str());
    f.precision(precision);
    
    // Write header.
    f << "function A = " << basename << '\n';
    
    // write matrix in triplet format
    f << "data = [\n";
    int n = Entry::rows;
    int m = Entry::cols;
    
    for (auto row = A.begin(); row!=A.end(); ++row)
    {
      for (int i=0; i<n; ++i)
      {
        auto ridx = row.index()*n + i;
        for (auto col=row.begin(); col!=row.end(); ++col)
        {
          for (int j=0; j<m; ++j)
          {
            auto cidx = col.index()*m + j;
            f << ridx+1 << ' ' << cidx+1 << ' ' << (*col)[i][j] << std::endl;
          }
        }
      }
    }
    f << "];\nA = sparse(data(:,1),data(:,2),data(:,3)," << n*A.N() << "," << m*A.M() << ");\n";
  }
  
  template <class Entry, class Index, class VEntry>
  void writeToMatlab(NumaBCRSMatrix<Entry,Index> const& A, Dune::BlockVector<VEntry> const& b, std::string const& basename, int precision=16)
  {
    std::string fname = basename + ".m";
    std::ofstream f(fname.c_str());
    f.precision(precision);
    
    // Write vector.
    f << "function [A,b] = " << basename << '\n'
    << " b = [\n";
    for (size_t i=0; i<b.N(); ++i)
      f << b[i] << std::endl;
    f << "];\n";
    
    // write matrix in triplet format
    f << "data = [\n";
    int n = Entry::rows;
    int m = Entry::cols;
    
    for (auto row = A.begin(); row!=A.end(); ++row)
    {
      for (int i=0; i<n; ++i)
      {
        auto ridx = row.index()*n + i;
        for (auto col=row.begin(); col!=row.end(); ++col)
        {
          for (int j=0; j<m; ++j)
          {
            auto cidx = col.index()*m + j;
            f << ridx+1 << ' ' << cidx+1 << ' ' << (*col)[i][j] << std::endl;
          }
        }
      }
    }
    f << "];\nA = sparse(data(:,1),data(:,2),data(:,3)," << n*A.N() << "," << m*A.M() << ");\n";
  }
  
  template <class Assembler>
  void writeToMatlab(Assembler const& assembler, std::string const& basename, int precision=16)
  {
    std::string fname = basename + ".m";
    std::ofstream f(fname.c_str());
    f.precision(precision);
 
    f << "function [A,b] = " << basename << '\n'
        << " b = [\n";
 
    assembler.toSequence(0,Assembler::TestVariableSetDescription::noOfVariables,
        std::ostream_iterator<typename Assembler::Scalar>(f,"\n"));
    f << "\n];\ndata = [\n";
 
    typedef MatrixAsTriplet<typename Assembler::Scalar> Matrix;
 
    Matrix A = assembler.template get<Matrix>(false);
    
    for (size_t i=0; i<A.ridx.size(); ++i)
      f << A.ridx[i]+1 << ' ' << A.cidx[i]+1 << ' ' << A.data[i] << '\n';
    f << "];\nA = sparse(data(:,1),data(:,2),data(:,3)," << A.N() << "," << A.M() << ");\n";
  }
 
  template <class Assembler>
  void writeToMatlab(Assembler const& assembler, std::string const& path, std::string const& basename, int precision=16)
  {
    std::string fname = path + ".m";
    std::ofstream f(fname.c_str());
    f.precision(precision);
 
    f << "function [A,b] = " << basename << '\n'
        << " b = [\n";
 
    assembler.toSequence(0,Assembler::TestVariableSetDescription::noOfVariables,
        std::ostream_iterator<typename Assembler::Scalar>(f,"\n"));
    f << "\n];\ndata = [\n";
 
    typedef MatrixAsTriplet<typename Assembler::Scalar> Matrix;
 
    Matrix A = assembler.template get<Matrix>(false);
    
    for (size_t i=0; i<A.ridx.size(); ++i)
      f << A.ridx[i]+1 << ' ' << A.cidx[i]+1 << ' ' << A.data[i] << '\n';
    f << "];\nA = sparse(data(:,1),data(:,2),data(:,3)," << A.N() << "," << A.M() << ");\n";
  }
}
 
 
#endif