ILUPACK toolbox for MATLAB

  1. Getting started
  2. Trouble shooting
  3. What's behind the toolbox
  4. Overview options
  5. Tools
  6. Automatic structure detection
  7. Function reference list

Getting started

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.


Trouble shooting

o
Your iterative solver does not converge
This happens if an error message is displayed on exit from
AMGsolver.

reduce options.droptol, delete preconditioner by

>> 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);



o
When calling AMGfactor, ILUPACK tells you that your memory is insufficient
Starting with ILUPACK V2.2 this error message should no longer occur. ILUPACK should automatically adapt param.elbow appropriately. Please report failures to us.



o
For the iterative solver you want to provide your own initial guess

Call AMGsolver as follows:
>> [PREC,options]=AMGsolver(A,PREC,options,b,x0);



o
The iterative solver AMGsolver terminates with an error message telling you that it did not converge within a certain number of steps, BUT you would accept more iteration steps.

increase options.maxit to your desired maximum number of iteration steps before using AMGsolver again.
In this case you could recycle a former computed solution, say x, from a previous call of AMGsolver and use it as initial guess for the next run, e.g.

>> options.maxit=10000;
>> [x,options]=AMGsolver(A, PREC, options, b,x);



o
You would like to have a more accurate solution

reduce options.restol and restart the iterative solver AMGsolver.
In this case you could recycle a former computed solution, say x, from a previous call of AMGsolver and use it as initial guess for the next run, e.g.

>> options.restol=1e-14;
>> [x,options]=AMGsolver(A, PREC, options,b,x);

What's behind the toolbox

The main philosophy of ILUPACK can be summarized in three major components.

  1. Maximum weight matching to improve the (block) diagonal dominance of the system
  2. Symmetric reorderings to reduce the fill independent on the drop tolerance
  3. Inverse-based pivoting to keep the error propagation small that is caused by dropping small entries
  4. multilevel approach as wrapper around these three aspects

Maximum weight matchings


Symmetric reorderings

ILUPACK supports several well-estabished symmetric reorderings, e.g.

Inverse-based pivoting

Watch the impact of using a different k in the multilevel method. Roughly speaking, a small k will produce more levels but you obtain more robustness and possibly less fill-in.

Overview options

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.

Tools

Beside the four routines that already have been described, ILUPACK offers some handy additional functions.

Automatic structure detection

ILUPACK offers many special purpose drivers for The ILUPACK toolbox for MATLAB automatically detects

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);


Function reference list

The following functions can be used as part of the ILUPACK toolbox.
  1. AMGinit
  2. AMGfactor
  3. AMGsolver, AMGhsolver, AMGtsolver
  4. AMGdelete
  5. AMGsol, AMGhsol, AMGtsol
  6. AMGnnz
  7. AMGspy
  8. AMGconvert
  9. loadhbo / savehbo
  10. mwm / mwmamd / mwmmetisn / mwmmetise / mwmmmd / mwmrcm
  11. symmwmamd / symmwmmetisn / symmwmmetise / symmwmmmd / symmwmrcm
  12. symmetisn / symmetise
  13. sqmr



  1. AMGinit

     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
    

  2. AMGfactor

     [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
    

  3. AMGsolver, AMGhsolver, AMGtsolver

    [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
    

  4. AMGdelete

     PREC=AMGdelete(PREC)
     
     delete preconditioner, in particular release the associated memory
    

  5. AMGsol, AMGhsol, AMGtsol

    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
    

  6. AMGnnz

     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
    

  7. AMGspy

     AMGspy(PREC)
    
     display multilevel preconditioner PREC
    
     AMGspy(A,PREC)
    
     display remapped original matrix A associated with the sequence of 
     reorderings given by PREC
    

  8. AMGconvert

     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
    

  9. loadhbo / savehbo

     [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
    
    

  10. mwm / mwmamd / mwmmetisn / mwmmetise / mwmmmd / mwmrcm

    [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
    

  11. symmwmamd / symmwmmetisn / symmwmmetise / symmwmmmd / symmwmrcm

     [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
    

  12. symmetisn / symmetise

    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
    

  13. sqmr

    [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).
    


m.bollhoefer@tu-bs.de

Last modified: February 10, 2012