Direct solvers for sparse matrices, including Gaussian
Elimination, frontal, multifrontal and supemodal, involve much
more complicated algorithms than for dense matrices. The
main complication is due to the need for efficiently handling
the fill-in in the factors L and U. A typical sparse solver
consists of four distinct steps as opposed to two in the dense
case:
1. An ordering step that reorders the rows and columns such
t that the factors suffer little fill, or that the matrix has
special structure, such as block-triangular form.
2. An analysis step or symbolic factorization that determines
the nonzero structures of the factors and creates suitable
data structures for the factors.
3. Numerical factorization that computes the L and U factors.
4. A solve step that performs forward and back substitution
using the factors.
4.1 Intel PARDISO 9.0
Intel® Math Kernel Library (Intel® MKL) offers highly
optimized, extensively threaded math routines for scientific,
engineering, and financial applications that require maximum
performance.
Multiprocessors with the PARDISO Direct Sparse Solver - an
easy-to-use, thread-safe, high-performance, and memory-
efficient software library licensed from the University of Basel.
Intel MKL also includes a Conjugate Gradient iterative solver
with a flexible reverse communication interface.
The package PARDISO is a thread-safe, high-performance,
robust, memory efficient and easy to use software for solving
large sparse symmetric and unsymmetric linear systems of
equations on shared memory multiprocessors. The solver uses a
combination of left- and right-looking Level-3 BLAS
supemode techniques to exploit pipelining parallelism and
memory hierarchies. In order to improve sequential and
parallel sparse numerical factorization performance, the
algorithms are based on a Level-3 BLAS update and pipelining
parallelism is exploited with a combination of left- and right
looking supemode techniques.
For sufficiently large problem sizes numerical experiments
demonstrate that the scalability of the parallel algorithm is
nearly independent of the shared-memory multiprocessing
architecture and speedup of up to seven using eight processors
have been observed. The approach is based on OpenMP
directives and has been successfully tested on the following
parallel computers: COMPAQ AlphaServer, AMD Opteron,
SGI Origin 2000, SUN Enterprise Server, Intel Pentium
IV/Itanium(R) with RedHat LINUX, NEC SX-5, and Cray J90
and IBM Regatta SMPs.
4.2 UMFPACK
UMFPACK is a set of routines for solving unsymmetric sparse
linear systems, Ax = b, using the Unsymmetric MultiFrontal
method and direct sparse LU factorization. It is written in
ANSI/ISO C, with a MATLAB interface. UMFPACK relies on
the Level-3 Basic Linear Algebra Subprograms (dense matrix
multiply) for its performance. This code works on Windows
and many versions of Unix (Sun Solaris, Red Hat Linux, IBM
AIX, SGI IRIX, and Compaq Alpha).
UMFPACK was developed by Timothy A. Davis from the
University of Florida, Gainesville, USA. Distributed under the
GNU LGPL license.
The famous math software MATLAB also use UMFPACK as
its solver. It represent today’s highest level and many other
solvers consider it as standard to compare.
4.3 GRUS SPARSE SOLVER
GSS(GRUS SPARSE SOLVER) is a software for large sparse
matrix solving. It is written in ANSI/ISO C,can be work well
with many sort of data and most of operating system platform.
GSS uses new models and many new algorithms at the
symbolic-analysis phase. GSS uses multifrontal method to
factor A to LU, it supports complete pivoting, partial
pivoting ,rook-pivoting and static pivoting . The tests by
comparing the time cost of computing large mounts of matrices
show that its average capability even better than UMFPACK in
many field such as the max matrix order number they could
work and their computing speed, the time GSS only cost one
third of that UMFPACK used.
/ >jyd xoiciou x/- an itijup: itfiflwXOK'GE3 /lda
There is no single algorithm or software that is best for all
types of linear systems. Some software is targeted for special
matrices such as symmetric and positive definite, some is
targeted for the most general cases.
5 APPLICATION IN PHOTOGRAMMETRY
Matrix computation are widely used in Photogrammetry such
as Complicated coordinate translation, Adjustment
calculation ,image enhance, image transform, image analysis
and other image processing techniques. For a image data can
naturally form an int matrix, image process can easily be
translated to matrix calculation.
Based on Meschach library the following example is core
content of a large block indirect adjustment programme in
photogrammetry filed to solve for the unknowns vector X of
adjustment values in the error equations:
V=AX-L weight P
The coefficient matrix A is a large band sparse matrix whose
order number may excess 10,000. L is the constant vector and
V is error vector. A,L,P are in parameters, the others are out
parameters.
/Ism Li to 3V)to3vjo:;b srfr , T 2o ¿aituAoo. adi qu coisrn k v k"io
#define MAXM 10
sp Adj ustLS( S PM AT *A, VEC *L, SPMAT *P, VEC
*Xout,VEC *V, double &m0)
{
SPMAT *AT,*Temp,*N;
VEC *b,*u,*ll;
int numsteps;
b-v_get(A->n);
u=v_get(A->m);
AT=sp_get( A->n, A->m,M AXM);
for(int i=0;i
m;i++)
for(int j=0;jn;j-H-)
sp set_val(AT,j ,i, sp_get_val(A,i,j));//AT=A A T
Temp=sp_get(A->n,A->m,MAXM);
N=sp_get( A->n, A->n,M AXM);
sp_mltadd(AT,P,1.0,Temp);
sp_mltadd(Temp,A,l .0,N);
sp mv mlt(Temp,L,b);
T emp=sp_resize(T emp, A->n, A->n);
Temp = sp_copy(N);
spICHfactor(T emp);
Xout = iter_spcg(N,Temp,b,le-7,Xout,1000,&num_steps);
//Xout out
11 = v_get(L->dim);
sv_mlt(-l,L,ll);
sp_mv_mlt(A,Xout,V);//V out
