ofact
Purpose
Factored matrix object.
Syntax
ofact
ofact('method MethodName');
kd=ofact(k); q = kd\b; ofact('clear',kd);
kd=ofact(k,'MethodName')
Description
The factored matrix object ofact is designed to let users write code
that is independent of the library used to solve static problems of
the form [K]{q}={F}. For FEM applications, choosing the
appropriate library for that purpose is crucial. Depending on the case
you may want to use full, skyline, or sparse solvers. Then whithin
each library you may want to specify options (direct, iterative,
in-core, out-of-core, parallel, ... ).
Using the ofact object in your code, lets you specify method at run
time rather than when writing the code. Typical steps are
ofact('method spfmex'); % choose method
kd = ofact(k); % create object and factor
static = kd\b % solve
ofact('clear',kd) % clear factor when done
For single solves static=ofact(k,b) performs the three steps (factor, solve clear) in a single pass.
The first step of method selection provides an open architecture that lets users introduce new solvers with no need to rewrite functions that use ofact objects. Currently available methods are listed simply by typing
>> ofact
Available factorization methods for OFACT object
-> spfmex : SDT sparse LDLt solver
sp_util : SDT skyline solver
lu : MATLAB sparse LU solver
mtaucs : TAUCS sparse solver
pardiso : PARDISO sparse solver
chol : MATLAB sparse Cholesky solver
*psldlt : SGI sparse solver (NOT AVAILABLE ON THIS MACHINE)
and the method used can be selected with ofact('method MethodName'). SDtools maintains pointers to pre-compiled solvers at http://www.sdtools.com/faq/FE_ofact.html.
The factorization kd = ofact(k); and resolution steps static = kd\b can be separated to allow multiple solves with a
single factor. Multiple solves are essential for eigenvalue and
quasi-newton solvers. static = ofact(k)\b is of course also correct.
The clearing step is needed when the factors are not stored as
MATLAB variables. They can be stored in another memory pile, in an
out-of-core file, or on another computer/processor. Since for large problems, factors require a lot of memory. Clearing them is an important step.
Historically the object was called skyline. For backward
compatibility reasons, a skyline function is provided.
umfpack
To use UMFPACK as an ofact solver you need to install it on your
machine. This code is availlable at www.cise.ufl.edu/research/sparse/umfpack.
pardiso
Based on the Intel MKL (Math Kernel Library), you should use version 8 and after.
By default the pardiso call used in the ofact object is set for symmetric matrices. For non-symmetric matrices, you have to complement the ofact standard command for factorization with the character string 'nonsym'. Moreover, when you pass a matrix from Matlab to PARDISO, you must transpose it in order to respect the PARDISO sparse matrix format.
Assuming that k is a real non-symmetric matrix and b a real vector, the solution q of the system k.q=b is computed by the following sequence of commands:
ofact pardiso % select PARDISO solver
kd = ofact('fact nonsym',k'); % factorization
q=kd\b; % solve
ofact('clear',kd); % clear ofact object
The factorization is composed of two steps: symbolic and numerical factorization.
For the first step the solver needs only the sparse matrix structure (i.e. non-zeros location), whereas the actual data stored in the matrix are required in the second step only. Consequently, for a problem with a unique factorization, you can group the steps. This is done with the standard command ofact('fact',...).
In case of multiple factorizations with a set of matrices having the same sparse structure, only the second step should be executed for each factorization, the first one is called just for the first factorization. This is possible using the commands 'symbfact' and 'numfact' instead of 'fact' as follows:
kd = ofact('symbfact',k); % just one call at the beginning
...
kd = ofact('numfact',k,kd); % at each factorization
q=kd\b; %
...
ofact('clear',kd); % just one call at the end
You can extend this to non-symmetric systems as described above.
spfmex
spfmex is a sparse multi-frontal solver based on spooles a compiled version is provided with SDT distributions.
sp_util
The skyline matrix storage is a traditional form to store the
sparse symmetric matrices corresponding to FE models. For a full
symmetric matrix kfull
kfull=[1 2
10 5 8 14
6 0 1
9 7
sym. 11 19
20]
The non-zero elements of each column that are above the diagonal
are stored sequentially in the data field k.data from the diagonal up
(this is known as the reverse Jenning's representation) and the indices
of the elements of k corresponding to diagonal elements of the full
matrix are stored in an index field k.ind. Here
k.data = [1; 10; 2; 6; 5; 9; 0; 8; 11; 7; 1; 14; 20; 19; 0]
k.ind = [1; 2; 4; 6; 9; 13; 15];
For easier manipulations and as shown above, it is assumed in the that the index field k.ind has one more element
than the number of columns of
kfull whose value is the index of a zero which is added at the
end of the data field k.data.
If you have imported the ind and data fields from an external
code, ks = ofact (data, ind) will create the ofact object. You can
then go back to the MATLAB sparse format using sparse(ks) (but
this needs to be done before the matrix is factored when solving a static problem).
Your solver
To add your own solver simply add a file called MySolver_utils.m
in the @ofact directory. This function must accept the commands detailed below.
Your object can use the fields .ty used to monitor what is stored in the object (0 unfactored ofact, 1 factored ofact, 2 LU, 3
Cholesky, 5 other), .ind, .data used to store the matrix or factor in true ofact format, .dinv inverse of diagonal (currently unused), .l L factor in lu decomposition or transpose of Cholesky factor, .u U factor in lu decomposition or Cholesky factor, .method other free format information used by the object method.
method
Is used to define defaults for what the solver does.
fact
This is the callback that is evaluated when ofact initializes a new matrix.
solve
This is the callback that is evaluated when ofact overloads the matrix left division (\)
clear
clear is used to provide a clean up method when factor information is not stored within the ofact object itself. For example, in persistent memory, in another process or on an another computer on the network.
See also
fe_eig, fe_reduc
©1991-2007 by SDTools
