dense_matrix_operations.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://www.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_DENSE_MATRIX_OPERATIONS_H
#define __HERMES_COMMON_DENSE_MATRIX_OPERATIONS_H

#include "common.h"
#include "util/compat.h"
#include "util/memory_handling.h"
#include "exceptions.h"

namespace Hermes
{
  namespace Algebra
  {
    namespace DenseMatrixOperations
    {
      template<typename T>
      T **new_matrix(unsigned int m, unsigned int n = 0)
      {
        if (!n) n = m;
        T **vec = (T**)calloc_with_check_direct_size<char>(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 transpose(T *matrix, unsigned int m, unsigned int n, unsigned int size)
      {
        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 * size + j], matrix[j * size + i]);
          }
        }

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

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

      template<typename T>
      void change_sign(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];
      }
      template<typename T>
      void change_sign(T *matrix, unsigned int m, unsigned int n, unsigned int size)
      {
        for (unsigned int i = 0; i < m; i++)
          for (unsigned int j = 0; j < n; j++)
          {
          int local_matrix_index_array = i * size + j;
          matrix[local_matrix_index_array] = -matrix[local_matrix_index_array];
          }
      }

      template<typename T, typename Int>
      HERMES_API void ludcmp(T **a, Int n, Int *indx, double *d);

      template<typename T, typename S, typename Int>
      HERMES_API void lubksb(T **a, Int n, Int *indx, S *b);

      template<typename T, typename Int>
      HERMES_API void choldc(T **a, Int n, T p[]);

      template<typename T, typename Int>
      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];
        }
      }
    }
  }
}
#endif