Factorization

Table of Contents

• Procedure for factorization
• Reuse symbolic factorization
• Triangular solve and accuracy check

Procedure for factorization

Before the selected inversion step, the matrix saved in DistSparseMatrix format must first be factorized. In principle this can be done with any \(LDL^T\) factorization or \(LU\) factorization routines. In the current version of PEXSI, SuperLU_DIST v3.3 is used for the \(LU\) factorization.

Note

To avoid conflict with other routines in PEXSI, the SuperLU_DIST routines are encapsulated in superlu_dist_interf.cpp. Access to SuperLU_DIST routines are made through a wrapper class SuperLUMatrix.

The basic steps for factorization include:

• Convert a DistSparseMatrix into the native format (SuperMatrix in SuperLU_DIST) of the factorization routine.
• Symbolic factorization.
• Numerical factorization.

Related structures and subroutines

SuperLUGrid

A thin interface for the mpi grid strucutre in SuperLU.

SuperLUOptions

A thin interface for passing parameters to set the SuperLU options.

SuperLUMatrix::DistSparseMatrixToSuperMatrixNRloc

Convert a distributed sparse matrix in compressed sparse column format into the SuperLU compressed row format. The output is saved in the current SuperLUMatrix.

Note

Although LU factorization is used, the routine assumes that the matrix is strictly symmetric, and therefore the compressed sparse row (CSR) format, used by SuperLU_DIST, gives exactly the same matrix as formed by the compresed sparse column format (CSC).

SuperLUMatrix::SymbolicFactorize

This routine factorizes the superlu matrix symbolically. Symbolic factorization contains three steps.

• Permute the matrix to reduce fill-in.
• Symbolic factorize the matrix.
• Distribute the matrix into 2D block cyclic format.

This routine is controlled via SuperLUOptions. In particular, the permutation strategy is controlled by SuperLUOptions::ColPerm.

SuperLUMatrix::NumericalFactorize

Performs LU factorization numerically.

Example

#include "ppexsi.hpp"
{
  ...;
  // Construct AMat
  DistSparseMatrix<Complex>  AMat;
  ...;

  // Setup SuperLU
  SuperLUGrid g( comm, nprow, npcol );
  SuperLUOptions luOpt;
  luOpt.ColPerm = "MMD_AT_PLUS_A";
  SuperLUMatrix luMat( g );

  // Matrix conversion
  luMat.DistSparseMatrixToSuperMatrixNRloc( AMat );

  // Symbolic factorization
  luMat.SymbolicFactorize();

  // Numerical factorization
  luMat.NumericalFactorize();

  ...;
}

Reuse symbolic factorization

In SuperLU_DIST, the same symbolic factorization can be reused for factorizing different matrices. To reuse the symbolic factorization, one should follow the steps below.

(After symbolic factorization)

• Destroy the SuperMatrix.
• Convert another DistSparseMatrix into the native format (SuperMatrix in SuperLU_DIST) of the factorization routine.
• Redistribute the matrix into 2D block cyclic format.
• Numerical factorization.

Related structures and subroutines

SuperLUMatrix::DestroyAOnly

Releases the data in A but keeps other data, such as LUstruct.

This allows one to perform factorization of matrices of the same pattern, such as the option

fact = SamePattern_SameRowPerm

in SuperLU_DIST.

SuperLUMatrix::Distribute

Distribute redistrbutes the SuperMatrix in parallel so that it is ready for the numerical factorization.

Example

#include "ppexsi.hpp"
{
  ...;
  // Construct AMat
  DistSparseMatrix<Complex>  AMat;
  ...;

  // Setup SuperLU
  SuperLUGrid g( comm, nprow, npcol );
  SuperLUOptions luOpt;
  luOpt.ColPerm = "MMD_AT_PLUS_A";
  SuperLUMatrix luMat( g );

  // Matrix conversion
  luMat.DistSparseMatrixToSuperMatrixNRloc( AMat );

  // Symbolic factorization
  luMat.SymbolicFactorize();

  // Destroy the SuperMatrix saved in luMat.
  luMat.DestroyAOnly();


  // Construct another matrix BMat with the same sparsity pattern as A.
  DistSparseMatrix<Complex>  BMat; 
  ...;
  // Matrix conversion
  luMat.DistSparseMatrixToSuperMatrixNRloc( BMat );
  // Redistribute into 2D block cyclic format.
  luMat.Distribute();


  // Numerical factorization
  luMat.NumericalFactorize();

  ...;
}

Triangular solve and accuracy check

The triangular solve routines provided by SuperLU_DIST can be used to check the accuracy of the factorization as well as the selected inversion.

(After numericl factorization)

• Construct the distributed right hand sides.
• Solve \(Ax=b\). Multiple right hand sides can be solved simultaneously.

Related structures and subroutines

SuperLUMatrix::SolveDistMultiVector

Solve A x = b with b overwritten by x for distributed multivector.

SuperLUMatrix::CheckErrorDistMultiVector

Print out the error by direct comparison with the true solution in distributed format.

Example

The following example performs factorization, solves for a series of right hand sides and compare the accuracy.

#include "ppexsi.hpp"
{
  ...;
  // Construct AMat
  DistSparseMatrix<Complex>  AMat;
  ...;

  // Setup SuperLU
  SuperLUGrid g( comm, nprow, npcol );
  SuperLUOptions luOpt;
  luOpt.ColPerm = "MMD_AT_PLUS_A";
  SuperLUMatrix luMat( g );

  // Matrix conversion
  luMat.DistSparseMatrixToSuperMatrixNRloc( AMat );

  // Symbolic factorization
  luMat.SymbolicFactorize();

  // Numerical factorization
  luMat.NumericalFactorize();

  // Construct a global matrix (for error checking)
  SuperLUMatrix A1( g ), GA( g );
  A1.DistSparseMatrixToSuperMatrixNRloc( AMat );
  A1.ConvertNRlocToNC( GA );

  // Construct the distributed right hand sides and the exact solution.
  CpxNumMat xTrueGlobal(n, nrhs), bGlobal(n, nrhs);
  CpxNumMat xTrueLocal, bLocal;
  UniformRandom( xTrueGlobal );
  GA.MultiplyGlobalMultiVector( xTrueGlobal, bGlobal );
  A1.DistributeGlobalMultiVector( xTrueGlobal, xTrueLocal );
  A1.DistributeGlobalMultiVector( bGlobal,     bLocal );


  // Solve and check the error.
  luMat.SolveDistMultiVector( bLocal, berr );
  luMat.CheckErrorDistMultiVector( bLocal, xTrueLocal );

  ...;
}