After compiling the mex interface drivers from the ILUPACK directory, add the associated system path, e.g., using the MATLAB command
>> addpath '...../ilupack'
After that you will able to use the ILUPACK toolbox for MATLAB.Suppose you would like to solve a system Ax=b, using the ILUPACK toolbox in the simplest case you would do this as follows
>> [PREC,options]=AMGfactor(A);
Then the preconditioner is built using the default options.Since ILUPACK has only built a preconditioner, to solve Ax=b we still need to invoke an iterative solver. This is done by
>> [x, options] = AMGsolver(A, PREC, options, b);
If you do not need the preconditioner anymore discard it using AMGdelete.
>> PREC = AMGdelete(PREC);
This will release any (also internal) memory that has been allocated for using PREC.
>> PREC = AMGdelete(PREC);
and rebuild the preconditioner with the revised options>> [PREC,options]=AMGfactor(A,options);
and try again AMGsolver, e.g.
>> PREC = AMGdelete(PREC);
>> options.droptol=1e-3;
>> [PREC,options]=AMGfactor(A,options);
>> [x,options]=AMGsolver(A,PREC,options,b);
>> options.maxit=10000;
>> [x,options]=AMGsolver(A, PREC, options, b,x);
>> options.restol=1e-14;
>> [x,options]=AMGsolver(A, PREC, options,b,x);
The main philosophy of ILUPACK can be summarized in three major components.
--> |
'amd' Approximate Minimum Degree (see MATLAB's symamd function)
'metisn', 'metise' Metis multilevel nested dissection by nodes, resp. edges. This is supported in MATLAB using symmetisn / symmetise interface drivers.
'rcm' Reverse Cuthill McKee (see MATLAB's symrcm function)
'mmd' Multiple Minimum Degree (MATLAB's outdated symmmd function)
ILUPACK offers many, many parameters that you may alter.
Some parameters you will find familiar, others may confuse you ...
Default Parameters
>> options=AMGinit(A)
options =
matching: 1
ordering: 'amd'
droptol: 0.0100
droptolS: 0.0010
droptolc: 1e-12
condest: 5
restol: 3.2927e-10
maxit: 500
elbow: 10
lfil: 156116
lfilS: 156116
typetv: 'none'
tv: [156115x1 double]
amg: 'ilu'
npresmoothing: 1
npostsmoothing: 1
ncoarse: 1
presmoother: 'gsf'
postsmoother: 'gsb'
FCpart: 'none'
typecoarse: 'ilu'
solver: 'gmres'
damping: 0.6667
nrestart: 30
ind: [156115x1 double]
isdefinite: 0
mixedprecision: 0
contraction: 0.5
coarsereduce: 1
decoupleconstraints: 0
matching | improve diagonal dominance using maximum weighted matchings |
ordering | preprocess the system by a symbolic reordering (e.g. 'amd', Approximate Minimum Degree) |
droptol | threshold to drop small entries during the factorization |
droptolS | threshold to drop small entries from the Schur complement (coarse grid system) |
droptolc | threshold for dropping small entries from the constraint part, if present |
condest | bound for the inverse triangular factors |
restol | stopping criterion for the iterative process (backwar error) |
maxit | maximum number of iteration steps |
elbow | memory space to keep the preconditioner (relative to the number of nonzeros of A) |
lfil | number of nonzeros per row in the approximate factorization |
lfilS | number of nonzeros per row in the approximate Schur complement |
typetv | kind of test vector (default 'none', alternatively 'static' for providing a vector for which the multilevel ILU is exact) |
tv | test vector such that the factorization is exact for this vector |
amg | kind of multilevel method (default: 'ilu' no inner iteration on the Schur complement system, alternatively 'amli' for inner iteration or 'mg' for multigrid-like set up |
npresmoothing | number of pre smoothing steps (if param.amg='mg' is chosen) |
npostsmoothing | number of post smoothing steps (if param.amg='mg' is chosen) |
ncoarse | number of coarse grid correction steps (only if param.amg='amli' or 'mg' is chosen) |
presmoother | type of pre smoother (e.g. 'gsf' Gauss-Seidel forward, only if param.amg='mg' is chosen) |
postsmoother | type of post smoother (e.g. 'gsb' Gauss-Seidel backward, only if param.amg='mg' is chosen) |
FCpart | decide whether nodes should be grouped a priori into some 'fine' grid nodes or 'coarse' grid nodes (default 'none') |
typecoarse | type of coarse grid system (define how the approximation is computed, e.g. 'ilu' is the easiest and fastest but least accurate coarse grid system) |
solver | which iterative solver is used (e.g. 'gmres') |
damping | damping factor for Jacobi smoothing |
nrestart | number of steps before GMRES is restarted |
ind | indicator array to indicate saddle point structure |
isdefinite | tell ILUPACK in advance that you matrix is positive definite |
mixedprecision | compute preconditioner in single precision |
contraction | contraction factor of the residual for inner flexible solver when AMLI or classical multigrid is used |
coarsereduce | If different from zero, then the L21 and the U12 block are discarded solving with L,U is done implicitly via L11,U11 and A21 (resp. A12). |
decoupleconstraints | This allows for saddle point type problems to explictly decouple the connections between the constraint part and the free part. |
Display Multilevel ILU
The multilevel ILU that has been computed by ILUPACK is passed to PREC. In PREC you can find the reorderings for any level, the partial incomplete factorizations as well as the matrices that interact between two neighbouring levels.
If you would like display the multilevel ILU you simply use AMGspy, e.g.
>> [PREC,options]=AMGfactor(A,options);
>> AMGspy(PREC);
The L factor is displayed using green colors, for U blue
colors are used and finally the transfer blocks have read colors.
If you also would like to see how the original system behaves associated
with the reordering and and the multilevel history, you may use
>> AMGspy(A,PREC);
Number of nonzeros consumed by the multilevel ILU
>> AMGnnz(PREC)
Solve a single preconditioned system LUx=b. This might be useful if you plan to create you own iterative solver.
>> x=AMGsol(PREC,b);
Load / save matrix in Harwell-Boeing format.
>> [A,rhs,rhstyp]=loadhbo(filename);
>> savehbo(filename, A);
This routine is also available for right side, initial guess, ...
However, ILUPACK asks YOU to provide whether the system is positive definite or not. You can do that by setting options.isdefinite before you factor your matrix, e.g.
>> options.isdefinite=1;
>> options=AMGinit(A,options);
>> [PREC,options]=AMGfactor(A,options);
ILUPACK offers YOU a tool to convert a symmetric preconditioner into a positive definite one using AMGconvert. This has to be done, after the preconditioner is built, e.g.
>> [PREC,options]=AMGfactor(A,options);
>> PREC = AMGconvert(PREC);
options = AMGinit(A) options = AMGinit(A,options) init structure of options to their default values for a given nxn matrix A input ----- A nxn matrix options optional input to indicate via options.isdefinite that your matrix is positive definite output ------ options structure with default parameters ------------------------------------------------------------------------ possible options: 1. options.matching -------------------- if different from zero then maximum weight matching will be used, for the real symmetric, complex Hermitian or complex symmetric the associated symmetric maximum weight matching will be used. `maximum weight matching' is a technique to reorder and rescale the matrix such that it is becomes more diagonally dominant default value: 1 2. options.ordering -------------------- several reorderings based on |A|+|A|' are offered. Unsymmetric patterns treated as if they were symmetric. The orderings are repeated on any coarser level. `amd' default, Approximate Minimum Degree `metisn' METIS multilevel nested dissection by NODES `metise' METIS multilevel nested dissection by EDGES `rcm' Reverse Cuthill-McKee `mmd' Minimum Degree `amf' Approximate Minimum Fill any other no reordering 3. options.droptol ------------------- threshold for dropping small entries during the computation of the incomplete LU decomposition default: 1e-2 4. options.droptolS -------------------- threshold for dropping small entries from the Schur complement default: 1e-2 5. options.droptolc -------------------- threshold for dropping small entries from the constraint part, if present (as indicated by negative entries in options.ind) default: 0 6. options.condest ------------------- bound for the inverse triangular factors from the incomplete LU decomposition default: 100 7. options.restol ------------------ bound for the accuracy of the approximate solution of Ax=b for a given right hand side after using ILUPACK-preconditioned iterative solution This tolerance refers to the BACKWARD ERROR, in the symmetric(Hermitian) positive definite case to the relative error in the energy norm default: sqrt(eps) 8. options.maxit ----------------- maximum number of iteration steps, before the iterative ILUPACK- preconditioned solvers terminates default: 500 10. options.elbow ----------------- elbow space for memory of the ILUPACK multilevel preconditioner. Since the core part of the code is FORTRAN 77, no dynamic memory allocation is available, one has to estimate the memory requirement in advance by multiples of the fill of the initial matrix default: 10 (ten times as much fill as the initial matrix) 11. options.lfil --------------- restrict the number of nonzeros per column in L (resp. per row in U) hard to at most `lfil' entries. default: n+1 12. options.lfilS ---------------- restrict the number of nonzeros per row in the approximate Schur complement hard to at most `lfilS' entries. default: n+1 13. options.typetv ------------------ define, whether to include a test vector into the computations If used, then (a) the test vector is also included to estimate the norm of the inverse triangular factors (b) it is used to obtain a refined fine/coarse grid partitioning (if switch is set) (c) diagonal compensation and off-diagonal lumping are added to improve the factorization default: 'none' alternatives: 'static' for a fixed test vector 'function_name' for a dynamically generated test vector the user has to provide a custom routine to generate the test vector format: tvd=function_name(mat,tv_old) On every level the associate matrix 'mat' is passed to this routine, also the coarse grid projection of the previous test vector 'tv_old' from the finer grid is passed. On exit, 'tvd' is the dynamic test vector that should be used 13. options.tv -------------- static test vector. Ignored if options.typetv=='none' If options.typetv=='function_name', then on the initial finest level, options.tv is passed to 'function_name' as initial guess for 'tv_old' 14. options.amg --------------- type of AMG preconditioner default: 'ilu' multilevel ILU alternatives: 'amli' multilevel ILU, where on each coarser grid an inner iteration is used to solve the coarse grid system. The number of inner interation steps is prescribed in options.ncoarse 'mg' classical multigrid with pre and post smoothing 15. options.npresmoothing ------------------------- number of pre smoothing steps (only needed if options.amg=='mg') default: 1 16. options.npostsmoothing -------------------------- number of post smoothing steps (only needed if options.amg=='mg') default: 1 17. options.ncoarse ------------------- number of coarse grid correction steps (only needed if options.amg=='mg' or options.amg=='amli') default: 1 if a negative value is used, then a flexible solver is used. 18. options.presmoother ------------------------- type of pre smoother default: 'gsf' Gauss-Seidel forward 'gsb' Gauss-Seidel backward 'j' Jacobi 'ilu' partial incomplete ILU on the F-nodes 'function_name' custom routine for smoothing d=function_name(mat,r) Given the matrix 'mat' on each level and the associated residual 'r(=b-mat*x_old)', the custom smoother is asked to compute an approximate defect 'd' such that x_new=x_old+d 19. options.postsmoother ------------------------- type of post smoother default: 'gsb' Gauss-Seidel backward 'gsf' Gauss-Seidel forward 'j' Jacobi 'ilu' partial incomplete ILU on the F-nodes 'function_name' custom routine for smoothing d=function_name(mat,r) Given the matrix 'mat' on each level and the associated residual 'r(=b-mat*x_old)', the custom smoother is asked to compute an approximate defect 'd' such that x_new=x_old+d 20. options.FCpart ------------------ preselect a partitioning into fine grid and coarse grid nodes default: 'none' 'yes' if F/C partioning is desired 21. options.typecoarse ---------------------- type of coarse grid system default: 'ilu' coarse grid system is computed from the approximate incomplete factorization by ignoring that entries have been discarded 'amg' Use the associated approximate interpolation operator P and the restriction operator R from the underlying inverse triangular factors to generate the Galerkin type coarse grid matrix R A P 22. options.damping ---------------------- damping factor if Jacobi smoothing is chosen 23. options.isdefinite ---------------------- if given on input then the matrix is assumed to be symmetric (Hermitian) positive definite and the parameters will be 24. options.ind --------------- indicate by negative signs which parts of the system refer to the second (typically) zero block of a saddle point system 25. options.mixedprecision --------------- if different from zero, then single precision for preconditioning is used 26. options.contraction ----------------------- contraction of the residual for inner flexible solver when AMLI or classical multigrid is used and options.ncoarse<0 (flexible coarse grid solver) 27. options.coarsereduce ------------------------ default: 1. If different from zero, then the L21 and the U12 block are discarded solving with L,U is done implicitly via L11,U11 and A21 (resp. A12). If set to zero, then L21, U12 are kept 28. options.decoupleconstraints ------------------------------- default: 1. This allows for saddle point type problems to explictly decouple the connections between the constraint part and the free part. Applied on every level, this allows for smaller coarse grid matrices. If set to zero, then the additional decoupling is not applied
[PREC, options] = AMGfactor(A, options); [PREC, options] = AMGfactor(A); Computes ILUPACK preconditioner PREC according to the given options. For details concerning `options' see `AMGinit' input ----- A nxn nonsingular matrix options parameters. If `options' is not passed then the default options from `AMGinit' will be used output ------ PREC ILUPACK multilevel preconditioner PREC is a structure of `nlev=length(PREC)' elements indicating the number of levels. For every level l we have PREC(l).n size of level l PREC(l).nB size of the leading block of level l PREC(l).L (block) lower triangular matrix PREC(l).D (block) diagonal matrix PREC(l).U if present: (block) upper triangular matrix for symmetrically structured matrices only L is present. L*D^{-1}*U is the approximate LU decomposition of the leading block of A after rescaling and reordering PREC(l).E except for l
[x, options] = AMGsolver(A, PREC, options, b, x0) [x, options] = AMGsolver(A, PREC, options, b) [x, options] = AMGsolver(A, PREC, b, x0) [x, options] = AMGsolver(A, PREC, b) Solves Ax=b using ILUPACK preconditioner PREC according to the given options [x, options] = AMGhsolver(A, PREC, options, b, x0) [x, options] = AMGhsolver(A, PREC, options, b) [x, options] = AMGhsolver(A, PREC, b, x0) [x, options] = AMGhsolver(A, PREC, b) Solves A'x=b, where A' is the (conjugate) transposed matrix w.r.t. A. This is done using ILUPACK preconditioner PREC according to the given options [x, options] = AMGtsolver(A, PREC, options, b, x0) [x, options] = AMGtsolver(A, PREC, options, b) [x, options] = AMGtsolver(A, PREC, b, x0) [x, options] = AMGtsolver(A, PREC, b) Solves A.'x=b, where A.' is the transposed matrix w.r.t. A. This is done using ILUPACK preconditioner PREC according to the given options
PREC=AMGdelete(PREC) delete preconditioner, in particular release the associated memory
x = AMGsol(PREC,b) Solves Ax=b using one step of the computed ILUPACK preconditioner PREC x = AMGhsol(PREC,b) For general nonsymmetric and non-Hermitian matrices A, this routine solves A'x=b where A' is the conjugate transposed matrix w.r.t. A. This is done using one step of the computed ILUPACK preconditioner PREC x = AMGtsol(PREC,b) For general nonsymmetric and non-Hermitian matrices A, this routine solves A.'x=b, where A.' is the transposed matrix w.r.t. A. This is done using one step of the computed ILUPACK preconditioner PREC
nz=AMGnnz(PREC) total number of nonzeros of the multilevel ILU to be consistent with MATLAB, the preconditioner is treated as if it were unsymmetric
AMGspy(PREC) display multilevel preconditioner PREC AMGspy(A,PREC) display remapped original matrix A associated with the sequence of reorderings given by PREC
PREC = AMGconvert(PREC) If the preconditioner is real symmetric and indefinite or complex Hermitian and indefinite, then this routine turns PREC into a positive definite preconditioner
[A,rhs,rhstyp]=loadhbo(filename) load matrix A and optionally right hand side b, initial guess x0 and exact solution x Input ----- filename name of the file without extension Output ------ A mxn sparse matrix rhs vector of right hand side(s), initial guess(es), exact solution(s) depending on `rhstyp' rhstyp rhstyp(1)=='F' or rhstyp=='f' => dense right hand side(s) rhstyp(2)=='G' or rhstyp=='g' => dense initial guess(es) rhstyp(2)=='X' or rhstyp=='g' => dense exact solution(s) depending on rhstyp it will be clear how many columns of `rhs' correspond to the hand side(s), the initial guess(es) or the exact solution(s) savehbo(filename, A) savehbo(filename, A,b) savehbo(filename, A,b,x0) savehbo(filename, A,b,x0,x) save matrix A and optionally right hand side b, initial guess x0 and exact solution x
[pl,pr,Dl,Dr] = mwm(A) reorder and rescale a given nxn matrix A using maximum weight matching input ----- A n x n matrix output ------ pl,pr left and right permutation vectors Dl, Dr left and right scaling matrices On exit B=Dl*A(pl,pr)*Dr would refer to the reordered and rescaled system such that in theory |B(i,i)|=1 and |B(i,j)|<=1. In practice, powers of 2 are used for scaling matrices Dl,Dr to avoid rounding errors. For this reason, the scaled entries fulfill the constraint only within the range that can be achieved by the nearest power of 2 for Dl,Dr [pl,pr,Dl,Dr] = mwmamd(A) reorder and rescale a given nxn matrix A using maximum weight matching followed by approximate minimum degree [pl,pr,Dl,Dr] = mwmmetisn(A) reorder and rescale a given nxn matrix A using maximum weight matching followed by Metis Nested Dissection by nodes [pl,pr,Dl,Dr] = mwmmetise(A) reorder and rescale a given nxn matrix A using maximum weight matching followed by Metis Nested Dissection by edges [pl,pr,Dl,Dr] = mwmmmd(A) reorder and rescale a given nxn matrix A using maximum weight matching followed by Minimum Degree [pl,pr,Dl,Dr] = mwmrcm(A) reorder and rescale a given nxn matrix A using maximum weight matching followed by reverse Cuthill-McKee
[p,D] = symmwmamd(A) reorder and rescale a given nxn SYMMETRIC/HERMITIAN matrix A using symmetric maximum weight matching followed by approximate minimum degree input ----- A nxn matrix output ------ p permutation vector. On exit D*A(p,p)*D refers to the reordered and rescaled system [p,D] = symmwmmetisn(A) reorder and rescale a given nxn SYMMETRIC/HERMITIAN matrix A using symmetric maximum weight matching followed by METIS multilevel nested dissection by nodes [p,D] = symmwmmetise(A) reorder and rescale a given nxn SYMMETRIC/HERMITIAN matrix A using symmetric maximum weight matching followed by METIS multilevel nested dissection by edges [p,D] = symmwmmmd(A) reorder and rescale a given nxn SYMMETRIC/HERMITIAN matrix A using symmetric maximum weight matching followed by minimum degree [p,D] = symmwmrcm(A) reorder and rescale a given nxn SYMMETRIC/HERMITIAN matrix A using symmetric maximum weight matching followed by Reverse Cuthill-McKee
p = symmetisn(A) reorder a given nxn matrix A using METIS multilevel nested dissection by nodes input ----- A nxn matrix output ------ p permutation vector. On exit A(p,p) refers to the reordered system p = symmetise(A) reorder a given nxn matrix A using METIS multilevel nested dissection by edges
[x,flag,iter,resvec]=sqmr(A,b,tol,maxit,M1,M2,x0) sQMR simplified Quasi-Minimal Residual Method for symmetric matrices A symmetric preconditioning. x = sqmr(A,b) attempts to solve the system of linear equations A*x=b for x. The n-by-n coefficient matrix A must be square and real symmetric (resp. complex Hermitian or complex symmetric) and the right hand side column vector b must have length n. x = sqmr(A,b,tol) specifies the tolerance of the method. If tol not specified, then sqmr uses the default, 1e-6. x = sqmr(A,b,tol,maxit) specifies the maximum number of iterations. If maxit is not specified then sqmr uses the default, min(N,500). x = sqmr(A,b,tol,maxit,M) and x = sqmr(A,b,tol,maxit,M1,M2) use preconditioners M or M=M1*M2 and effectively solve the system inv(M)*A*x = inv(M)*b for x. M may be a function handle mfun such that mfun(x) returns M\x. x = sqmr(A,b,tol,maxit,M1,M2,x0) and x=sqmr(A,b,tol,maxit,M,x0) specifies the initial guess. If x0 is not specified then sqmr uses the zero vector. [x,flag] = sqmr(A,b,...) also returns a convergence flag: 0 sqmr converged to the desired tolerance tol within maxit iterations. As stopping criterion the backward error ||Ax-b||/(||A|| ||x||+||b||) is used 1 sqmr iterated maxit times but did not converge. 2 break down -m unknown error code [x,flag,iter] = sqmr(A,b,...) also returns the iteration number at which x was computed: 0 <= iter <= maxit. [x,flag,iter,resvec] = sqmr(A,b,...) also returns a vector of the estimated residual norms at each iteration, including norm(b-A*x0).