sparsemat.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.googlecode.com.
//
// 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_SPARSEMAT
#define MFEM_SPARSEMAT

// Data types for sparse matrix

#include "../general/mem_alloc.hpp"
#include "../general/table.hpp"
#include "densemat.hpp"

class RowNode
{
public:
   RowNode *Prev;
   int Column;
   double Value;
};

class SparseMatrix : public Matrix
{
private:

   int *I, *J, width;

   double *A;

   RowNode **Rows;

   int current_row;
   union { int *J; RowNode **Node; } ColPtr;

#ifdef MFEM_USE_MEMALLOC
   MemAlloc <RowNode, 1024> NodesMem;
#endif

   inline void SetColPtr(const int row);
   inline void ClearColPtr();
   inline double &SearchRow(const int col);
   inline void _Add_(const int col, const double a)
   { SearchRow(col) += a; }
   inline void _Set_(const int col, const double a)
   { SearchRow(col) = a; }
   inline double _Get_(const int col);

   inline double &SearchRow(const int row, const int col);
   inline void _Add_(const int row, const int col, const double a)
   { SearchRow(row, col) += a; }
   inline void _Set_(const int row, const int col, const double a)
   { SearchRow(row, col) = a; }

public:
   explicit SparseMatrix(int nrows, int ncols = 0);

   SparseMatrix(int *i, int *j, double *data, int m, int n)
      : Matrix (m), I(i), J(j), width(n), A(data)
   { Rows = NULL; ColPtr.J = NULL; }

   inline int *GetI() const { return I; }
   inline int *GetJ() const { return J; }
   inline double *GetData() const { return A; }
   inline int Width() const { return width; }
   int RowSize(int i);

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

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

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

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

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

   void AddMult(const Vector &x, Vector &y, const double a = 1.0) const;

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

   void AddMultTranspose(const Vector &x, Vector &y,
                         const double a = 1.0) const;

   void PartMult(const Array<int> &rows, const Vector &x, Vector &y);

   double InnerProduct(const Vector &x, const Vector &y) const;

   void GetRowSums(Vector &x) const;

   virtual MatrixInverse *Inverse() const;

   void EliminateRow(int row, const double sol, Vector &rhs);
   void EliminateRow(int row);
   void EliminateCol(int col);
   void EliminateCols(Array<int> &cols, Vector *x = NULL, Vector *b = NULL);

   void EliminateRowCol(int rc, const double sol, Vector &rhs, int d = 0);

   void EliminateRowColMultipleRHS(int rc, const Vector &sol,
                                   DenseMatrix &rhs, int d = 0);

   void EliminateRowCol(int rc, int d = 0);
   // Same as above + save the eliminated entries in Ae so that
   // (*this) + Ae is the original matrix
   void EliminateRowCol(int rc, SparseMatrix &Ae, int d = 0);

   void SetDiagIdentity();
   void EliminateZeroRows();

   void Gauss_Seidel_forw(const Vector &x, Vector &y) const;
   void Gauss_Seidel_back(const Vector &x, Vector &y) const;

   double GetJacobiScaling() const;
   void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const;

   void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
                double sc = 1.0) const;

   virtual void Finalize(int skip_zeros = 1);

   int Finalized() { return (A != NULL); }

   void GetBlocks(Array2D<SparseMatrix *> &blocks) const;

   void GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
                     DenseMatrix &subm);

   void Set(const int i, const int j, const double a);
   void Add(const int i, const int j, const double a);

   void SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
                     const DenseMatrix &subm, int skip_zeros = 1);

   void SetSubMatrixTranspose(const Array<int> &rows, const Array<int> &cols,
                              const DenseMatrix &subm, int skip_zeros = 1);

   void AddSubMatrix(const Array<int> &rows, const Array<int> &cols,
                     const DenseMatrix &subm, int skip_zeros = 1);

   void AddRow(const int row, const Array<int> &cols, const Vector &srow);

   void ScaleRow(const int row, const double scale);

   SparseMatrix &operator+=(SparseMatrix &B);

   SparseMatrix &operator=(double a);

   void Print(ostream &out = cout, int width = 4) const;

   void PrintMatlab(ostream &out = cout) const;

   void PrintMM(ostream &out = cout) const;

   void PrintCSR(ostream &out) const;

   void PrintCSR2(ostream &out) const;

   int Walk(int &i, int &j, double &a);

   double IsSymmetric() const;

   void Symmetrize();

   int NumNonZeroElems() const;

   double MaxNorm() const;

   int CountSmallElems(double tol);

   void LoseData() { I=0; J=0; A=0; }

   virtual ~SparseMatrix();
};

void SparseMatrixFunction(SparseMatrix &S, double (*f)(double));


SparseMatrix *Transpose(SparseMatrix &A);

SparseMatrix *Mult(SparseMatrix &A, SparseMatrix &B,
                   SparseMatrix *OAB = NULL);


SparseMatrix *RAP(SparseMatrix &A, SparseMatrix &R,
                  SparseMatrix *ORAP = NULL);

SparseMatrix *Mult_AtDA(SparseMatrix &A, Vector &D,
                        SparseMatrix *OAtDA = NULL);


// Inline methods

inline void SparseMatrix::SetColPtr(const int row)
{
   if (Rows)
   {
      if (ColPtr.Node == NULL)
      {
         ColPtr.Node = new RowNode *[width];
         for (int i = 0; i < width; i++)
            ColPtr.Node[i] = NULL;
      }
      for (RowNode *node_p = Rows[row]; node_p != NULL; node_p = node_p->Prev)
      {
         ColPtr.Node[node_p->Column] = node_p;
      }
   }
   else
   {
      if (ColPtr.J == NULL)
      {
         ColPtr.J = new int[width];
         for (int i = 0; i < width; i++)
            ColPtr.J[i] = -1;
      }
      for (int j = I[row], end = I[row+1]; j < end; j++)
      {
         ColPtr.J[J[j]] = j;
      }
   }
   current_row = row;
}

inline void SparseMatrix::ClearColPtr()
{
   if (Rows)
      for (RowNode *node_p = Rows[current_row]; node_p != NULL;
           node_p = node_p->Prev)
      {
         ColPtr.Node[node_p->Column] = NULL;
      }
   else
      for (int j = I[current_row], end = I[current_row+1]; j < end; j++)
      {
         ColPtr.J[J[j]] = -1;
      }
}

inline double &SparseMatrix::SearchRow(const int col)
{
   if (Rows)
   {
      RowNode *node_p = ColPtr.Node[col];
      if (node_p == NULL)
      {
#ifdef MFEM_USE_MEMALLOC
         node_p = NodesMem.Alloc();
#else
         node_p = new RowNode;
#endif
         node_p->Prev = Rows[current_row];
         node_p->Column = col;
         node_p->Value = 0.0;
         Rows[current_row] = ColPtr.Node[col] = node_p;
      }
      return node_p->Value;
   }
   else
   {
      const int j = ColPtr.J[col];
      if (j == -1)
         mfem_error("SparseMatrix::SearchRow : entry is not allocated!");
      return A[j];
   }
}

inline double SparseMatrix::_Get_(const int col)
{
   if (Rows)
   {
      RowNode *node_p = ColPtr.Node[col];
      return (node_p == NULL) ? 0.0 : node_p->Value;
   }
   else
   {
      const int j = ColPtr.J[col];
      return (j == -1) ? 0.0 : A[j];
   }
}

inline double &SparseMatrix::SearchRow(const int row, const int col)
{
   if (Rows)
   {
      RowNode *node_p;

      for (node_p = Rows[row]; 1; node_p = node_p->Prev)
         if (node_p == NULL)
         {
#ifdef MFEM_USE_MEMALLOC
            node_p = NodesMem.Alloc();
#else
            node_p = new RowNode;
#endif
            node_p->Prev = Rows[row];
            node_p->Column = col;
            node_p->Value = 0.0;
            Rows[row] = node_p;
            break;
         }
         else if (node_p->Column == col)
         {
            break;
         }
      return node_p->Value;
   }
   else
   {
      int *Ip = I+row, *Jp = J;
      for (int k = Ip[0], end = Ip[1]; k < end; k++)
         if (Jp[k] == col)
            return A[k];
      mfem_error("SparseMatrix::SearchRow(...)");
   }
   return A[0];
}

#endif