diff options
Diffstat (limited to 'modules/sparse/macros/%bicgstab.sci')
| -rwxr-xr-x | modules/sparse/macros/%bicgstab.sci | 202 |
1 files changed, 202 insertions, 0 deletions
diff --git a/modules/sparse/macros/%bicgstab.sci b/modules/sparse/macros/%bicgstab.sci new file mode 100755 index 000000000..fe90693ed --- /dev/null +++ b/modules/sparse/macros/%bicgstab.sci @@ -0,0 +1,202 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier +// +// 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 + +// +// bicgstab -- +// BICG solves the linear system %Ax=b using the BiConjugate Gradient Stabilized method. +// If M is given, it is used as a preconditionning matrix. +// If both M and M2 are given, the matrix M * M2 is used as a preconditionning +// matrix. +// +// input %A REAL matrix or a function y=Ax(x) which computes y=%A*x for a given x +// b REAL right hand side vector +// tol, optional REAL error tolerance (default: 1e-8) +// maxIter, optional INTEGER maximum number of iterations (default: size(%b)) +// %M, optional REAL preconditioner matrix (default: none) +// %M2, optional REAL preconditioner matrix (default: none) +// x0, optional REAL initial guess vector (default: the zero vector) +// verbose, optional INTEGER set to 1 to enable verbose logging (default : 0) +// +// output x REAL solution vector +// resNorm REAL final relative norm of the residual +// iter INTEGER number of iterations performed +// resVec REAL residual vector +// +// References +// +// "Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems" +// by Henk van der Vorst. +// +// http://epubs.siam.org/doi/abs/10.1137/0913035 +// http://dl.acm.org/citation.cfm?id=131916.131930&coll=DL&dl=GUIDE&CFID=372773884&CFTOKEN=56630250 +// http://mathworld.wolfram.com/BiconjugateGradientStabilizedMethod.html +// +// Notes +// This script was originally a matlab > scilab translation of the bicgstab.m +// script from http://www.netlib.org/templates/matlab +// +// The input / output arguments of this command are the same as Matlab's bicgstab command. +// + +function [x, resNorm, iter, resVec] = %bicgstab(%A, %b, tol, maxIter, %M, %M2, x0, verbose ) + + // Initialization + bnrm2 = norm(%b); + if (verbose==1) then + printf(gettext("Norm of right-hand side : %s\n"), string(bnrm2)); + end + if (bnrm2 == 0) then + if (verbose==1) then + printf(gettext("Special processing where the right-hand side is zero.\n")); + end + // When rhs is 0, there is a trivial solution : x=0 + x = zeros(%b); + resNorm = 0; + resVec = resNorm; + else + x = x0; + // r = %b - %A*x; + if (matrixType ==1), + r = %b - %A*x; + r2 = r; + else + r = %b - %A(x,Aargs(:)); + r2 = r; + end + resNorm = norm(r) / bnrm2; + resVec = resNorm; + end + if (verbose==1) then + printf(gettext(" Type of preconditionning #1 : %d\n"),precondType); + printf(gettext(" Type of preconditionning #2 : %d\n"),precondBis); + end + // Begin iteration + // Distinguish the number of iterations processed from the currentiter index + iter = 0 + for currentiter = 1:maxIter + if (resNorm <= tol) then + if (verbose==1) then + printf(gettext(" New residual = %s < tol = %s => break\n"),string(resNorm),string(tol)); + end + break; + end + iter = iter + 1 + if (verbose==1) then + printf(gettext(" Iteration #%s/%s residual : %s\n"),string(currentiter),string(maxIter),string(resNorm)); + printf(" x=\n"); + disp(x); + end + rho = r2'*r; + if (rho == 0) then + break; + end + if (currentiter > 1) then + bet = (rho/rho_old)*(alp/ome); + p = r + bet*(p-ome*v); + else + p = r; + end + // Solve M' M2' P = p + if %M == [] & %M2 == [] then + P = p; + elseif %M2 == [] then + // Compute r so that M' P = p + if (precondType == 1) then + P = %M \ p; + elseif (precondType == 2) then + P = %M(p,Margs(:)); + else + P = p; + end + else + // Compute r so that M' M2' P = p + if (precondBis == 1) then + P = %M \ p; + P = %M2 \ p; + elseif (precondBis == 2) then + P = %M(p,Margs(:)); + P = %M2(p,M2args(:)); + else + P = p; + end + end + // v = %A*P; + if (matrixType == 1), + v = %A*P; + else + v = %A(P); + end + alp = rho / (r2'*v); + s = r - alp*v; + // Check for convergence + if (norm(s) <= tol) then + x = x+alp*P; + resNorm = norm(s) / bnrm2; + if (verbose==1) then + printf(gettext(" New residual = %s < tol = %s => break\n"),string(resNorm),string(tol)); + end + resVec = [resVec;resNorm]; + break; + end + // Solve M M2 S = s + if %M == [] & %M2 == [] then + S = s; + elseif %M2 == [] then + // Compute r so that M S = s + if (precondType == 1) then + S = %M \ s; + elseif (precondType == 2) then + S = %M(s,Margs(:)); + else + S = s; + end + else + // Compute r so that M M2 S = s + if (precondBis == 1) then + S = %M \ s; + S = %M2 \ s; + elseif (precondBis == 2) then + S = %M(s,Margs(:)); + S = %M2(s,M2args(:)); + else + S = s; + end + end + // t = %A*S; + if (matrixType == 1), + t = %A*S; + else + t = %A(S); + end + ome = (t'*s)/(t'*t); + x = x + alp*P+ome*S; + r = s - ome*t; + resNorm = norm(r) / bnrm2; + // Caution : transform the scalar resVec into vector resVec ! + resVec = [resVec;resNorm]; + rho_old = rho; + end + // Test for convergence + if (resNorm > tol) then + if (verbose==1) then + printf(gettext("Final residual = %s > tol =%s\n"),string(resNorm),string(tol)); + printf(gettext("Algorithm fails\n")); + end + flag = 1; + if (lhs < 2) then + warning(msprintf(gettext("%s: Convergence error.\n"),"bicgstab")); + end + else + flag = 0; + if (verbose==1) then + printf(gettext("Algorithm pass\n")); + end + end + +endfunction |