<< splspc_gmresba Splspc splspc_grevgmresba >>

Splspc >> Splspc > splspc_grevgmresab

splspc_grevgmresab

AB-GMRES iterative solver with a GREVILLE preconditioner

Calling Sequence

x=splspc_grevgmresab(A, b)
x=splspc_grevgmresab(A, b, tgrev)       
x=splspc_grevgmresab(A, b, tgrev, tsgrev)        
x=splspc_grevgmresab(A, b, tgrev, tsgrev, restart)
x=splspc_grevgmresab(A, b, tgrev, tsgrev, restart, maxiter)
x=splspc_grevgmresab(A, b, tgrev, tsgrev, restart, maxiter, tol)
[x, nbiter]= splsc_grevgmresab(...)
[x, nbiter, res]= splsc_grevgmresab(...)

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)

tgrev

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

tsgrev

matrix of double (size 1 x 1) : a switching tolerance parameter for the detection of linearly dependants colums in the matrix A (default = 1.e-6)

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'*A*x

Description

This function solves the sparse linear least squares problem min |Ax-b|, with the AB-GMRES method with Greville's 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 AB-GMRES method consist in solving the least square problem min(b-ABx), ie finding x such as the norm of b-ABx is minimal. This method is mainly useful for underdetermined 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 m x m.

The Greville's method is a preconditioning technique applied before an iterative method (such as GMRES or CGLS) for solving least square problems. This methods seeks to calculate an approximation of A^(-1). Indeed, it gives an upper triangular matrix Z, a diagonal matrix D and a matrix V such as A^-1 ~= Z(D^-1)V. If A is full column rank, then the implementation of the preconditioner B is as follow : we have C=Z(D^-1)Z and B=C(A^T). When A is not full column rank, a switching tolerance parameter is used to detect with a chosen accuracy the linear dependency of the columns of A. The accuracy and speed of the method are strongly dependant of the choice of the dropping and switching tolerance parameters. 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,"lp_scfxm1.mtx");
A = mmread(filename);
A = A';
e = ones(size(A,2),1);
b = A * e;
// With default settings
x=splspc_grevgmresab(A,b);
norm(A*x-b)/norm(b)
// Configure the options
tgrev=1.e-3;
tsgrev=1.e-6;
restart=1000;
max_iter=3000;
tol=1.e-6;
[x,nbiter,res]=splspc_grevgmresab(A,b,tgrev,tsgrev,restart,max_iter,tol);
// Configure only max_iter
x=splspc_grevgmresab(A,b,[],[],[],max_iter);

Authors

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

2011 - National Institute of Informatics - K. Hayami, X. Cui and J.-F. Yin

2011 - DIGITEO - Michael Baudin

2011 - National Institute of Informatics - Benoit Goepfert

See Also

<< splspc_gmresba Splspc splspc_grevgmresba >>