matrix.h

Go to the documentation of this file.

// This file is part of HermesCommon
//
// Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
// Email: hpfem-group@unr.edu, home page: http://hpfem.org/.
//
// Hermes2D 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 2 of the License,
// or (at your option) any later version.
//
// Hermes2D 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 Hermes2D; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#ifndef __HERMES_COMMON_MATRIX_H
#define __HERMES_COMMON_MATRIX_H

#include "common.h"
#include "vector.h"
#include "exceptions.h"
#include "mixins.h"

namespace Hermes
{
  enum MatrixSolverType
  {
    SOLVER_UMFPACK = 0,
    SOLVER_PETSC,
    SOLVER_MUMPS,
    SOLVER_SUPERLU,
    SOLVER_AMESOS,
    SOLVER_AZTECOO
  };

  namespace Algebra
  {
    namespace DenseMatrixOperations
    {
      template<typename T>
      T **new_matrix(unsigned int m, unsigned int n = 0)
      {
        if(!n) n = m;
        T **vec = (T **) new char[sizeof(T *) * m + sizeof(T) * m * n];
        memset(vec, 0, sizeof(T *) * m + sizeof(T) * m * n);
        T *row = (T *) (vec + m);
        for (unsigned int i = 0; i < m; i++, row += n) vec[i] = row;
        return vec;
      }

      template<typename T>
      T **new_matrix_malloc(unsigned int m, unsigned int n = 0)
      {
        if(!n) n = m;
        T **vec = (T **) malloc(sizeof(T *) * m + sizeof(T) * m * n);
        memset(vec, 0, sizeof(T *) * m + sizeof(T) * m * n);
        T *row = (T *) (vec + m);
        for (unsigned int i = 0; i < m; i++, row += n) vec[i] = row;
        return vec;
      }

      template<typename T>
      void copy_matrix(T** dest, T** src, unsigned int m, unsigned int n = 0)
      {
        if(n == 0) n = m;
        for(unsigned int i = 0; i < m; i++)
        {
          memcpy(dest[i], src[i], n*sizeof(T));
        }
      }

      template<typename T>
      void save_matrix_octave(const std::string& matrix_name, T** matrix, unsigned int m, unsigned int n = 0, const std::string& filename = std::string())
      {
        if(n == 0) n = m;

        //create filename
        std::string fname = filename;
        if(fname.empty())
          fname = matrix_name + ".mat";

        //open file
        std::ofstream fout(fname.c_str());
        if(!fout.is_open())
        {
          throw Hermes::Exceptions::Exception("Unable to save a matrix to a file \"%s\"", fname.c_str());
          return;
        }

        //write header
        fout << std::string("# name: ") << matrix_name << std::endl;
        fout << std::string("# type: matrix") << std::endl;
        fout << std::string("# rows: ") << m << std::endl;
        fout << std::string("# columns: ") << n << std::endl;

        //write contents
        for(unsigned int i = 0; i < m; i++)
        {
          for(unsigned int k = 0; k < n; k++)
            fout << ' ' << matrix[i][k];
          fout << std::endl;
        }

        //finish
        fout.close();
      }

      template<typename T>
      void save_sparse_matrix_octave(const std::string& matrix_name, const T* Ax, const int* Ap, const int* Ai,
        unsigned int m, const std::string& filename = std::string())
      {
        // create filename
        std::string fname = filename;
        if(fname.empty())
          fname = matrix_name + ".mat";

        // open file
        std::ofstream fout(fname.c_str());
        if(!fout.is_open())
        {
          throw Hermes::Exceptions::Exception("Unable to save a matrix to a file \"%s\"", fname.c_str());
          return;
        }

        // write header
        fout << std::string("# name: ") << matrix_name << std::endl;
        fout << std::string("# type: sparse matrix") << std::endl;
        fout << std::string("# nnz: ") << Ap[m] << std::endl;
        fout << std::string("# rows: ") << m << std::endl;
        fout << std::string("# columns: ") << m << std::endl;

        // write contents
        for (int j = 0; j < m; j++)
          for (int i = Ap[j]; i < Ap[j + 1]; i++)
            fout << j + 1 << " " << Ai[i] + 1 << " " << Ax[i] << std::endl;

        // finish
        fout.close();
      }

      template<typename T>
      void transpose(T **matrix, unsigned int m, unsigned int n)
      {
        unsigned int min = std::min(m, n);
        for (unsigned int i = 0; i < min; i++)
          for (unsigned int j = i + 1; j < min; j++)
            std::swap(matrix[i][j], matrix[j][i]);

        if(m < n)
          for (unsigned int i = 0; i < m; i++)
            for (unsigned int j = m; j < n; j++)
              matrix[j][i] = matrix[i][j];
        else if(n < m)
          for (unsigned int i = n; i < m; i++)
            for (unsigned int j = 0; j < n; j++)
              matrix[j][i] = matrix[i][j];
      }

      template<typename T>
      void chsgn(T **matrix, unsigned int m, unsigned int n)
      {
        for (unsigned int i = 0; i < m; i++)
          for (unsigned int j = 0; j < n; j++)
            matrix[i][j] = -matrix[i][j];
      }

      HERMES_API void ludcmp(double **a, int n, int *indx, double *d);

      template<typename T>
      void lubksb(double **a, int n, int *indx, T *b)
      {
        int i, ip, j;
        T sum;

        for (i = 0; i < n; i++)
        {
          ip = indx[i];
          sum = b[ip];
          b[ip] = b[i];
          for (j = 0; j < i; j++) sum -= a[i][j]*b[j];
          b[i] = sum;
        }
        for (i = n-1; i >= 0; i--)
        {
          sum = b[i];
          for (j = i + 1; j < n; j++) sum -= a[i][j]*b[j];
          b[i] = sum / a[i][i];
        }
      }

      HERMES_API void choldc(double **a, int n, double p[]);

      template<typename T>
      void cholsl(double **a, int n, double p[], T b[], T x[])
      {
        int i, k;
        T sum;

        for (i = 0; i < n; i++)
        {
          sum = b[i];
          k = i;
          while (--k >= 0) sum -= a[i][k] * x[k];
          x[i] = sum / p[i];
        }

        for (i = n-1; i >= 0; i--)
        {
          sum = x[i];
          k = i;
          while (++k < n) sum -= a[k][i] * x[k];
          x[i] = sum / p[i];
        }
      }
    }

    enum EMatrixDumpFormat
    {
      DF_MATLAB_SPARSE, 




      DF_PLAIN_ASCII,
      DF_HERMES_BIN,
      DF_MATRIX_MARKET 
    };

    template<typename Scalar>
    class HERMES_API Matrix : public Hermes::Mixins::Loggable
    {
    public:
      unsigned int get_size() { return this->size;};

      Matrix(unsigned int size) { this->size = size;};

      virtual ~Matrix() {};

      Matrix() { this->size = 0;};

      virtual void alloc() = 0;

      virtual void free() = 0;

      virtual Scalar get(unsigned int m, unsigned int n) = 0;

      virtual void zero() = 0;

      virtual void add_to_diagonal(Scalar v) = 0;

      virtual void add(unsigned int m, unsigned int n, Scalar v) = 0;

      virtual void add(unsigned int m, unsigned int n, Scalar **mat, int *rows, int *cols) = 0;

      virtual bool dump(FILE *file, const char *var_name, EMatrixDumpFormat fmt = DF_MATLAB_SPARSE, char* number_format = "%lf") = 0;

      virtual unsigned int get_matrix_size() const = 0;

    protected:

      unsigned int size;  
    };

    template<typename Scalar>
    class HERMES_API SparseMatrix : public Matrix<Scalar> {
    public:
      SparseMatrix();
      SparseMatrix(unsigned int size);
      virtual ~SparseMatrix();

      virtual void prealloc(unsigned int n);

      virtual void pre_add_ij(unsigned int row, unsigned int col);

      virtual void finish() { }

      virtual unsigned int get_size() { return this->size; }

      virtual void add_sparse_matrix(SparseMatrix* mat)
      {
        throw Hermes::Exceptions::Exception("add_sparse_matrix() undefined.");
      };

      virtual void add_sparse_to_diagonal_blocks(int num_stages, SparseMatrix<Scalar>* mat)
      {
        throw Hermes::Exceptions::Exception("add_sparse_to_diagonal_blocks() undefined.");
      };

      virtual int get_num_row_entries(unsigned int row) { return -1; }

      virtual void extract_row_copy(unsigned int row, unsigned int len,
        unsigned int &n_entries, double *vals,
        unsigned int *idxs) { }

      virtual int get_num_col_entries(unsigned int col) { return -1; }

      virtual void extract_col_copy(unsigned int col, unsigned int len,
        unsigned int &n_entries, double *vals,
        unsigned int *idxs) { }

      virtual void multiply_with_vector(Scalar* vector_in, Scalar* vector_out) {
        throw Hermes::Exceptions::Exception("multiply_with_vector() undefined.");
      };

      virtual void multiply_with_Scalar(Scalar value) {
        throw Hermes::Exceptions::Exception("multiply_with_Scalar() undefined.");
      };

      virtual SparseMatrix* duplicate() { return (SparseMatrix*)NULL;};

      virtual double get_fill_in() const = 0;

      unsigned row_storage:1; 
      unsigned col_storage:1; 

      virtual unsigned int get_nnz() const {
        throw Hermes::Exceptions::Exception("get_nnz() undefined.");
        return 0;
      }

    protected:
      static const int PAGE_SIZE = 62;

      struct Page {
        int count;
        int idx[PAGE_SIZE];
        Page *next;
      };

      Page **pages;

      int sort_and_store_indices(Page *page, int *buffer, int *max);
      int get_num_indices();

      int mem_size;
    };

    template<typename Scalar>
    class HERMES_API Vector : public Hermes::Mixins::Loggable
    {
    public:
      virtual ~Vector() { }

      virtual void alloc(unsigned int ndofs) = 0;
      virtual void free() = 0;
      virtual void finish() { }

      virtual Scalar get(unsigned int idx) = 0;

      virtual void extract(Scalar *v) const = 0;

      virtual void zero() = 0;

      virtual void change_sign() = 0;

      virtual void set(unsigned int idx, Scalar y) = 0;

      virtual void add(unsigned int idx, Scalar y) = 0;

      virtual void add_vector(Vector<Scalar>* vec) = 0;
      virtual void add_vector(Scalar* vec) = 0;

      virtual void add(unsigned int n, unsigned int *idx, Scalar *y) = 0;

      unsigned int length() const {return this->size;}

      virtual bool dump(FILE *file, const char *var_name,
        EMatrixDumpFormat fmt = DF_MATLAB_SPARSE, char* number_format = "%lf") = 0;

    protected:
      unsigned int size;
    };

    template<typename Scalar> HERMES_API
      Vector<Scalar>* create_vector();

    template<typename Scalar> HERMES_API
      SparseMatrix<Scalar>*  create_matrix();
  }
}
#endif