densemat.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#ifndef MFEM_DENSEMAT
#define MFEM_DENSEMAT

#include "../config/config.hpp"
#include "matrix.hpp"

namespace mfem
{

class DenseMatrix : public Matrix
{
   friend class DenseTensor;
   friend class DenseMatrixInverse;

private:
   double *data;
   int capacity; // zero or negative capacity means we do not own the data.

   void Eigensystem(Vector &ev, DenseMatrix *evect = NULL);

public:
   DenseMatrix();

   DenseMatrix(const DenseMatrix &);

   explicit DenseMatrix(int s);

   DenseMatrix(int m, int n);

   DenseMatrix(const DenseMatrix &mat, char ch);

   DenseMatrix(double *d, int h, int w) : Matrix(h, w)
   { data = d; capacity = -h*w; }

   void UseExternalData(double *d, int h, int w)
   { data = d; height = h; width = w; capacity = -h*w; }

   void ClearExternalData() { data = NULL; height = width = 0; capacity = 0; }

   int Size() const { return Width(); }

   void SetSize(int s) { SetSize(s, s); }

   void SetSize(int h, int w);

   inline double *Data() const { return data; }

   inline double &operator()(int i, int j);

   inline const double &operator()(int i, int j) const;

   double operator*(const DenseMatrix &m) const;

   double Trace() const;

   virtual double &Elem(int i, int j);

   virtual const double &Elem(int i, int j) const;

   void Mult(const double *x, double *y) const;

   virtual void Mult(const Vector &x, Vector &y) const;

   void MultTranspose(const double *x, double *y) const;

   virtual void MultTranspose(const Vector &x, Vector &y) const;

   void AddMult(const Vector &x, Vector &y) const;

   void AddMult_a(double a, const Vector &x, Vector &y) const;

   // y += a * A^t x
   void AddMultTranspose_a(double a, const Vector &x, Vector &y) const;

   double InnerProduct(const double *x, const double *y) const;

   void LeftScaling(const Vector & s);
   void InvLeftScaling(const Vector & s);
   void RightScaling(const Vector & s);
   void InvRightScaling(const Vector & s);
   void SymmetricScaling(const Vector & s);
   void InvSymmetricScaling(const Vector & s);

   double InnerProduct(const Vector &x, const Vector &y) const
   { return InnerProduct((const double *)x, (const double *)y); }

   virtual MatrixInverse *Inverse() const;

   void Invert();

   double Det() const;

   double Weight() const;

   void Add(const double c, const DenseMatrix &A);

   DenseMatrix &operator=(double c);

   DenseMatrix &operator=(const double *d);

   DenseMatrix &operator=(const DenseMatrix &m);

   DenseMatrix &operator+=(DenseMatrix &m);

   DenseMatrix &operator-=(DenseMatrix &m);

   DenseMatrix &operator*=(double c);

   void Neg();

   void Norm2(double *v) const;

   double MaxMaxNorm() const;

   double FNorm() const;

   void Eigenvalues(Vector &ev)
   { Eigensystem(ev); }

   void Eigenvalues(Vector &ev, DenseMatrix &evect)
   { Eigensystem(ev, &evect); }

   void Eigensystem(Vector &ev, DenseMatrix &evect)
   { Eigensystem(ev, &evect); }

   void SingularValues(Vector &sv) const;
   int Rank(double tol) const;

   double CalcSingularvalue(const int i) const;

   void CalcEigenvalues(double *lambda, double *vec) const;

   void GetColumn(int c, Vector &col);

   void GetColumnReference(int c, Vector &col)
   { col.SetDataAndSize(data + c * height, height); }

   void GetDiag(Vector &d);
   void Getl1Diag(Vector &l);
   void GetRowSums(Vector &l) const;

   void Diag(double c, int n);
   void Diag(double *diag, int n);

   void Transpose();
   void Transpose(DenseMatrix &A);
   void Symmetrize();

   void Lump();

   void GradToCurl(DenseMatrix &curl);
   void GradToDiv(Vector &div);

   void CopyRows(const DenseMatrix &A, int row1, int row2);
   void CopyCols(const DenseMatrix &A, int col1, int col2);
   void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco);
   void CopyMN(const DenseMatrix &A, int row_offset, int col_offset);
   void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset);
   void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco,
               int row_offset, int col_offset);
   void CopyMNDiag(double c, int n, int row_offset, int col_offset);
   void CopyMNDiag(double *diag, int n, int row_offset, int col_offset);

   void AddMatrix(DenseMatrix &A, int ro, int co);
   void AddMatrix(double a, DenseMatrix &A, int ro, int co);

   void AddToVector(int offset, Vector &v) const;
   void GetFromVector(int offset, const Vector &v);
   void AdjustDofDirection(Array<int> &dofs);

   void SetRow(int row, double value);
   void SetCol(int col, double value);

   void Threshold(double eps);

   int CheckFinite() const { return mfem::CheckFinite(data, height*width); }

   virtual void Print(std::ostream &out = std::cout, int width_ = 4) const;
   virtual void PrintMatlab(std::ostream &out = std::cout) const;
   virtual void PrintT(std::ostream &out = std::cout, int width_ = 4) const;

   void TestInversion();

   virtual ~DenseMatrix();
};

void Add(const DenseMatrix &A, const DenseMatrix &B,
         double alpha, DenseMatrix &C);

void Add(double alpha, const DenseMatrix &A,
         double beta,  const DenseMatrix &B, DenseMatrix &C);

void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);

void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);

void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja);

void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat);

void CalcInverse(const DenseMatrix &a, DenseMatrix &inva);

void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva);

void CalcOrtho(const DenseMatrix &J, Vector &n);

void MultAAt(const DenseMatrix &a, DenseMatrix &aat);

void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt);

void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt);

void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);

void MultADBt(const DenseMatrix &A, const Vector &D,
              const DenseMatrix &B, DenseMatrix &ADAt);

void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);

void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB);

void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);

void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);

void MultVVt(const Vector &v, DenseMatrix &vvt);

void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);

void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);

void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
                   DenseMatrix &VWt);

void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);


class LUFactors
{
public:
   double *data;
   int *ipiv;
#ifdef MFEM_USE_LAPACK
   static const int ipiv_base = 1;
#else
   static const int ipiv_base = 0;
#endif

   LUFactors() { }

   LUFactors(double *data_, int *ipiv_) : data(data_), ipiv(ipiv_) { }

   void Factor(int m);

   void Mult(int m, int n, double *X) const;

   void LSolve(int m, int n, double *X) const;

   void USolve(int m, int n, double *X) const;

   void Solve(int m, int n, double *X) const;

   void GetInverseMatrix(int m, double *X) const;

   static void SubMult(int m, int n, int r, const double *A21,
                       const double *X1, double *X2);

   void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const;

   void BlockForwSolve(int m, int n, int r, const double *L21,
                       double *B1, double *B2) const;

   void BlockBackSolve(int m, int n, int r, const double *U12,
                       const double *X2, double *Y1) const;
};


class DenseMatrixInverse : public MatrixInverse
{
private:
   const DenseMatrix *a;
   LUFactors lu;

public:
   DenseMatrixInverse() : a(NULL), lu(NULL, NULL) { }

   DenseMatrixInverse(const DenseMatrix &mat);

   DenseMatrixInverse(const DenseMatrix *mat);

   int Size() const { return Width(); }

   void Factor();

   void Factor(const DenseMatrix &mat);

   virtual void SetOperator(const Operator &op);

   virtual void Mult(const Vector &x, Vector &y) const;

   void Mult(const DenseMatrix &B, DenseMatrix &X) const;

   void GetInverseMatrix(DenseMatrix &Ainv) const
   {
      Ainv.SetSize(width);
      lu.GetInverseMatrix(width, Ainv.Data());
   }

   void TestInversion();

   virtual ~DenseMatrixInverse();
};


class DenseMatrixEigensystem
{
   DenseMatrix &mat;
   Vector      EVal;
   DenseMatrix EVect;
   Vector ev;
   int n;

#ifdef MFEM_USE_LAPACK
   double *work;
   char jobz, uplo;
   int lwork, info;
#endif

public:

   DenseMatrixEigensystem(DenseMatrix &m);
   void Eval();
   Vector &Eigenvalues() { return EVal; }
   DenseMatrix &Eigenvectors() { return EVect; }
   double Eigenvalue(int i) { return EVal(i); }
   const Vector &Eigenvector(int i)
   {
      ev.SetData(EVect.Data() + i * EVect.Height());
      return ev;
   }
   ~DenseMatrixEigensystem();
};


class DenseMatrixSVD
{
   Vector sv;
   int m, n;

#ifdef MFEM_USE_LAPACK
   double *work;
   char jobu, jobvt;
   int lwork, info;
#endif

   void Init();
public:

   DenseMatrixSVD(DenseMatrix &M);
   DenseMatrixSVD(int h, int w);
   void Eval(DenseMatrix &M);
   Vector &Singularvalues() { return sv; }
   double Singularvalue(int i) { return sv(i); }
   ~DenseMatrixSVD();
};

class Table;

class DenseTensor
{
private:
   DenseMatrix Mk;
   double *tdata;
   int nk;

public:
   DenseTensor() { nk = 0; tdata = NULL; }

   DenseTensor(int i, int j, int k)
      : Mk(NULL, i, j)
   { nk = k; tdata = new double[i*j*k]; }

   int SizeI() const { return Mk.Height(); }
   int SizeJ() const { return Mk.Width(); }
   int SizeK() const { return nk; }

   void SetSize(int i, int j, int k)
   {
      delete [] tdata;
      Mk.UseExternalData(NULL, i, j);
      nk = k;
      tdata = new double[i*j*k];
   }

   DenseMatrix &operator()(int k) { Mk.data = GetData(k); return Mk; }
   const DenseMatrix &operator()(int k) const
   { return const_cast<DenseTensor&>(*this)(k); }

   double &operator()(int i, int j, int k)
   { return tdata[i+SizeI()*(j+SizeJ()*k)]; }
   const double &operator()(int i, int j, int k) const
   { return tdata[i+SizeI()*(j+SizeJ()*k)]; }

   double *GetData(int k) { return tdata+k*Mk.Height()*Mk.Width(); }

   double *Data() { return tdata; }

   void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const;

   ~DenseTensor() { delete [] tdata; }
};


// Inline methods

inline double &DenseMatrix::operator()(int i, int j)
{
#ifdef MFEM_DEBUG
   if ( data == 0 || i < 0 || i >= height || j < 0 || j >= width )
   {
      mfem_error("DenseMatrix::operator()");
   }
#endif

   return data[i+j*height];
}

inline const double &DenseMatrix::operator()(int i, int j) const
{
#ifdef MFEM_DEBUG
   if ( data == 0 || i < 0 || i >= height || j < 0 || j >= width )
   {
      mfem_error("DenseMatrix::operator() const");
   }
#endif

   return data[i+j*height];
}

}

#endif