diff options
Diffstat (limited to 'modules/sparse/macros/gmres.sci')
| -rwxr-xr-x | modules/sparse/macros/gmres.sci | 308 |
1 files changed, 308 insertions, 0 deletions
diff --git a/modules/sparse/macros/gmres.sci b/modules/sparse/macros/gmres.sci new file mode 100755 index 000000000..ec0943042 --- /dev/null +++ b/modules/sparse/macros/gmres.sci @@ -0,0 +1,308 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) XXXX-2008 - INRIA +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +// [x, flag, resNorm, iter, resVec] = gmres( A, b, x, M, restrt, max_it, tol ) +// +// GMRES solves the linear system Ax=b +// using the Generalized Minimal RESidual ( GMRES ) method with restarts . +// +// input A REAL nonsymmetric positive definite matrix or function +// x REAL initial guess vector +// b REAL right hand side vector +// M REAL preconditioner matrix or function +// restrt INTEGER number of iterations between restarts +// max_it INTEGER maximum number of iterations +// tol REAL error tolerance +// +// output x REAL solution vector +// flag INTEGER: 0 = solution found to tolerance +// 1 = no convergence given max_it +// resNorm REAL final residual norm +// iter INTEGER number of iterations performed +// resVec REAL residual vector + +// Details of this algorithm are described in +// +// "Templates for the Solution of Linear Systems: Building Blocks +// for Iterative Methods", +// Barrett, Berry, Chan, Demmel, Donato, Dongarra, Eijkhout, +// Pozo, Romine, and Van der Vorst, SIAM Publications, 1993 +// (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps). +// +// "Iterative Methods for Sparse Linear Systems, Second Edition" +// Saad, SIAM Publications, 2003 +// (ftp ftp.cs.umn.edu; cd dept/users/saad/PS; get all_ps.zip). + +function [x, flag, resNorm, iter, resVec] = gmres(A, varargin) + + // ----------------------- + // Parsing input arguments + // ----------------------- + + [lhs,rhs]=argn(0); + if ( rhs < 2 ), + error(msprintf(gettext("%s: Wrong number of input argument: At least %d expected.\n"),"gmres",2)); + end + + // Parsing the matrix A et the right hand side vector b + select type(A) + case 1 then + matrixType = 1; + case 5 then + matrixType = 1; + case 13 then + matrixType = 0; + end + // If A is a matrix (full or sparse) + if (matrixType == 1), + if (size(A,1) ~= size(A,2)), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Square matrix expected.\n"),"gmres",1)); + end + end + b=varargin(1); + if (size(b,2) ~= 1), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Column vector expected.\n"),"gmres",2)); + end + if (matrixType==1), + if (size(b,1) ~= size(A,1)), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"gmres",2,1)); + end + end + + // Number of iterations between restarts + if (rhs >= 3), + restrt=varargin(2); + if (size(restrt) ~= [1 1]), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"),"gmres",3)); + end + else + restrt=20; + end + + // Error tolerance tol + if (rhs >= 4), + tol=varargin(3); + if (size(tol) ~= [1 1]); + error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"),"gmres",4)); + end + else + tol = 1e-6; + end + + // Maximum number of iterations max_it + if (rhs >= 5), + max_it=varargin(4); + if (size(max_it) ~= [1 1]), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"),"gmres",5)); + end + else + max_it=size(b,1); + end + + // Parsing of the preconditioner matrix M + if (rhs >= 6), + M = varargin(5); + select type(M) + case 1 then + precondType = 1; + case 5 then + precondType = 1; + case 13 then + precondType = 0; + end + if (precondType == 1), + if (size(M,1) ~= size(M,2)), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Square matrix expected.\n"),"gmres",4)); + end + if (size(M,1) == 0), + precondType = 2; // no preconditionning + elseif ( size(M,1) ~= size(b,1) ), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"gmres",4,2)); + end + end + else + precondType = 2; // no preconditionning + end + + // Parsing of the initial vector x + if (rhs >= 7), + x=varargin(6); + if (size(x,2) ~= 1), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Column vector expected.\n"),"gmres",3)); + end + if ( size(x,1) ~= size(b,1) ), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"gmres",3,2)); + end + else + x=zeros(b); + end + + if (rhs > 7), + error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"gmres",2,7)); + end + + // ------------ + // Computations + // ------------ + + j = 0; + flag = 0; + it2 = 0; + + bnrm2 = norm(b); + if (bnrm2 == 0.0), + x = zeros(b); + resNorm = 0; + iter = 0; + resVec = resNorm; + flag = 0; + return + end + + // r = M \ ( b-A*x ); + if (matrixType == 1), + r = b - A*x; + else + r = b - A(x); + end + if (precondType == 1), + r = M \ r; + elseif (precondType == 0), + r = M(r); + end + resNorm = norm(r)/bnrm2; + resVec = resNorm; + if (resNorm < tol), + iter=0; + return; + end + + n = size(b,1); + m = restrt; + V(1:n,1:m+1) = zeros(n,m+1); + H(1:m+1,1:m) = zeros(m+1,m); + cs(1:m) = zeros(m,1); + sn(1:m) = zeros(m,1); + e1 = zeros(n,1); + e1(1) = 1.0; + + for j = 1:max_it + // r = M \ ( b-A*x ); + if (matrixType == 1), + r = b - A*x; + else + r = b - A(x); + end + if (precondType == 1), + r = M \ r; + elseif (precondType == 0), + r = M(r); + end + + V(:,1) = r / norm( r ); + s = norm( r )*e1; + for i = 1:m // construct orthonormal + it2 = it2 + 1; // basis using Gram-Schmidt + // w = M \ (A*V(:,i)); + if (matrixType == 1), + w = A*V(:,i); + else + w = A(V(:,i)); + end + if (precondType == 1), + w = M \ w; + elseif (precondType == 0), + w = M(w); + end + + for k = 1:i + H(k,i)= w'*V(:,k); + w = w - H(k,i)*V(:,k); + end + H(i+1,i) = norm( w ); + V(:,i+1) = w / H(i+1,i); + for k = 1:i-1 // apply Givens rotation + temp = cs(k)*H(k,i) + sn(k)*H(k+1,i); + H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i); + H(k,i) = temp; + end + // form i-th rotation matrix + [tp1,tp2] = rotmat( H(i,i), H(i+1,i) ); + cs(i) = tp1; + sn(i) = tp2; + temp = cs(i)*s(i); + s(i+1) = -sn(i)*s(i); + s(i) = temp; + H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i); + H(i+1,i) = 0.0; + resNorm = real(abs(s(i+1))) / bnrm2; + resVec = [resVec;resNorm]; + if ( resNorm <= tol ), + y = H(1:i,1:i) \ s(1:i); + x = x + V(:,1:i)*y; + break; + end + end + if (resNorm <= tol), + iter = j-1+it2; + break; + end + y = H(1:m,1:m) \ s(1:m); + // update approximation + x = x + V(:,1:m)*y; + // r = M \ ( b-A*x ) + if (matrixType == 1), + r = b - A*x; + else + r = b - A(x); + end + if (precondType == 1), + r = M \ r; + elseif (precondType == 0), + r = M(r); + end + s(j+1) = norm(r); + resNorm = real(s(j+1)) / bnrm2; + resVec = [resVec; resNorm]; + + if ( resNorm <= tol ), + iter = j+it2; + break; + end + if ( j== max_it ), + iter=j+it2; + end + end + if ( resNorm > tol ), + flag = 1; + if (lhs < 2), + warning(msprintf(gettext("%s: Did not converge.\n"),"gmres")); + end + end +endfunction //GMRES + + +// +// Compute the Givens rotation matrix parameters for a and b. +// +function [ c, s ] = rotmat( a, b ) + if ( b == 0.0 ), + c = 1.0; + s = 0.0; + elseif ( abs(b) > abs(a) ), + temp = a / b; + s = 1.0 / sqrt( 1.0 + temp^2 ); + c = temp * s; + else + temp = b / a; + c = 1.0 / sqrt( 1.0 + temp^2 ); + s = temp * c; + end +endfunction //rotmat |