dense_matrix_operations.cpp

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.
#include "dense_matrix_operations.h"
#include "exceptions.h"
#include "util/memory_handling.h"

namespace Hermes
{
  namespace Algebra
  {
    namespace DenseMatrixOperations
    {
      template<typename T, typename Int>
      void ludcmp(T **a, Int n, Int *indx, double *d)
      {
        Int i, imax = 0, j, k;
        T big, dum, sum, temp;
        T *vv = malloc_with_check<T>(n);

        *d = 1.0;
        for (i = 0; i < n; i++)
        {
          big = 0.0;
          for (j = 0; j < n; j++)
          {
            temp = a[i][j];
            if (std::abs(temp) > std::abs(big))
              big = temp;
          }
          if (big == 0.0)
          {
            free_with_check(vv);
            throw Exceptions::Exception("Singular matrix in routine LUDCMP!");
          }
          vv[i] = 1.0 / big;
        }
        for (j = 0; j < n; j++)
        {
          for (i = 0; i < j; i++)
          {
            sum = a[i][j];
            for (k = 0; k < i; k++) sum -= a[i][k] * a[k][j];
            a[i][j] = sum;
          }
          big = 0.0;
          for (i = j; i < n; i++)
          {
            sum = a[i][j];
            for (k = 0; k < j; k++) sum -= a[i][k] * a[k][j];
            a[i][j] = sum;
            dum = vv[i] * std::abs(sum);
            if (std::abs(dum) >= std::abs(big))
            {
              big = dum;
              imax = i;
            }
          }
          if (j != imax)
          {
            for (k = 0; k < n; k++)
            {
              dum = a[imax][k];
              a[imax][k] = a[j][k];
              a[j][k] = dum;
            }
            *d = -(*d);
            vv[imax] = vv[j];
          }
          indx[j] = imax;
          if (a[j][j] == 0.0) a[j][j] = 1.0e-20;
          if (j != n - 1)
          {
            dum = 1.0 / (a[j][j]);
            for (i = j + 1; i < n; i++) a[i][j] *= dum;
          }
        }
        free_with_check(vv);
      }

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

        for (i = 0; i < n; i++)
        {
          ip = indx[i];
          sum = b[ip];
          b[ip] = b[i];
          for (j = 0; j < i; j++)
            sum = 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];
        }
      }

      template<typename T, typename Int>
      void choldc(T **a, Int n, T p[])
      {
        Int i, j, k;
        for (i = 0; i < n; i++)
        {
          for (j = i; j < n; j++)
          {
            T sum = a[i][j];
            k = i;
            while (k > 0)
            {
              sum -= a[i][k - 1] * a[j][k - 1];
              k--;
            }
            if (i == j)
            {
              if ((std::complex<double>(sum)).real() <= 0.0)
                throw Exceptions::Exception("CHOLDC failed!");
              else p[i] = sqrt(sum);
            }
            else a[j][i] = sum / p[i];
          }
        }
      }

      template HERMES_API void ludcmp<double, int>(double **a, int n, int *indx, double *d);
      template HERMES_API void ludcmp<std::complex<double>, int>(std::complex<double> **a, int n, int *indx, double *d);

      template HERMES_API void ludcmp<double, unsigned short>(double **a, unsigned short n, unsigned short *indx, double *d);
      template HERMES_API void ludcmp<std::complex<double>, unsigned short>(std::complex<double> **a, unsigned short n, unsigned short *indx, double *d);

      template HERMES_API void ludcmp<double, unsigned char>(double **a, unsigned char n, unsigned char *indx, double *d);
      template HERMES_API void ludcmp<std::complex<double>, unsigned char>(std::complex<double> **a, unsigned char n, unsigned char *indx, double *d);

      template HERMES_API void lubksb<double, std::complex<double>, int>(double **a, int n, int *indx, std::complex<double> *d);
      template HERMES_API void lubksb<double, double, int>(double **a, int n, int *indx, double *d);
      template HERMES_API void lubksb<std::complex<double>, std::complex<double>, int>(std::complex<double> **a, int n, int *indx, std::complex<double> *d);

      template HERMES_API void lubksb<double, std::complex<double>, unsigned short>(double **a, unsigned short n, unsigned short *indx, std::complex<double> *d);
      template HERMES_API void lubksb<double, double, unsigned short>(double **a, unsigned short n, unsigned short *indx, double *d);
      template HERMES_API void lubksb<std::complex<double>, std::complex<double>, unsigned short>(std::complex<double> **a, unsigned short n, unsigned short *indx, std::complex<double> *d);

      template HERMES_API void lubksb<double, std::complex<double>, unsigned char>(double **a, unsigned char n, unsigned char *indx, std::complex<double> *d);
      template HERMES_API void lubksb<double, double, unsigned char>(double **a, unsigned char n, unsigned char *indx, double *d);
      template HERMES_API void lubksb<std::complex<double>, std::complex<double>, unsigned char>(std::complex<double> **a, unsigned char n, unsigned char *indx, std::complex<double> *d);

      template HERMES_API void choldc<double, int>(double **a, int n, double p[]);
      template HERMES_API void choldc<std::complex<double>, int>(std::complex<double> **a, int n, std::complex<double> p[]);

      template HERMES_API void choldc<double, unsigned short>(double **a, unsigned short n, double p[]);
      template HERMES_API void choldc<std::complex<double>, unsigned short>(std::complex<double> **a, unsigned short n, std::complex<double> p[]);

      template HERMES_API void choldc<double, unsigned char>(double **a, unsigned char n, double p[]);
      template HERMES_API void choldc<std::complex<double>, unsigned char>(std::complex<double> **a, unsigned char n, std::complex<double> p[]);
    }
  }
}