DirectSolver Class Reference

#include <DirectSolver.h>

Inheritance diagram for DirectSolver:

Collaboration diagram for DirectSolver:

Public Member Functions
	DirectSolver ()
virtual	~DirectSolver ()
virtual void	cleanup ()
virtual MFRPtr	solve (LinearSystem &ls)
virtual void	smooth (LinearSystem &ls)
Public Member Functions inherited from LinearSolver
	LinearSolver ()

Private Member Functions
	DirectSolver (const DirectSolver &)

Additional Inherited Members
Public Attributes inherited from LinearSolver
int	nMaxIterations
int	verbosity
double	relativeTolerance
double	absoluteTolerance

Detailed Description

Solve a linear system using UMFPack

Definition at line 22 of file DirectSolver.h.

Constructor & Destructor Documentation

DirectSolver::DirectSolver

(

)

Definition at line 11 of file DirectSolver.cpp.

{}

DirectSolver::~DirectSolver ( )

virtual

Definition at line 14 of file DirectSolver.cpp.

{}

DirectSolver::DirectSolver ( const DirectSolver &  )

private

Member Function Documentation

void DirectSolver::cleanup ( )

virtual

Implements LinearSolver.

Definition at line 18 of file DirectSolver.cpp.

{}

void DirectSolver::smooth ( LinearSystem &  ls )

virtual

Implements LinearSolver.

Definition at line 108 of file DirectSolver.cpp.

{
  throw CException("cannot use DirectSolver as preconditioner");
}

MFRPtr DirectSolver::solve ( LinearSystem &  ls )

virtual

Implements LinearSolver.

Definition at line 23 of file DirectSolver.cpp.

References MultiFieldMatrix::computeResidual(), MultiField::getArrayIndex(), LinearSystem::getB(), Matrix::getConnectivity(), ArrayBase::getData(), Array< T >::getData(), ArrayBase::getDataSize(), LinearSystem::getDelta(), MultiField::getLength(), ArrayBase::getLength(), LinearSystem::getMatrix(), MultiFieldMatrix::getMatrix(), CRConnectivity::getMultiTranspose(), MultiField::getOneNorm(), LinearSystem::getResidual(), Matrix::setFlatMatrix(), and LinearSolver::verbosity.

{

  typedef CRMatrix<double,double,double>  FlatMatrix;
  
  const MultiFieldMatrix& mfMatrix = ls.getMatrix();
  const MultiField& bField = ls.getB();
  MultiField& deltaField = ls.getDelta();

  if (bField.getLength() != 1)
    throw CException("direct solver only works for single matrix");

  const MultiField::ArrayIndex rowIndex = bField.getArrayIndex(0);

  const Matrix& origMatrix = mfMatrix.getMatrix(rowIndex,rowIndex);
  const CRConnectivity& origConn = origMatrix.getConnectivity();

  const ArrayBase& origB = bField[rowIndex];
  const int blockSize = origB.getDataSize()/(sizeof(double)*origB.getLength());

  shared_ptr<CRConnectivity> flatConnPtr = origConn.getMultiTranspose(blockSize);

  FlatMatrix flatMatrix(*flatConnPtr);

  origMatrix.setFlatMatrix(flatMatrix);

  Array<double>& flatCoeffs = flatMatrix.getOffDiag();
  const Array<int>& flatRow = flatConnPtr->getRow();
  const Array<int>& flatCol = flatConnPtr->getCol();
  const ArrayBase& flatB = origB;
  
  ArrayBase& flatDelta = deltaField[rowIndex];
  
  const int nFlatEqs = flatConnPtr->getRowDim();
    
  // original system is in delta form

  mfMatrix.computeResidual(deltaField,bField,ls.getResidual());

  MFRPtr rNorm0(ls.getResidual().getOneNorm());

  if (verbosity >0)
    cout << 0 << ": " << *rNorm0 << endl;
  
  
  double *null = (double *) NULL ;
  
  void *Symbolic, *Numeric ;
  (void) umfpack_di_symbolic (nFlatEqs, nFlatEqs,
                              (int*)flatRow.getData(),
                              (int*)flatCol.getData(),
                              (double*) flatCoeffs.getData(),
                              &Symbolic, null, null) ;
  (void) umfpack_di_numeric ( (int*)flatRow.getData(),
                              (int*)flatCol.getData(),
                              (double*) flatCoeffs.getData(),
                              
                              Symbolic, &Numeric, null, null) ;
  umfpack_di_free_symbolic (&Symbolic) ;
  (void) umfpack_di_solve (UMFPACK_A,
                           (int*)flatRow.getData(),
                           (int*)flatCol.getData(),
                           (double*) flatCoeffs.getData(),
                           (double*) flatDelta.getData(),
                           (double*) flatB.getData(),
                           Numeric, null, null) ;
   umfpack_di_free_numeric (&Numeric) ;

   // umfpack solves ax=b, we want ax+b=0;
   double *flatDeltaData = (double*) flatDelta.getData();
   for(int n=0; n<nFlatEqs; n++)
     flatDeltaData[n] *= -1.0;
   
  mfMatrix.computeResidual(deltaField,bField,ls.getResidual());

  MFRPtr rNormN(ls.getResidual().getOneNorm());

  if (verbosity >0)
    cout <<  "Final : " << *rNormN << endl;
  

  return rNorm0;
}

---

The documentation for this class was generated from the following files:

• DirectSolver.h
• DirectSolver.cpp