<< splspc_rifgmresab Splspc

Splspc >> Splspc > splspc_rifgmresba

splspc_rifgmresba

BA-GMRES iterative solver with a RIF preconditioner

Calling Sequence

x=splspc_rifgmresba(A, b)
x=splspc_rifgmresba(A, b, trif)        
x=splspc_rifgmresba(A, b, trif, restart)
x=splspc_rifgmresba(A, b, trif, restart, maxiter)
x=splspc_rifgmresba(A, b, trif, restart, maxiter, tol)
[x, nbiter]= splsc_rifgmresba(...)
[x, nbiter, res]= splsc_rifgmresba(...)

Parameters

A

sparse matrix (size m x n) : the matrix A in the least square problem min(b-Ax)

b

matrix of double (size m x 1) : the vector b in the least square problem min(b-Ax)

trif

matrix of double (size 1 x 1) : a dropping tolerance parameter for the eliminations of the smaller elements during the algorithm (default = 1.e-1)

restart

matrix of integer (size 1 x 1) : the number of iterations before restarting the algorithm (default = min(m,n))

maxiter

matrix of integer (size 1 x 1) : the maximum number of iterations in the algorithm (default = min(m,n))

tol

matrix of integer (size 1 x 1) : the threshold of tolerance for stopping the algorithm (default = 1.e-6)

x

matrix of double (size n x 1) : solution vector x in the least square problem min(b-Ax)

nbiter

matrix of integer (size 1 x 1) : total number of necessary iterations to find the solution vector

res

matrix of double (size 1 x 1) : norm of the computed residual r=||A'*(b-A*x)||/||A'*b||

Description

This function solves the sparse linear least squares problem min |Ax-b|, with the BA-GMRES method with RIF preconditionning.

Any optional argument equal to the empty matrix [] is replaced by its default value.

If the algorithm does not converge, a warning is produced.

The GMRES method consist in solving the least square problem min(b-Ax), for A square matrix. But when A is not square, GMRES can't be directly used and need to be slightly modified.

The BA-GMRES method consist in solving the least square problem min(Bb-BAx), ie finding x such as the norm of Bb-BAx is minimal. This method is mainly useful for overdetermined problems, ie m > n for A of size m x n, because the matrix obtained by doing the product of B and A is of size n x n.

The Robust Incomplete Factorization (RIF) is a method applied before an iterative method (such as GMRES or CGLS). Depending on the larger dimension of the matrix A, this methods seeks to calculate an approximation of A^(T)A or AA^(T). Indeed, it gives an upper triangular matrix Z such as A^(T)A (or AA^(T)) is almost equal to (Z^-T)D(Z^-1) with D diagonal. Thus, the preconditioner is computed in this way : if m > n, we have C={Z(D^-1)(Z^T)} and the preconditioner is B=C(A^T). Alternatively, if m < n, we have C={Z(D^-1)(Z^T)} and the preconditioner is B=(A^T)C. The accuracy and speed of the method is strongly dependant of the choice of the dropping tolerance parameter. If m < n, A^T is automatically used rather than A to calculate the preconditioner and then B^T is taken as a preconditioner in the iterative method.

Restart can be used to save memory but is not adviced in this case. The initial vector x0 for the iterative process has been set up to a vector of zeros.

Examples

path = fullfile(splspc_getpath(),"tests","matrices");
filename=fullfile(path,"illc1033.mtx");
A = mmread(filename);
b=ones(size(A,1),1);
// With default settings
x=splspc_rifgmresba(A,b);
norm(A*x-b)/norm(b)
// Configure the options
trif=1.e-1;
restart=1000;
max_iter=3000;
tol=1.e-6;
[x,nbiter,res]=splspc_rifgmresba(A,b,trif,restart,max_iter,tol);
// Configure only max_iter
x=splspc_rifgmresba(A,b,[],[],max_iter);

Authors

1994 - 2008 - Sparselib++ v1.7 - R. Pozo, K. Remington, and A. Lumsdaine

2010 - National Institute of Informatics - K. Hayami, T. Ito and J.-F. Yin

2011 - DIGITEO - Michael Baudin

2011 - National Institute of Informatics - Benoit Goepfert

See Also

<< splspc_rifgmresab Splspc