PDF Index|
Contents
Functions
PDF Index |
Purpose
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 within 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.
To use UMFPACK as a ofact solver you need to install it on your machine. This code is available at www.cise.ufl.edu/research/sparse/umfpack.
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 is a sparse multi-frontal solver based on spooles a compiled version is provided with SDT distributions.
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 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.
Is used to define defaults for what the solver does.
This is the callback that is evaluated when ofact initializes a new matrix.
This is the callback that is evaluated when ofact overloads the matrix left division (\)
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 another computer on the network.
See also
