public class PardisoSolver extends java.lang.Object implements DirectSolver
Matrix M; // matrix to be solved VectorNd x, b; // solution vector and righthand side PardisoSolver solver = new PardisoSolver(); solver.analyze (M, M.rowSize(), Matrix.SYMMETRIC); // symbolic factorization solver.factor(); // numeric factorization solver.solve (x, b); // solution using factorization solver.dispose(); // release resources when we are doneIt is very important to call
dispose()
when the solver
is no longer required, in order to release internal native resources that
have been allocated for the Pardiso native code.
It is not necessary to call the analyze()
and
factor()
methods every time a solution is
required. analyze()
needs to be called only when a matrix is
first presented to the solver or when its sparsity structure changes.
factor()
needs to be called only when the numeric values of the
matrix change. For a given set of numeric values, once factor()
has been called, solve()
can be called as many times as desired
to generate solutions for different righthand sides.
It is also possible to avoid using a Matrix
object, and
instead call analyze()
and factor()
with
compressed row storage (CRS) data structures directly, as in
solver.analyze (vals, colIdxs, rowOffs, size, matrixType); solver.factor (vals);Here
vals
, colIdxs
, and rowOffs
descrive the sparse matrix structure using the CRS format as described
in the documentation for
Matrix.setCRSValues
.
Pardiso also supports iterative solving, using preconditioned CGS iteration.
The preconditioner is supplied by the most recent matrix factorization. If
the current matrix values are close to those associated with the
factorization, then the resulting iterative solution time can be
considerably faster than the alternative combination of a factor() and a
solve(). There are several iterativeSolve()
methods, which
obtain the current matrix values either directly as an input argument, or
from a Matrix
object supplied through the
analyze()
methods. When an iterativeSolve()
method
is successful, it returns a positive number indicating the number of
iterations performed. A possible call sequence is as follows:
solver.analyze (M, M.rowSize(), Matrix.SYMMETRIC); // symbolic factorization
solver.factor(); // numeric factorization
while (computing) {
... update matrix values and righthand side b ...;
if (solver.iterativeSolve(x, b) <= 0) {
// iterative solve failed; do explicit factor and solve
solver.factor();
solver.solve (x, b);
}
}
A more sophisticated version of the above code will also call
factor()
and solve()
when the compute time
associated with iterativeSolve()
exceeds a certain threshold.
After refactorization, the time required by iterativeSolve()
should be reduced since the next set of matrix values will again
(presumably) be close to those associated with the factorization. This
functionality is provided by the autoFactorAndSolve()
methods:
solver.analyze (M, M.rowSize(), Matrix.SYMMETRIC); // symbolic factorization solver.factor(); // numeric factorization while (computing) { ... update matrix values and righthand side b ...; // automatically choose between iterative and direct solving solver.autoFactorAndSolve (x, b); }
Modifier and Type  Class and Description 

static class 
PardisoSolver.ReorderMethod
Describes the reorder methods that can be used during the analyze phase
to reduce factorization fillin.

Modifier and Type  Field and Description 

static int 
ANALYZED
Indicates that a matrix has been set and analyzed for this solver.

static boolean 
DEFAULT_SHOW_PERTURBED_PIVOTS 
static int 
FACTORED
Indicates that a matrix has been set, analyzed, and numerically factored
for this solver.

static boolean 
printThreadInfo 
static boolean 
supportsMultipleRhs 
static int 
UNSET
Indicates that no matrix is currently set for this solver.

Constructor and Description 

PardisoSolver()
Creates a new PardisoSolver object.

Modifier and Type  Method and Description 

void 
analyze(double[] vals,
int[] colIdxs,
int[] rowOffs,
int size,
int type)
Sets the matrix associated with this solver and performs
symbolic analysis on it.

void 
analyze(Matrix M,
int size,
int type)
Sets the matrix associated with this solver and performs
symbolic analysis on it.

void 
analyzeAndFactor(Matrix M)
Convenience method that sets the matrix associated with this solver,
performs symbolic analysis on it, and factors it.

void 
autoFactorAndSolve(double[] x,
double[] b,
int tolExp)
Implementation of
autoFactorAndSolve(maspack.matrix.VectorNd,maspack.matrix.VectorNd,int)
that uses double[] objects
to store to the result and righthand side. 
void 
autoFactorAndSolve(VectorNd x,
VectorNd b,
int tolExp)
Calls
factor() and solve(double[],double[]) together,
or, if tolExp is positive, automatically determines
when to call
iterativeSolve(maspack.matrix.VectorNd,maspack.matrix.VectorNd,int)
with the specific tolExp instead,
depending on whether the matrix is factored and if it is estimated
that iterativeSolve will save time. 
void 
dispose()
Releases the native resources used by this solver.

void 
factor()
Performs a numeric factorization of the matrix associated with this
solver, using the current numeric values contained within
the matrix that was supplied by a previous call to
analyze(Matrix,int,int) or
analyzeAndFactor(Matrix) . 
void 
factor(double[] vals)
Performs a numeric factorization of the most recently analyzed matrix
solver using the supplied numeric values.

void 
finalize() 
int 
getAnalysisMemoryUsage()
Returns the amount of permanent memory (in kbytes) that was allocated
during the most recent
analyze() call. 
boolean 
getApplyScaling()
Returns whether or not matrix scaling is enabled.

boolean 
getApplyWeightedMatchings()
Returns whether or not weighted matchings are enabled.

static int 
getDefaultNumThreads()
Returns the default number of threads that Pardiso is assigned when a
PardisoSolver is created. 
java.lang.String 
getErrorMessage()
Returns a message describing the reason for failure of
the most recent call to any of the
analyze() ,
factor() , or iterativeSolve() methods,
or null if the method succeeded. 
int 
getFactorSolveMemoryUsage()
Returns the total memory consumption (in kbytes) required for the
factor() and solve() calls. 
static java.lang.String 
getInitErrorMessage()
Returns a message describing an error that occurred during
initialization, or
null if no error occurred. 
boolean 
getMatrixChecking()
Returns true if matrix checking is enabled.

int 
getMaxRefinementSteps()
Returns the maximum number of iterative refinement steps that Pardiso
should perform after a solve.

int 
getMessageLevel()
Returns the message level for the Pardiso native code.

int 
getNumNegEigenvalues()
Returns the number of negative eigenvalues that were detected during the
most recent numeric factorization of a symmetric indefinite matrix (i.e.,
during the last
factor() call). 
int 
getNumNonZerosInFactors()
Returns the number of nonzeros elements in the factorization.

int 
getNumPerturbedPivots()
Returns the number of pivot perturbations that were required
during the last recent numeric factorization (i.e., during
the last
factor() call). 
int 
getNumPosEigenvalues()
Returns the number of positive eigenvalues that were detected during the
most recent numeric factorization of a symmetric indefinite matrix (i.e.,
during the last
factor() call). 
int 
getNumRefinementSteps()
Returns the number of iterative refinement steps that Pardiso
actually performed during the most call to
solve() . 
int 
getNumThreads()
Returns the number of threads that Pardiso should use.

int 
getPeakAnalysisMemoryUsage()
Returns the peak amount of memory (in kbytes) that was used
during the most recent
analyze() call. 
int 
getPivotPerturbation()
Returns the size of the perturbation that should be used to resolve
zeropivots, expressed as a negative poweroften exponent.

PardisoSolver.ReorderMethod 
getReorderMethod()
Gets the reorder method that is used during the analyze phase to
reduced factorization fillin.

static boolean 
getShowPerturbedPivots()
Queries whether the "num perturbed pivots" message is enabled.

int 
getSPDZeroPivot()
If the solver detects a zero or negative pivot for an SPD matrix,
this method returns the row number where the first zero or negative
pivot was detected.

int 
getState()
Returns the current stateface for this solver.

boolean 
getUse2x2Pivoting()
Returns whether or not 2 x 2 pivoting is enabled.

boolean 
hasAutoIterativeSolving()
Returns true if this solver supports automatic iterative solving using a
recent directlyfactored matrix as a preconditioner.

static boolean 
isAvailable()
Returns true if the Pardiso solver is available.

int 
iterativeSolve(double[] vals,
double[] x,
double[] b,
int tolExp)
Attempts to use preconditioned CGS iteration to solve M x = b for a given
righthand side
b . 
int 
iterativeSolve(double[] x,
double[] b,
int tolExp)
Attempts to use preconditioned CGS iteration to solve M x = b for a given
righthand side
b . 
int 
iterativeSolve(VectorNd x,
VectorNd b,
int tolExp)
Implementation of
iterativeSolve(double[],double[],int)
that uses VectorNd objects
to store to the result and righthand side. 
static void 
main(java.lang.String[] args) 
double 
residual(int[] rowOffs,
int[] colIdxs,
double[] vals,
int size,
double[] x,
double[] b,
boolean symmetric)
Computes the norm of the residual

void 
setApplyScaling(int enable)
Sets whether or not Pardiso should apply matrix scaling to its
factorizations.

void 
setApplyWeightedMatchings(int enable)
Sets whether or not Pardiso should apply weighted matching to its
factorizations.

static void 
setDefaultNumThreads(int num)
Sets the default number of threads that Pardiso is assigned when a
PardisoSolver is created. 
void 
setMatrixChecking(boolean enable)
Enables or disables checking of the integrity of the matrix data
structures passed to Pardiso.

void 
setMaxRefinementSteps(int nsteps)
Sets the maximum number of iterative refinement steps that Pardiso should
perform after a solve.

void 
setMessageLevel(int level)
Sets the message level for the Pardiso native code.

void 
setNumThreads(int num)
Sets the number of threads that Pardiso should use.

void 
setPivotPerturbation(int n)
Sets the size of the perturbation that should be used to resolve
zeropivots, expressed as a negative poweroften exponent.

void 
setReorderMethod(PardisoSolver.ReorderMethod method)
Sets the reorder method that should be used during the analyze phase to
reduced factorization fillin.

static void 
setShowPerturbedPivots(boolean enable)
Enables/disables the "num perturbed pivots" message (which usually
indicates an illconditioned solve).

void 
setUse2x2Pivoting(int enable)
Sets whether or not Pardiso should use 2 x 2 Bunch and Kaufman
pivoting (in addition to regular 1 x 1 pivoting) for symmetric
indefinite matrices.

void 
solve(double[] x,
double[] b)
Solves the matrix associated with this solver for x, given a
specific righthand side b.

void 
solve(double[] X,
double[] B,
int nrhs)
Solves the matrix associated with this solver for a set of vectors X,
given a set of righthand sides B.

void 
solve(VectorNd x,
VectorNd b)
Solves the matrix associated with this solver for x, given a
specific righthand side b.

void 
systemExit(int code)
This method is a hook that gives us access to _exit(), which is
is in turn needed to exit ArtiSynth on the MacBook Pro 8.2
version of Ubuntu, since otherwise the JVM exit process encounters
a SEGV in XQueryExtension.

public static boolean printThreadInfo
public static boolean supportsMultipleRhs
public static boolean DEFAULT_SHOW_PERTURBED_PIVOTS
public static final int UNSET
public static final int ANALYZED
public static final int FACTORED
public int getState()
public static boolean isAvailable()
public static java.lang.String getInitErrorMessage()
null
if no error occurred.public void analyze(Matrix M, int size, int type)
ANALYZED
.
Normally the matrix is simply supplied by the argument M
,
unless the size
arugment is less than M.rowSize()
,
in which case the matrix is taken to be the topleft diagonal submatrix
of the indicated size. This solver retains a pointer to M
until the next call to analyze()
or
analyzeAndFactor()
.
The type of the matrix is given by type
:
analyze
in interface DirectSolver
M
 supples the matrix to be analyzedsize
 size of the matrix to be analyzedtype
 type of the matrix to be analyzedjava.lang.IllegalArgumentException
 if the matrix is not square or if
size
is out of boundsNumericalException
 if the matrix cannot be analyzed for numeric
reasons.public void factor()
analyze(Matrix,int,int)
or
analyzeAndFactor(Matrix)
.
After calling this method, the solver's state is set to
FACTORED
.factor
in interface DirectSolver
ImproperStateException
 if not preceded by a call to
analyze(Matrix,int,int)
or
analyzeAndFactor(Matrix)
NumericalException
 if the matrix cannot be factored for numeric
reasons.public void factor(double[] vals)
FACTORED
.vals
 nonzero matrix element valuesImproperStateException
 if this solver's state is
UNSET
or if the number of supplied values is
less that the number of nonzero elements in the analyzed matrix.NumericalException
 if the matrix cannot be factored for numeric
reasons.public void analyzeAndFactor(Matrix M)
Matrix.INDEFINITE
, meaning
that Pardiso will produce an L U factorization. After calling this
method, the solver's state is set to ANALYZED
.analyzeAndFactor
in interface DirectSolver
M
 matrix to factorNumericalException
 if the matrix cannot be analyzed or factored
for numeric reasons.public void autoFactorAndSolve(VectorNd x, VectorNd b, int tolExp)
factor()
and solve(double[],double[])
together,
or, if tolExp
is positive, automatically determines
when to call
iterativeSolve(maspack.matrix.VectorNd,maspack.matrix.VectorNd,int)
with the specific tolExp
instead,
depending on whether the matrix is factored and if it is estimated
that iterativeSolve
will save time.
It is assumed that a matrix was supplied to the solver using a previous
call to
analyze(Matrix,int,int)
or
analyzeAndFactor(Matrix)
.autoFactorAndSolve
in interface DirectSolver
x
 returns the solution valueb
 supplies the righthand sidetolExp
 if positive, enables iterative solving and provides
the exponent of the stopping criterionjava.lang.IllegalArgumentException
 if the dimensions of x
or
b
are incompatible with the matrix size, or if
topExp
is negative.ImproperStateException
 if not preceded by a call to
analyze(Matrix,int,int)
or
analyzeAndFactor(Matrix)
NumericalException
 if the matrix cannot be factored for numeric
reasonspublic void autoFactorAndSolve(double[] x, double[] b, int tolExp)
autoFactorAndSolve(maspack.matrix.VectorNd,maspack.matrix.VectorNd,int)
that uses double[]
objects
to store to the result and righthand side.x
 returns the solution valueb
 supplies the righthand sidetolExp
 if positive, enables iterative solving and provides
the exponent of the stopping criterionjava.lang.IllegalArgumentException
 if the dimensions of x
or
b
are incompatible with the matrix size, or if
topExp
is negative.ImproperStateException
 if not preceded by a call to
analyze(Matrix,int,int)
or
analyzeAndFactor(Matrix)
NumericalException
 if the matrix cannot be factored for numeric
reasonspublic static void setShowPerturbedPivots(boolean enable)
enable
 if true
, enables the messagepublic static boolean getShowPerturbedPivots()
true
if the message is enabledpublic static void setDefaultNumThreads(int num)
PardisoSolver
is created. The results are undefined if this
number exceeds the maximum number of threads available on the
system. Setting num
to a value <=
0 will reset the
number of threads to the default used by OpenMP, which is typically the
value stored in the environment variable OMP_NUM_THREADS
.num
 default number of threads to usegetDefaultNumThreads()
public static int getDefaultNumThreads()
PardisoSolver
is created.setDefaultNumThreads(int)
public void setNumThreads(int num)
num
to a value <=
0 will
reset the number of threads to the default used by OpenMP, which is
typically the value stored in the environment variable
OMP_NUM_THREADS
.
Note: under the MKL version of Pardiso, changes to the thread number are applied globally to all instances of Pardiso running in the same process. Therefore, this method is of limited utility and does not allow different numbers of threads to be used by different PardisoSolver instances. Moreover, the thread number should not be changed in between the analyze, factor and solve phases.
num
 number of threads to usegetNumThreads()
public int getNumThreads()
OMP_NUM_THREADS
.setNumThreads(int)
public void setMaxRefinementSteps(int nsteps)
nsteps
 maximum number of iterative refinement stepsgetMaxRefinementSteps()
,
getNumRefinementSteps()
public int getMaxRefinementSteps()
setMaxRefinementSteps(int)
,
getNumRefinementSteps()
public int getNumRefinementSteps()
solve()
.getMaxRefinementSteps()
,
setMaxRefinementSteps(int)
public PardisoSolver.ReorderMethod getReorderMethod()
setReorderMethod(maspack.solvers.PardisoSolver.ReorderMethod)
public void setReorderMethod(PardisoSolver.ReorderMethod method)
DEFAULT
will cause Pardiso to choose its default value.method
 reorder methodgetReorderMethod()
public void setPivotPerturbation(int n)
norm(M) 10^{n}where
n
is the argument to this method and
norm(M)
is a norm of the matrix. Setting n
this to 1 will cause Pardiso to choose a default value appropriate to
the matrix type.n
 perturbation exponentgetPivotPerturbation()
public int getPivotPerturbation()
setPivotPerturbation(int)
public void setApplyScaling(int enable)
setApplyWeightedMatchings()
) are also
selected. Scaling is controlled by enable
as follows:
enable > 0
enables scaling, enable = 0
disables
scaling, and enable < 0
causes the solver to choose a
default value appropriate to the matrix type.enable
 enables/disables matrix scalinggetApplyScaling()
public boolean getApplyScaling()
setApplyScaling(int)
public void setApplyWeightedMatchings(int enable)
setApplyScaling()
) is recommended for highly indefinite symmetric
systems (such as saddle point problems) and is the default for symmetric
matrices. Weighted matchings are controlled by enable
as
follows: enable > 0
enables them, enable = 0
disables them, and enable < 0
causes the solver to choose a
default value appropriate to the matrix type.enable
 enables/disables weight matchingsgetApplyWeightedMatchings()
public boolean getApplyWeightedMatchings()
setApplyWeightedMatchings(int)
public void setUse2x2Pivoting(int enable)
enable
as
follows: enable > 0
enables it, enable = 0
disables it, and enable < 0
causes the solver to choose a
default value.enable
 enables 2 x 2 pivotinggetUse2x2Pivoting()
public boolean getUse2x2Pivoting()
setUse2x2Pivoting(int)
public void setMatrixChecking(boolean enable)
enable
 enables/disables matrix checkinggetMatrixChecking()
public boolean getMatrixChecking()
setMatrixChecking(boolean)
public void setMessageLevel(int level)
level
 native code message levelgetMessageLevel()
public int getMessageLevel()
setMessageLevel(int)
public int getNumNonZerosInFactors()
analyze()
call.public int getNumNegEigenvalues()
factor()
call). For matrices that are
not symmetric indefinite, 1 is returned.public int getNumPosEigenvalues()
factor()
call). For matrices that are
not symmetric indefinite, 1 is returned.public int getNumPerturbedPivots()
factor()
call). Pivot perturbation
generally indicates a singular, or very nearly singular, matrix.public int getSPDZeroPivot()
public int getPeakAnalysisMemoryUsage()
analyze()
call.public int getAnalysisMemoryUsage()
analyze()
call.public int getFactorSolveMemoryUsage()
factor()
and solve()
calls. This number is
determined in the most recent factor()
call.public void analyze(double[] vals, int[] colIdxs, int[] rowOffs, int size, int type)
ANALYZED
.
The matrix structure and its initial values are described using a
compressed row storage (CRS) format. See Matrix.setCRSValues()
for a detailed
description of this format. The matrix type is the same as that supplied
to analyze(Matrix,int,int)
. It is not possible to call factor()
after calling this version of analyze
because it
does not supply a matrix that can be used to obtain values from.
vals
 values of the nonzero matrix elements. These may be used to
assist the symbolic factorization, but will not used in any actual
numeric factorization.colIdxs
 1based column indices of the nonzero matrix elements.rowOffs
 1based row start offsets into vals
and
colIdxs
, corresponding to CRS format.size
 size of the matrix to be analyzedtype
 type of the matrix to be analyzedjava.lang.IllegalArgumentException
 if the CRS data structures
are inconsistent.NumericalException
 if the matrix cannot be analyzed for numeric
reasons.public void solve(VectorNd x, VectorNd b)
FACTORED
.solve
in interface DirectSolver
x
 returns the solution valueb
 supplies the righthand sidejava.lang.IllegalStateException
 if this solver's state is not
FACTORED
java.lang.IllegalArgumentException
 if the dimensions of x
or
b
are incompatible with the matrix size.public void solve(double[] x, double[] b)
FACTORED
.x
 returns the solution valueb
 supplies the righthand sidejava.lang.IllegalStateException
 if this solver's state is not
FACTORED
java.lang.IllegalArgumentException
 if the dimensions of x
or
b
are incompatible with the matrix size.public void solve(double[] X, double[] B, int nrhs)
nrhs
. Both X
and B
should be stored in
the column major ordering used by FORTRAN. It is assumed that the matrix
has been factored and that this solver's state is
FACTORED
.X
 returns the solutions in column major orderB
 supplies the righthand sides in column major ordernrhs
 number of righthand sides to solve forjava.lang.IllegalStateException
 if this solver's state is not
FACTORED
java.lang.IllegalArgumentException
 if the dimensions of X
or
B
are incompatible with the matrix size and rhs
.public double residual(int[] rowOffs, int[] colIdxs, double[] vals, int size, double[] x, double[] b, boolean symmetric)
M x  bfor a given values of M, x, and b. The values of
M
are given in compressed row storage (CRS) format.rowOffs
 matrix row offsets (CRS format)colIdxs
 nonzero element column indices (CRS format)vals
 nonzero element value (CRS format)size
 size of the matrixx
 supplies the solution valueb
 supplies the righthand sidesymmetric
 if true
, assumes that the arguments
define only the upper triangular portion of a symmetric matrix.java.lang.IllegalArgumentException
 if the dimensions of x
or
b
are incompatible with the matrix size.public int iterativeSolve(VectorNd x, VectorNd b, int tolExp)
iterativeSolve(double[],double[],int)
that uses VectorNd
objects
to store to the result and righthand side.x
 returns the solution valueb
 supplies the righthand sidetolExp
 exponent for the stopping criterionImproperStateException
 if not preceded by a call to
analyze(Matrix,int,int)
or
analyzeAndFactor(Matrix)
,
or if the matrix has not previously been factored.java.lang.IllegalArgumentException
 if the dimensions of x
or
b
are incompatible with the matrix size, or if
topExp
is negative.public int iterativeSolve(double[] x, double[] b, int tolExp)
b
. Current numeric values for M are obtained
from the matrix that was supplied by a previous call to analyze(Matrix,int,int)
or analyzeAndFactor(Matrix)
. The
most recent numeric factorization of this matrix will be used as a
preconditioner for the CGS iteration, with a relative stopping criterion
given by 10^tolExp
. It is assumed that the matrix has
been previously factored with a call to factor()
.
If the CGS iteration is successful, this method returns a positive
value giving the number of iterations required. If unsuccessful,
it returns a nonpostive value giving the negative of the number
of iterations that were actually performed, and
getErrorMessage()
can be used to determine
the underlying error.
x
 returns the solution valueb
 supplies the righthand sidetolExp
 exponent for the stopping criterionImproperStateException
 if not preceded by a call to
analyze(Matrix,int,int)
or
analyzeAndFactor(Matrix)
,
or if the matrix has not previously been factored.java.lang.IllegalArgumentException
 if the dimensions of x
or
b
are incompatible with the matrix size, or if
topExp
is negative.public int iterativeSolve(double[] vals, double[] x, double[] b, int tolExp)
b
. Current numeric values for M are
supplied by the argument vals
. Otherwise this
method behaves identically to
iterativeSolve(double[],double[],int)
.vals
 supplied the current matrix valuesx
 returns the solution valueb
 supplies the righthand sidetolExp
 exponent for the stopping criterionImproperStateException
 if the matrix has not previously been
factored.java.lang.IllegalArgumentException
 if there are insufficient values
specified by vals
, the dimensions of x
or
b
are incompatible with the matrix size, or if
topExp
is negative.public void dispose()
dispose
in interface DirectSolver
public void finalize()
finalize
in class java.lang.Object
public boolean hasAutoIterativeSolving()
autoFactorAndSolve
method.hasAutoIterativeSolving
in interface DirectSolver
public static void main(java.lang.String[] args)
public java.lang.String getErrorMessage()
analyze()
,
factor()
, or iterativeSolve()
methods,
or null
if the method succeeded.public void systemExit(int code)