diff options
Diffstat (limited to 'modules/sparse/macros')
35 files changed, 2295 insertions, 0 deletions
diff --git a/modules/sparse/macros/%bicg.bin b/modules/sparse/macros/%bicg.bin Binary files differnew file mode 100755 index 000000000..ce2399d50 --- /dev/null +++ b/modules/sparse/macros/%bicg.bin diff --git a/modules/sparse/macros/%bicg.sci b/modules/sparse/macros/%bicg.sci new file mode 100755 index 000000000..800fd319c --- /dev/null +++ b/modules/sparse/macros/%bicg.sci @@ -0,0 +1,189 @@ +// 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 + +// +// bicg -- +// BICG solves the linear system %Ax=b using the BiConjugate Gradient 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 +// +// "Numerical Recipes: The Art of Scientific Computing." (third ed.) +// by William Press, Saul Teukolsky, William Vetterling, Brian Flannery. +// +// http://apps.nrbook.com/empanel/index.html?pg=87 +// http://dl.acm.org/citation.cfm?doid=1874391.187410 +// http://mathworld.wolfram.com/BiconjugateGradientMethod.html +// +// Notes +// This script was originally a matlab > scilab translation of the bicg.m +// script from http://www.netlib.org/templates/matlab +// +// The input / output arguments of this command are the same as Matlab's bicg command. +// + +function [x, resNorm, iter, resVec] = %bicg(%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 + // Solve M M2 r = r + if %M == [] & %M2 == [] then + z = r; + elseif %M2 == [] then + // Compute r so that M r = r + if (precondType == 1) then + z = %M \ r; + elseif (precondType == 2) then + z = %M(r,Margs(:)); + else + z = r; + end + else + // Compute r so that M M2 r = r + if (precondBis == 1) then + z = %M \ r; + z = %M2 \ r; + elseif (precondBis == 2) then + z = %M(r,Margs(:)); + z = %M2(r,M2args(:)); + else + z = r; + end + end + // Solve M' M2' r2 = r2 + if %M == [] & %M2 == [] then + z2 = r2; + elseif %M2 == [] then + // Compute r so that M' r = r + if (precondType == 1) then + z2 = %M' \ r2; + elseif (precondType == 2) then + z2 = %M(r2,Margs(:)); + else + z2 = r2; + end + else + // Compute r so that M' M2' r = r + if (precondBis == 1) then + z2 = %M' \ r2; + z2 = %M2' \ r2; + elseif (precondBis == 2) then + z2 = %M(r2,Margs(:)); + z2 = %M2(r2,M2args(:)); + else + z2 = r2; + end + end + rho = r'*z2; + if (rho == 0) then + break; + end + if (currentiter > 1) then + bet = rho / rho_old; + p = z + bet*p; + p2 = z2 + bet*p2; + else + p = z; + p2 = z2; + end + // q = %A*p; q2 = %A'*p2; + if (matrixType == 1), + q = %A*p; + q2 = %A'*p2; + else + q = %A(p); + q2 = %A(p2); + end + alp = rho / (p2'*q); + x = x + alp*p; + r = r - alp*q; + r2 = r2 - alp*q2; + 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"),"bicg")); + end + else + flag = 0; + if (verbose==1) then + printf(gettext("Algorithm pass\n")); + end + end + +endfunction diff --git a/modules/sparse/macros/%bicgstab.bin b/modules/sparse/macros/%bicgstab.bin Binary files differnew file mode 100755 index 000000000..51b70eaaf --- /dev/null +++ b/modules/sparse/macros/%bicgstab.bin 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 diff --git a/modules/sparse/macros/%cgs.bin b/modules/sparse/macros/%cgs.bin Binary files differnew file mode 100755 index 000000000..fee108f61 --- /dev/null +++ b/modules/sparse/macros/%cgs.bin diff --git a/modules/sparse/macros/%cgs.sci b/modules/sparse/macros/%cgs.sci new file mode 100755 index 000000000..bec33e08b --- /dev/null +++ b/modules/sparse/macros/%cgs.sci @@ -0,0 +1,194 @@ +// 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 + +// +// cgs -- +// CGS solves the linear system %Ax=b using the Conjugate Gradient Squared 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 +// +// "CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems" +// by Peter Sonneveld +// +// http://epubs.siam.org/doi/abs/10.1137/0910004 +// http://dl.acm.org/citation.cfm?id=64888&preflayout=flat +// http://mathworld.wolfram.com/ConjugateGradientSquaredMethod.html +// +// Notes +// This script was originally a matlab > scilab translation of the cgs.m +// script from http://www.netlib.org/templates/matlab +// +// The input / output arguments of this command are the same as Matlab's cgs command. +// + +function [x, resNorm, iter, resVec] = %cgs(%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; + u = r + bet*q; + p = u + bet*(q+bet*p); + else + u = r; + p = u; + end + // Solve M M2 P = p + if %M == [] & %M2 == [] then + P = p; + elseif %M2 == [] then + // Compute P 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 P 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); + q = u - (alp*v); + // Solve M M2 u = u+q + uq = u + q; + if %M == [] & %M2 == [] then + U = uq; + elseif %M2 == [] then + // Compute Q so that M U = u+q + if (precondType == 1) then + U = %M \ uq; + elseif (precondType == 2) then + U = %M(uq,Margs(:)); + else + U = uq; + end + else + // Compute z so that M M2 U = u+q + if (precondBis == 1) then + U = %M \ uq; + U = %M2 \ uq; + elseif (precondBis == 2) then + U = %M(uq,Margs(:)); + U = %M2(uq,M2args(:)); + else + U = uq; + end + end + x = x + alp*U; + // U = %A*U; + if (matrixType ==1), + U = %A*U; + else + U = %A(U); + end + r = r - alp*U; + 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"),"cgs")); + end + else + flag = 0; + if (verbose==1) then + printf(gettext("Algorithm pass\n")); + end + end + +endfunction diff --git a/modules/sparse/macros/%pcg.bin b/modules/sparse/macros/%pcg.bin Binary files differnew file mode 100755 index 000000000..b36f6439d --- /dev/null +++ b/modules/sparse/macros/%pcg.bin diff --git a/modules/sparse/macros/%pcg.sci b/modules/sparse/macros/%pcg.sci new file mode 100755 index 000000000..be3162f87 --- /dev/null +++ b/modules/sparse/macros/%pcg.sci @@ -0,0 +1,149 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2009 - INRIA - Michael Baudin +// Copyright (C) 2008 - INRIA - Michael Baudin +// Copyright (C) 2006 - INRIA - Serge Steer +// Copyright (C) 2005 - IRISA - Sage Group +// +// 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 + +// +// pcg -- +// PCG solves the symmetric positive definite linear system %Ax=b +// using the Preconditionned Conjugate Gradient. +// 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 symmetric positive definite 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 +// +// "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). +// +// Golub and Van Loan, Matrix Computations +// +// Notes +// This script was originally a matlab > scilab translation of the cg.m +// script from http://www.netlib.org/templates/matlab +// +// The input / output arguments of this command are the same as +// Matlab's cg command. +// + +function [x, resNorm, iter, resVec] = %pcg(%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; + else + r = %b - %A(x,Aargs(:)); + 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 + if %M == [] & %M2 == [] then + z = r; + elseif %M2 == [] then + // Compute z so that M z = r + if (precondType == 1) then + z = %M \ r; + elseif (precondType == 2) then + z = %M(r,Margs(:)); + else + z = r; + end + else + // Compute z so that M M2 z = r + if (precondBis == 1) then + z = %M \ r; + z = %M2 \ z; + elseif (precondBis == 2) then + z = %M(r,Margs(:)); + z = %M2(z,M2args(:)); + else + z = r; + end + end + rho = r'*z; + if (currentiter > 1) then + bet = rho / rho_old; + p = z + bet*p; + else + p = z; + end + // q = %A*p; + if (matrixType ==1), + q = %A*p; + else + q = %A(p); + end + alp = rho / (p'*q ); + x = x + alp*p; + r = r - alp*q; + resNorm = norm(r) / bnrm2; + // Caution : transform the scalar resVec into vector resVec ! + resVec = [resVec;resNorm]; + rho_old = rho; + end + +endfunction diff --git a/modules/sparse/macros/adj2sp.bin b/modules/sparse/macros/adj2sp.bin Binary files differnew file mode 100755 index 000000000..3a6cceefc --- /dev/null +++ b/modules/sparse/macros/adj2sp.bin diff --git a/modules/sparse/macros/adj2sp.sci b/modules/sparse/macros/adj2sp.sci new file mode 100755 index 000000000..7a2b66085 --- /dev/null +++ b/modules/sparse/macros/adj2sp.sci @@ -0,0 +1,97 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) XXXX-2008 - INRIA +// Copyright (C) 2010 - DIGITEO - Michael Baudin +// +// 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 + + +function A=adj2sp(varargin) + // adjacency to sparse conversion. + // See adj2sp.xml for help. + // + // Check number of arguments + [lhs, rhs] = argn() + if ( and(rhs <> [3 4]) ) then + lstr = gettext("%s: Wrong number of input arguments: %d to %d expected.\n") + error ( msprintf(lstr,"adj2sp",3,4)) + end + if ( and(lhs <> [0 1]) ) then + lstr = gettext("%s: Wrong number of output arguments: %d to %d expected.\n") + error ( msprintf(lstr,"adj2sp",0,1)) + end + // + // Get arguments + xadj = varargin(1) + iadj = varargin(2) + v = varargin(3) + // + // Get size of the matrix + if ( rhs == 3 ) then + m = max(iadj) + n = size(xadj,"*")-1 + mn = [m,n] + else + mn = varargin(4) + end + // + // Check type of arguments + if ( typeof(xadj) <> "constant" ) then + lstr = gettext("%s: Wrong type for input argument #%d.\n") + error ( msprintf(lstr,"adj2sp",1)) + end + if ( typeof(iadj) <> "constant" ) then + lstr = gettext("%s: Wrong type for input argument #%d.\n") + error ( msprintf(lstr,"adj2sp",2)) + end + if ( typeof(mn) <> "constant" ) then + lstr = gettext("%s: Wrong type for input argument #%d.\n") + error ( msprintf(lstr,"adj2sp",3)) + end + // + // Check size of arguments + if ( size(mn,"*") <> 2 ) then + lstr = gettext("%s: Wrong size for input argument #%d: %s has size %d, but %d is expected.\n") + error ( msprintf(lstr,"adj2sp",4,"mn",size(mn,"*"),2)) + end + // + // Check content of arguments + if ( mn(1) < 0 ) then + lstr = gettext("%s: Wrong value for input argument #%d: %s=%d, but a positive value is expected.\n") + error ( msprintf(lstr,"adj2sp",4,"mn(1)",mn(1))) + end + if ( mn(2) < 0 ) then + lstr = gettext("%s: Wrong value for input argument #%d: %s=%d, but a positive value is expected.\n") + error ( msprintf(lstr,"adj2sp",4,"mn(2)",mn(2))) + end + if ( int(mn(1)) <> mn(1) ) then + lstr = gettext("%s: Wrong value for input argument #%d: %s=%s, but a floating point integer is expected.\n") + error ( msprintf(lstr,"adj2sp",4,"mn(1)",string(mn(1)))) + end + if ( int(mn(2)) <> mn(2) ) then + lstr = gettext("%s: Wrong value for input argument #%d: %s=%s, but a floating point integer is expected.\n") + error ( msprintf(lstr,"adj2sp",4,"mn(2)",string(mn(2)))) + end + if ( mn(1) <> max(iadj) ) then + lstr = gettext("%s: Wrong value for input argument #%d: %s=%d, which does not match %s: %d.\n") + error ( msprintf(lstr,"adj2sp",4,"mn(1)",mn(1),"max(iadj)",max(iadj))) + end + if ( mn(2) <> size(xadj,"*")-1 ) then + lstr = gettext("%s: Wrong value for input argument #%d: %s=%d, which does not match size of %s: %d.\n") + error ( msprintf(lstr,"adj2sp",4,"mn(2)",mn(2),"xadj",size(xadj,"*"))) + end + // + // Make xadj, iadj, mn column matrices + xadj = xadj(:) + iadj = iadj(:) + mn = mn(:) + // + // Proceed... + nza = size(iadj,"*") + j = fadj2sp(xadj,mn(2),nza) + A = sparse([j,iadj],v,[mn(2),mn(1)])' +endfunction + diff --git a/modules/sparse/macros/buildmacros.bat b/modules/sparse/macros/buildmacros.bat new file mode 100755 index 000000000..5df334e1e --- /dev/null +++ b/modules/sparse/macros/buildmacros.bat @@ -0,0 +1,11 @@ +rem Scilab ( http://mwww.scilab.org/ ) - This file is part of Scilab +rem Copyright (C) INRIA +rem +rem This file must be used under the terms of the CeCILL. +rem This source file is licensed as described in the file COPYING, which +rem you should have received as part of this distribution. The terms +rem are also available at +rem http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +@..\..\..\bin\scilex -nwni -ns -e exec('buildmacros.sce');quit;
\ No newline at end of file diff --git a/modules/sparse/macros/buildmacros.sce b/modules/sparse/macros/buildmacros.sce new file mode 100755 index 000000000..1ff122048 --- /dev/null +++ b/modules/sparse/macros/buildmacros.sce @@ -0,0 +1,15 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2005-2008 - INRIA - Allan CORNET +// +// 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 + +if (isdef("genlib") == %f) then + exec(SCI+"/modules/functions/scripts/buildmacros/loadgenlib.sce"); +end +//------------------------------------ +genlib("sparselib","SCI/modules/sparse/macros",%f,%t); +//------------------------------------ diff --git a/modules/sparse/macros/chfact.bin b/modules/sparse/macros/chfact.bin Binary files differnew file mode 100755 index 000000000..bd0766b0f --- /dev/null +++ b/modules/sparse/macros/chfact.bin diff --git a/modules/sparse/macros/chfact.sci b/modules/sparse/macros/chfact.sci new file mode 100755 index 000000000..34e05d8e3 --- /dev/null +++ b/modules/sparse/macros/chfact.sci @@ -0,0 +1,132 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-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 + +function spcho=chfact(A) + //cholesky factors, returned in a tlist + //spcho = {xlnz, nnzl, xsuper, xlindx, lindx, snode, + // split, tmpsiz, perm, invp, lnz}. + // + // invp=spcho('invp');xlnz=spcho('xlnz'); + // xlindx=spcho('xlindx');lindx=spcho('lindx');lnz=spcho('lnz'); + // P=sparse([(1:m)',invp],ones(invp),[m,m]); + // adjncy = spcompack(xlnz,xlindx,lindx); + // R=adj2sp(xlnz,adjncy,lnz); + // =>P*R*R'*P'=A + m = size(A,1); + if size(A,2) ~= m | type(A)~=5 | A == [] then, + error("Matrix must be square, sparse and nonempty."); + end; + neqns=m; + [xadj,adjncy,v]=sp2adj(A-diag(diag(A))); + [perm,invp,nofsub]=ordmmd(xadj,adjncy,neqns); + cachsz = 16; + spcho=symfct(xadj,adjncy,perm,invp,cachsz,neqns); + [xadjf,adjncyf,v]=sp2adj(A); + spcho=inpnv(xadjf,adjncyf,v,spcho); + level=8; + spcho=blkfc1(spcho,level); +endfunction + +function [spcho]= blkfc1(spcho,level) + //retrieves Fortran variables (see sfinit.f,bfinit.f,symfct.f ) + //[xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp,lnz]=spcho(2:12); + xsuper=spcho("xsuper"); + nsuper=size(xsuper,1)-1; + snode=spcho("snode"); + neqns=size(snode,1); + lnz=spcho("lnz"); + nnzl=size(lnz,1); + iwsiz = 2*neqns+2*nsuper; + iwork = zeros(iwsiz,1); + tmpsiz = spcho("tmpsiz"); + tmpvec= zeros(tmpsiz,1); + iflag=0; + // calling blkfc1i primitive + [lnz,iflag]=blkfc1i(neqns,nsuper,xsuper,snode,spcho("split"),spcho("xlindx"),spcho("lindx"),spcho("xlnz"),lnz,iwsiz,iwork,tmpsiz,tmpvec,iflag,level); + // + if max(abs(lnz)) > 5d63 then + warning(gettext(" Possible matrix is not positive definite.")); + end; + + select iflag + case -1 then, + error(gettext("Non-positive diagonal encountered, matrix is not positive definite.")), + case -2 then, + error(gettext("Insufficient working storage in blkfct, temp(*).")), + case -3 then, + error(gettext("Insufficient working storage in blkfct, iwork(*).")), + end; + // + spcho("lnz")=lnz; +endfunction + +function rhs=blkslv(spcho,rhs) + // + //[xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp,lnz]=spcho(2:12); + xsuper=spcho("xsuper"); + nsuper=size(xsuper,1)-1; + neqns =size(rhs,1); + rhs=blkslvi(nsuper,xsuper,spcho("xlindx"),spcho("lindx"),spcho("xlnz"),... + spcho("lnz"),rhs); +endfunction + +function [spcho]=inpnv(xadjf,adjf,anzf,spcho) + // + //[xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp,lnz]=spcho(2:12); + // + xsuper=spcho("xsuper"); + neqns=size(xadjf,1)-1, + nsuper=size(xsuper,1)-1, + // + offset=zeros(neqns,1); + lnz=zeros(spcho("nnzl"),1); + // + lnz = inpnvi(neqns,xadjf,adjf,anzf,spcho("perm"),spcho("invp"),nsuper,... + xsuper,spcho("xlindx"),spcho("lindx"),spcho("xlnz"),lnz,offset); + // + spcho("lnz")=lnz; +endfunction + +function [spcho] = symfct(xadj,adjncy,perm,invp,cachsz,neqns) + // + // sfinit - input + // + nnza=size(adjncy,1); + iwsiz = 7*neqns+4; + iwork=zeros(iwsiz,1); + /// + if size(perm)~= [neqns,1] then, error(gettext(" SYMFCT requires PERM to be neqns x 1")), + end; + if size(invp)~= [neqns,1] then, error(gettext(" SYMFCT requires INVN to be neqns x 1")), + end; + if size(cachsz)~= [1,1] then, error(gettext(" SYMFCT requires CACHSZ to be 1 x 1")), + end; + // + [perm,invp,colcnt,nnzl,nsub,nsuper,snode,xsuper,iflag]=... + sfinit(neqns,nnza,xadj,adjncy,perm,invp,iwsiz,iwork); + // + if iflag == -1 then error(gettext(" Insufficient working storage in sfinit")),end; + // + bb=xsuper(1:nsuper+1,1); + xsuper=bb + iwsiz = 2*nsuper+2*neqns+1; + iwork=zeros(iwsiz,1); + // + // + [xlindx,lindx,xlnz,iflag]=symfcti(neqns,nnza,xadj,adjncy,perm,... + invp,colcnt,nsuper,xsuper,snode,nsub,iwsiz,iwork) + if iflag == -1 then error(" Insufficient working storage in symfct"),end; + if iflag == -2 then error(" Inconsistancy in the input in symfct"),end; + // + [tmpsiz,split]=bfinit(neqns,nsuper,xsuper,snode,xlindx,lindx,cachsz) + + spcho = tlist(["chol","xlnz","nnzl","xsuper","xlindx","lindx","snode","split",... + "tmpsiz","perm","invp","lnz"],... + xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp,[]) +endfunction diff --git a/modules/sparse/macros/chsolve.bin b/modules/sparse/macros/chsolve.bin Binary files differnew file mode 100755 index 000000000..a48c9d05c --- /dev/null +++ b/modules/sparse/macros/chsolve.bin diff --git a/modules/sparse/macros/chsolve.sci b/modules/sparse/macros/chsolve.sci new file mode 100755 index 000000000..21ecbd7a1 --- /dev/null +++ b/modules/sparse/macros/chsolve.sci @@ -0,0 +1,22 @@ + +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) ????-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 + +function sol=chsolve(spcho,rhs) + // Cholesky solver for A*sol=rhs (A is symmetric >0) + // 1rst step: spcho=chfact(A) + // 2nd step: sol=chsolve(spcho,rhs) + // Example: N=20; A=sprand(N,N,0.1); + // A=A*A'+speye(A); + // sol=(1:N)'; rhs=A*sol; + // spcho=chfact(A); sol=chsolve(spcho,rhs) + perm=spcho("perm"); + sol=blkslv(spcho,rhs(perm)); + sol(perm)=sol +endfunction diff --git a/modules/sparse/macros/cleanmacros.bat b/modules/sparse/macros/cleanmacros.bat new file mode 100755 index 000000000..768190a8b --- /dev/null +++ b/modules/sparse/macros/cleanmacros.bat @@ -0,0 +1,13 @@ +rem Scilab ( http://mwww.scilab.org/ ) - This file is part of Scilab +rem Copyright (C) INRIA +rem +rem This file must be used under the terms of the CeCILL. +rem This source file is licensed as described in the file COPYING, which +rem you should have received as part of this distribution. The terms +rem are also available at +rem http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt + + +@del *.bin 2>NUL +@del lib 2>NUL +@del names 2>NUL
\ No newline at end of file diff --git a/modules/sparse/macros/conjgrad.bin b/modules/sparse/macros/conjgrad.bin Binary files differnew file mode 100755 index 000000000..6c1566c7e --- /dev/null +++ b/modules/sparse/macros/conjgrad.bin diff --git a/modules/sparse/macros/conjgrad.sci b/modules/sparse/macros/conjgrad.sci new file mode 100755 index 000000000..07a7da992 --- /dev/null +++ b/modules/sparse/macros/conjgrad.sci @@ -0,0 +1,280 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: transformed into a gateway to +// propose pcg, cgs, bicg and bicgstab. +// Copyright (C) 2009 - INRIA - Michael Baudin +// Copyright (C) 2008 - INRIA - Michael Baudin +// Copyright (C) 2006 - INRIA - Serge Steer +// Copyright (C) 2005 - IRISA - Sage Group +// +// 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 + +// +// conjgrad -- +// This function regroups four methods from the "Conjugate Gradient family" to solve the linear system %Ax=b: +// - PCG (Preconditioned Conjugate Gradient): A must be symmetric positive definite, +// - CGS (preconditioned Conjugate Gradient Squared): A must be square, +// - BICG (preconditioned BiConjugate Gradient): A must be square, +// - BICGSTAB (preconditioned BiConjugate Gradient Stabilized): A must be square (default 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 +// flag INTEGER: 0 = solution found to tolerance +// 1 = no convergence given maxIter +// resNorm REAL final relative norm of the residual +// iter INTEGER number of iterations performed +// resVec REAL residual vector +// +// References +// +// PCG +// "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). +// +// Golub and Van Loan, Matrix Computations +// +// CGS +// "CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems" +// by Peter Sonneveld +// +// http://epubs.siam.org/doi/abs/10.1137/0910004 +// http://dl.acm.org/citation.cfm?id=64888&preflayout=flat +// http://mathworld.wolfram.com/ConjugateGradientSquaredMethod.html +// +// BICG +// "Numerical Recipes: The Art of Scientific Computing." (third ed.) +// by William Press, Saul Teukolsky, William Vetterling, Brian Flannery. +// +// http://apps.nrbook.com/empanel/index.html?pg=87 +// http://dl.acm.org/citation.cfm?doid=1874391.187410 +// http://mathworld.wolfram.com/BiconjugateGradientMethod.html +// +// BICGSTAB +// "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 +// +// The input / output arguments of this command are the same as +// Matlab's pcg, cgs, bicg and bicgstab commands, augmented with the 'method' argument +// + +function [x, flag, resNorm, iter, resVec] = conjgrad(%A, %b, method, tol, maxIter, %M, %M2, x0, verbose ) + + [lhs, rhs] = argn(0); + + if rhs < 2 then + error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"conjgrad",2,7)); + end + if rhs > 8 then + error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"conjgrad",2,7)); + end + if exists("method", "local") == 0 then + method = "bicgstab"; + end + if exists("tol", "local") == 0 then + tol = 1e-8 + end + if exists("maxIter", "local") == 0 then + maxIter = size(%b, 1) + end + if exists("%M", "local") == 0 then + %M = [] + end + if exists("%M2", "local") == 0 then + %M2 = [] + end + if exists("x0", "local") == 0 then + x0 = zeros(%b); + end + if exists("verbose", "local") == 0 then + verbose = 0; + end + if verbose == 1 then + printf(gettext("Arguments:\n")); + printf(" tol = "+string(tol)+"\n"); + printf(" maxIter = "+string(maxIter)+"\n"); + printf(" M = \n") + disp(%M) + printf(" M2 = \n"); + disp(%M2) + printf(" x0 = \n"); + disp(x0) + printf(" verbose = "+string(verbose)+"\n"); + end + // Compute matrixType + select type(%A) + case 1 then + matrixType = 1; + case 5 then + matrixType = 1; + case 13 then + matrixType = 0; + Aargs = list() + case 15 then + Aargs = list(%A(2:$)) + // Caution : modify the input argument %A ! + %A = %A(1); + matrixType = 0; + else + error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",1)); + end + // If %A is a matrix (dense or sparse) + if matrixType == 1 then + if size(%A,1) ~= size(%A,2) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",1)); + end + end + // Check right hand side %b + if type(%b) ~= 1 + error(msprintf(gettext("%s: Wrong type for input argument #%d: A matrix expected.\n"),"conjgrad",2)); + end + if size(%b,2) ~= 1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Column vector expected.\n"),"conjgrad",2)); + end + if matrixType == 1 then + if size(%b,1) ~= size(%A,1) then + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1)); + end + end + if matrixType == 1 then + if size(%b,1) ~= size(%A,1) then + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1)); + end + end + // Check method + if type(method) ~= 10 | size(method) ~= [1 1] + error(msprintf(gettext("%s: Wrong type for input argument #%d: Single String expected.\n"),"conjgrad",3)); + elseif and(method ~= ["pcg" "cgs" "bicg" "bicgstab"]), + error(msprintf(gettext("%s: Wrong value for input argument #%d: %s, %s, %s or %s expected.\n"),"conjgrad",3,"pcg","cgs","bicg","bicgstab")); + end + // Check type of the error tolerance tol + if or(size(tol) ~= [1 1]) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Scalar expected.\n"),"conjgrad",4)); + end + // Check the type of maximum number of iterations maxIter + if or(size(maxIter) ~= [1 1]) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Scalar expected.\n"),"conjgrad",5)); + end + // Compute precondType + select type(%M) + case 1 then + // M is a matrix + // precondType = 0 if the M is empty + // = 1 if the M is not empty + precondType = bool2s(size(%M,"*")>=1); + case 5 then + precondType = 1; + case 13 then + Margs = list() + precondType = 2; + case 15 then + Margs = list(%M(2:$)) + // Caution : modify the input argument %M ! + %M = %M(1); + precondType = 2; + else + error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",6)); + end + if precondType == 1 then + if size(%M,1) ~= size(%M,2) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",6)); + end + if size(%M,1) ~= size(%b,1) then + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",6,2)); + end + end + // Compute precondBis + select type(%M2) + case 1 then + // M2 is a matrix + // precondBis = 0 if the M2 is empty + // = 1 if the M2 is not empty + precondBis = bool2s(size(%M2,"*")>=1); + case 5 then + precondBis = 1; + case 13 then + M2args = list() + precondBis = 2; + case 15 then + M2args = list(%M2(2:$)) + // Caution : modify the input argument %M2 ! + %M2 = %M2(1); + // Caution : modify precondType again ! + precondType = 2; + else + error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",7)); + end + if precondBis == 1 then + if size(%M2,1) ~= size(%M2,2) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",7)); + end + if size(%M2,1) ~= size(%b,1) then + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",7,2)); + end + end + // Check size of the initial vector x0 + if size(x0,2) ~= 1 then + error(msprintf(gettext("%s: Wrong value for input argument #%d: Column vector expected.\n"),"conjgrad",8)); + end + if size(x0,1) ~= size(%b,1) then + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",8,2)); + end + + // ------------ + // Computations + // ------------ + select method + case "pcg" + [x, resNorm, iter, resVec] = %pcg(%A, %b, tol, maxIter, %M, %M2, x0, verbose ) + case "cgs" + [x, resNorm, iter, resVec] = %cgs(%A, %b, tol, maxIter, %M, %M2, x0, verbose ) + case "bicg" + [x, resNorm, iter, resVec] = %bicg(%A, %b, tol, maxIter, %M, %M2, x0, verbose ) + else // "bicgstab" + [x, resNorm, iter, resVec] = %bicgstab(%A, %b, tol, maxIter, %M, %M2, x0, verbose ) + 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"),"conjgrad")); + end + else + flag = 0; + if verbose == 1 then + printf(gettext("Algorithm pass\n")); + end + end + +endfunction diff --git a/modules/sparse/macros/gmres.bin b/modules/sparse/macros/gmres.bin Binary files differnew file mode 100755 index 000000000..af336c31c --- /dev/null +++ b/modules/sparse/macros/gmres.bin 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 diff --git a/modules/sparse/macros/issparse.bin b/modules/sparse/macros/issparse.bin Binary files differnew file mode 100755 index 000000000..cf8e972d7 --- /dev/null +++ b/modules/sparse/macros/issparse.bin diff --git a/modules/sparse/macros/issparse.sci b/modules/sparse/macros/issparse.sci new file mode 100755 index 000000000..c0a109ac8 --- /dev/null +++ b/modules/sparse/macros/issparse.sci @@ -0,0 +1,16 @@ +// 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 + + +function r=issparse(x) + r=0; + if or(type(x)==[5,7]) then + r=1 + end +endfunction diff --git a/modules/sparse/macros/lib b/modules/sparse/macros/lib Binary files differnew file mode 100755 index 000000000..59c004d55 --- /dev/null +++ b/modules/sparse/macros/lib diff --git a/modules/sparse/macros/names b/modules/sparse/macros/names new file mode 100755 index 000000000..edde6e8fd --- /dev/null +++ b/modules/sparse/macros/names @@ -0,0 +1,15 @@ +%bicg +%bicgstab +%cgs +%pcg +adj2sp +chfact +chsolve +conjgrad +gmres +issparse +qmr +sp2adj +speye +sprand +spzeros diff --git a/modules/sparse/macros/qmr.bin b/modules/sparse/macros/qmr.bin Binary files differnew file mode 100755 index 000000000..daf5a12ac --- /dev/null +++ b/modules/sparse/macros/qmr.bin diff --git a/modules/sparse/macros/qmr.sci b/modules/sparse/macros/qmr.sci new file mode 100755 index 000000000..739058377 --- /dev/null +++ b/modules/sparse/macros/qmr.sci @@ -0,0 +1,491 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) XXXX-2008 - INRIA +// Copyright (C) 2005 - IRISA - Sage Group + +// +// 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, err, iter, res] = qmr( A, Ap, b, x, M1, M1p, M2, M2p, max_it, tol ) +// +// QMR solves the linear system Ax=b using the +// Quasi Minimal Residual method with preconditioning. +// +// input A REAL matrix or function +// x REAL initial guess vector +// b REAL right hand side vector +// M1 REAL left preconditioner matrix +// M2 REAL right preconditioner matrix +// 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 +// breakdown: +// -1: rho +// -2: Beta +// -3: gam +// -4: delta +// -5: ep +// -6: xi +// err REAL final residual norm +// iter INTEGER number of iterations performed +// res 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). + +// Sage Group (IRISA, 2005) + +function [x, flag, err, iter, res] = qmr( A, varargin) + + // ----------------------- + // Parsing input arguments + // ----------------------- + [lhs,rhs]=argn(0); + if ( rhs < 2 ), + error(msprintf(gettext("%s: Wrong number of input arguments: At least %d expected.\n"),"qmr",2)); + end + + // Parsing the matrix A + select type(A) + case 1 then + cpt=1; + case 5 then + cpt=1; + case 13 then + cpt=0; + else + error(msprintf(gettext("%s: Wrong type for input argument #%d : A real or complex matrix or a sparse matrix or a function expected.\n"),"qmr",1)); + end + + // If A is a matrix (dense or sparse) + if (cpt==1), + if (size(A,1) ~= size(A,2)), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Square matrix expected.\n"),"qmr",1)); + end + fct=0; + cptmat = 1; + end + + // If A is a function + if (cpt==0), + fct = 0; + if rhs >= 2, + funcorvec=varargin(1); + if and(type(funcorvec) <> [1 5 13]) then + error(msprintf(gettext("%s: Wrong value for input argument #%d: Transpose of the function %s expected.\n"),"qmr",2,"A")); + // if the following input argument is a sparse or dense matrix + elseif or(type(funcorvec) == [1 5]) then + if size(getfield(1,macrovar(A)),"*") == 1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: A transpose of the function %s expected.\n"), "qmr", 2, "A")); + end + matvec = A; + cptmat = 1; + // if the following input argument is a function + else + matvec = A; + matvecp = funcorvec; + cptmat = 2; + fct = 1; + warning(msprintf(gettext("%s : Calling qmr(A,Ap) is deprecated. Please see qmr documentation for more details.\n"),"qmr")); + end + end + end + + // Parsing right hand side b + if ( rhs >= fct+2 ), + b=varargin(fct+1); + // if b is not constant or sparse + if and(type(b) <> [1 5]) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: A real or complex, full or sparse column vector expected.\n"), "qmr", fct+2)); + end + if ( size(b,2) ~= 1), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Column vector expected.\n"),"qmr",fct+2)); + end + end + + // Parsing initial vector x + if ( rhs >= fct+3), + x=varargin(fct+2); + // if x is not constant or sparse + if and(type(x) <> [1 5]) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: A real or complex, full or sparse column vector expected.\n"), "qmr", fct + 3)); + end + if (size(x,2) ~= 1), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Column vector expected.\n"),"qmr",fct+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"),"qmr",fct+3,fct+2)); + end + // By default + else + x=zeros(size(b,1),1); + end + + //-------------------------------------------------------- + // Parsing of the preconditioner matrix M1 + //-------------------------------------------------------- + if (rhs >=fct+4), + Prec_g=varargin(fct+3); + select type(Prec_g) + case 1 then + cpt=1; + case 5 then + cpt=1; + case 13 then + cpt=0; + else + error(msprintf(gettext("%s: Wrong type for input argument #%d: A real or complex, full or sparse, square matrix or a function expected.\n"), "qmr", fct + 4)); + end + + // if M1 is a matrix + if ( cpt==1 ), + if (size(Prec_g,1) ~= size(Prec_g,2)), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Square matrix expected.\n"),"qmr",fct+4)); + end + if (size(Prec_g,1)~=size(b,1)), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"qmr",fct+4,fct+2)); + end + cptmatM1 = 1; + end + + // if M1 is a function + if ( cpt==0 ), + precond_g = Prec_g; + cptmatM1 = 1; + if ( rhs >= fct+5 & size(getfield(1,macrovar(precond_g)),"*") == 1), + Precp_g = varargin(fct+4); + if (type(Precp_g) == 13 & size(getfield(1,macrovar(Precp_g)),"*")==1) then + //precond_g = Prec_g; + precondp_g = Precp_g; + cptmatM1 = 2; + fct = fct+1; + warning(msprintf(gettext("%s : Calling qmr(...,M1,M1p) is deprecated. Please see qmr documentation for more details.\n"),"qmr")); + else + error(msprintf(gettext("%s: Wrong type for input argument #%d: A transpose of the function expected.\n"), "qmr", fct + 5)); + end + elseif rhs < fct+5 & size(getfield(1, macrovar(precond_g)), "*") == 1 then + error(msprintf(gettext("%s: Wrong prototype of input argument #%d: If M1 is function, use the header M1(x,t) instead M1(x).\n"), "qmr", fct+4)); + end + end + // By default + else + deff("y=precond_g(x)","y=x"); + deff("y=precondp_g(x)","y=x"); + cptmatM1 = 2; + end + + //-------------------------------------------------------- + // Parsing of the preconditioner matrix M2 + //-------------------------------------------------------- + if (rhs >=fct+5), + Prec_d=varargin(fct+4); + select type(Prec_d) + case 1 then + cpt=1; + case 5 then + cpt=1; + case 13 then + cpt=0; + else + error(msprintf(gettext("%s: Wrong type for input argument #%d: A real or complex, full or sparse, square matrix or a function expected.\n"), "qmr", fct + 5)); + end + + // M2 matrix + if ( cpt==1 ), + + if (size(Prec_d,1) ~= size(Prec_d,2)), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Square matrix expected.\n"),"qmr",fct+5)); + end + if (size(Prec_d,1)~=size(b,1)), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"qmr",fct+5,fct+2)); + end + cptmatM2 = 1; + + end + + // M2 function + if ( cpt==0 ) + precond_d = Prec_d; + cptmatM2 = 1; + if ( rhs >= fct+6 & size(getfield(1,macrovar(precond_d)),"*") == 1), + Precp_d=varargin(fct+5); + if (type(Precp_d) == 13 & size(getfield(1,macrovar(Precp_d)),"*") == 1) then + precond_d = Prec_d; + precondp_d = Precp_d; + cptmatM2 = 2; + fct = fct+1; + warning(msprintf(gettext("%s : Calling qmr(...,M2,M2p) is deprecated. Please see qmr documentation for more details.\n"),"qmr")); + else + error(msprintf(gettext("%s: Wrong type for input argument #%d: A transpose of the function expected.\n"), "qmr", fct + 6)); + end + elseif rhs < fct+6 & size(getfield(1, macrovar(precond_d)), "*") == 1 then + error(msprintf(gettext("%s: Wrong prototype of input argument #%d: If M2 is function, use the header M2(x,t) instead M2(x).\n"), "qmr", fct+5)); + end + end + // By default + else + deff("y=precond_d(x)","y=x"); + deff("y=precondp_d(x)","y=x"); + cptmatM2 = 2; + end + + //-------------------------------------------------------- + // Parsing of the maximum number of iterations max_it + //-------------------------------------------------------- + if (rhs >= fct+6), + max_it=varargin(fct+5); + // if max_it is not constant + if type(max_it) <> 1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Scalar expected.\n"),"qmr",fct + 6)); + end + + if ~isreal(max_it) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: A real scalar expected.\n"),"qmr",fct + 6)); + end + + if (size(max_it,1) ~= 1 | size(max_it,2) ~=1), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"),"qmr",fct+6)); + end + // By default + else + max_it=size(b,1); + end + + //-------------------------------------------------------- + // Parsing of the error tolerance tol + //-------------------------------------------------------- + if (rhs == fct+7), + tol=varargin(fct+6); + // if tol is not constant + if type(tol) <> 1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Scalar expected.\n"),"qmr",fct + 7)); + end + + if ~isreal(tol) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: A real scalar expected.\n"),"qmr",fct + 7)); + end + + if (size(tol,1) ~= 1 | size(tol,2) ~=1), + error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"),"qmr",fct+7)); + end + // By default + else + tol=1000*%eps; + end + + //-------------------------------------------------------- + // test about input arguments number + //-------------------------------------------------------- + if (rhs > fct+7), + error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"qmr",2,fct+7)); + end + + // ------------ + // Computations + // ------------ + + // initialization + i = 0; + flag = 0; + iter = 0; + bnrm2 = norm( b ); + if (bnrm2 == 0.0), + bnrm2 = 1.0; + end + + // r = b - A*x; + if (cptmat == 1) then + r = b - matvec(x,"notransp"); + elseif (cptmat==2) then // If A is a function + r = b - matvec(x); + end + + err = norm( r ) / bnrm2; + res = err; + if ( err < tol ), return; end + + // [M1,M2] = lu( M ); + v_tld = r; + + // y = M1 \ v_tld; + if (cptmatM1 == 1) then + y = precond_g(v_tld,"notransp"); + elseif (cptmatM1==2) then + y = precond_g(v_tld); + end + + rho = norm( y ); + + w_tld = r; + // z = M2' \ w_tld; + if (cptmatM2 == 1) then + z = precond_d(w_tld,"transp"); + elseif (cptmatM2 == 2) then + z = precondp_d(w_tld); + end + + xi = norm( z ); + + gam = 1.0; + eta = -1.0; + theta = 0.0; + + for i = 1:max_it, // begin iteration + if ( rho == 0.0 | xi == 0.0 ), iter=i; break; end + v = v_tld / rho; + y = y / rho; + + w = w_tld / xi; + z = z / xi; + + delta = z'*y; + if ( delta == 0.0 ), iter=i; break; end + + // y_tld = M2 \ y; + if (cptmatM2 == 1) then + y_tld = precond_d(y,"notransp"); + elseif (cptmatM2 == 2) then + y_tld = precond_d(y); + end + + // z_tld = M1'\ z; + if (cptmatM1 == 1) then + z_tld = precond_g(z,"transp"); + elseif (cptmatM1 == 2) then + z_tld = precondp_g(z); + end + + if ( i > 1 ), // direction vector + p = y_tld - ( xi*delta / ep )*p; + q = z_tld - ( rho*delta / ep )*q; + else + p = y_tld; + q = z_tld; + end + + // p_tld = A*p; + if (cptmat == 1) then + p_tld = matvec(p,"notransp"); + elseif (cptmat == 2) then + p_tld = matvec(p); + end + + ep = q'*p_tld; + if ( ep == 0.0 ), iter=i; break; end + + Beta = ep / delta; + if ( Beta == 0.0 ), iter=i; break; end + + v_tld = p_tld - Beta*v; + + // y = M1 \ v_tld; + if (cptmatM1 == 1) then + y = precond_g(v_tld,"notransp"); + elseif (cptmatM1==2) then + y = precond_g(v_tld); + end + + rho_1 = rho; + rho = norm( y ); + + // w_tld = ( A'*q ) - ( Beta*w ); + if (cptmat == 1) then + w_tld = ( matvec(q,"transp") ) - ( Beta*w ); + elseif (cptmat == 2) then + w_tld = ( matvecp(q) ) - ( Beta*w ); + end + + // z = M2' \ w_tld; + if (cptmatM2 == 1) then + z = precond_d(w_tld,"transp"); + elseif (cptmatM2 == 2) then + z = precondp_d(w_tld); + end + + xi = norm( z ); + gamma_1 = gam; + theta_1 = theta; + theta = rho / ( gamma_1*Beta ); + gam = 1.0 / sqrt( 1.0 + (theta^2) ); + if ( gam == 0.0 ), iter=i; break; end + + eta = -eta*rho_1*(gam^2) / ( Beta*(gamma_1^2) ); + + if ( i > 1 ), // compute adjustment + d = eta*p + (( theta_1*gam )^2)*d; + s = eta*p_tld + (( theta_1*gam )^2)*s; + else + d = eta*p; + s = eta*p_tld; + end + x = x + d; // update approximation + + r = r - s; // update residual + err = norm( r ) / bnrm2; // check convergence + res = [res;err]; + + if ( err <= tol ), iter=i; break; end + + if ( i == max_it ), iter=i; end + end + + if ( err <= tol ), // converged + flag = 0; + elseif ( rho == 0.0 ), // breakdown + flag = -1; + elseif ( Beta == 0.0 ), + flag = -2; + elseif ( gam == 0.0 ), + flag = -3; + elseif ( delta == 0.0 ), + flag = -4; + elseif ( ep == 0.0 ), + flag = -5; + elseif ( xi == 0.0 ), + flag = -6; + else // no convergence + flag = 1; + end + +endfunction + +function y = matvec(x,t) + if (t=="notransp") then + y = A*x; + elseif (t=="transp") then + y = A'*x; + end +endfunction + +function y = precond_g(x,t) + if (t=="notransp") then + y = Prec_g*x; + elseif (t=="transp") then + y = Prec_g'*x; + end +endfunction + +function y = precond_d(x,t) + if (t=="notransp") then + y = Prec_d*x; + elseif (t=="transp") then + y = Prec_d'*x; + end +endfunction + diff --git a/modules/sparse/macros/sp2adj.bin b/modules/sparse/macros/sp2adj.bin Binary files differnew file mode 100755 index 000000000..b44ecd4d7 --- /dev/null +++ b/modules/sparse/macros/sp2adj.bin diff --git a/modules/sparse/macros/sp2adj.sci b/modules/sparse/macros/sp2adj.sci new file mode 100755 index 000000000..da86ffdba --- /dev/null +++ b/modules/sparse/macros/sp2adj.sci @@ -0,0 +1,50 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) XXXX-2008 - INRIA +// Copyright (C) 2010 - DIGITEO - Michael Baudin +// +// 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 + + +function [xadj,iadj,v]=sp2adj(A) + // sparse conversion to adjacency + // See sp2adj.xml for help. + // + // Check number of arguments + [lhs, rhs] = argn() + if ( rhs <> 1 ) then + lstr = gettext("%s: Wrong number of input arguments: %d to %d expected.\n") + error ( msprintf(lstr,"sp2adj",1,1)) + end + if ( and(lhs <> [0 1 2 3]) ) then + lstr = gettext("%s: Wrong number of output arguments: %d to %d expected.\n") + error ( msprintf(lstr,"sp2adj",0,3)) + end + // + // Check type of arguments + if ( typeof(A) <> "sparse" ) then + lstr = gettext("%s: Wrong type for input argument #%d.\n") + error ( msprintf(lstr,"sp2adj",1)) + end + // + // Check size of arguments + // Nothing to do + // + // Check content of arguments + // Nothing to do + // + [ij,v,n]=spget(A') + N=n(1) + if ( ij == [] ) then, + xadj=ones(n(2)+1,1) + iadj=[] + v=[] + else, + [xadj,la,iadj]=ta2lpd(ij(:,1)',ij(:,2)',N+1,N) + xadj=xadj(:) + iadj=iadj(:) + end +endfunction diff --git a/modules/sparse/macros/speye.bin b/modules/sparse/macros/speye.bin Binary files differnew file mode 100755 index 000000000..0732eca35 --- /dev/null +++ b/modules/sparse/macros/speye.bin diff --git a/modules/sparse/macros/speye.sci b/modules/sparse/macros/speye.sci new file mode 100755 index 000000000..56d62ebd9 --- /dev/null +++ b/modules/sparse/macros/speye.sci @@ -0,0 +1,16 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 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 + +function [sp]=speye(m,n) + + [lhs,rhs]=argn(0) + if rhs==1 then [m,n]=size(m),end + mn=min(m,n) + sp=sparse([(1:mn)' (1:mn)'],ones(mn,1),[m,n]) +endfunction diff --git a/modules/sparse/macros/sprand.bin b/modules/sparse/macros/sprand.bin Binary files differnew file mode 100755 index 000000000..17f3057d3 --- /dev/null +++ b/modules/sparse/macros/sprand.bin diff --git a/modules/sparse/macros/sprand.sci b/modules/sparse/macros/sprand.sci new file mode 100755 index 000000000..e42c3e5bb --- /dev/null +++ b/modules/sparse/macros/sprand.sci @@ -0,0 +1,74 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) INRIA +// Copyright (C) 2010 - DIGITEO - Michael Baudin +// +// 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 + +function a=sprand(m,n,density,typ) + + //---- Check arguments---------------------------------------------- + rhs = argn(2) + + if ( rhs < 3 | rhs > 4 ) then + error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"), "sprand" , 3 , 4 )); + end + + if ( rhs < 4 ) then + typ="def"; + end + + if type(typ)<>10 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: A string expected.\n"),"sprand",4)); + end + + if and(typ<>["u";"n";"uniform";"normal";"def";"nor"]) then + error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'' or ''%s'' expected.\n"),"sprand",4,"uniform","normal")); + end + if typ == "u" | typ == "uniform" then //"uniform" is the syntax for uniform distribution with rand, the equivalent with grand is "def" + typ = "def"; + elseif typ == "n" | typ == "normal" then //"normal" is the syntax for normal distribution with rand, the equivalent with grand is "nor" + typ = "nor"; + end + + density=max(min(density,1),0); + + nel=m*n*density; //the objective number of non zero elements + if nel==0 then + a=sparse([],[],[m,n]) + return + end + + //---- generate a sequence of increments---------------------------- + mdist = 1/density //the mean distance between to consecutive index + nel1 = (2.2-density)*nel; //generate more increments than requested nnz elements + ij = round(1+grand(nel1,1,"exp",(mdist-1))) + + //---- sum the increments to get the index-------------------------- + ij=cumsum(ij); + + //---- eliminate the index with exceed the maximum matrix index + ij(find(ij>m*n))=[]; + nel1=size(ij,"*"); + if nel1==0 then + a=sparse([],[],[m,n]); + return + end + //---- generate the row and column indices from the 1D index-------- + + ij=ind2sub([m,n],ij); + + //---- generates the random non zeros elements -------------------- + //according to the requested law and create the sparse matrix + if typ == "nor" then // Because of the syntax of grand, we have two cases, one with "def" and one with "nor" + a=sparse(ij,grand(nel1,1,typ,0,1),[m,n]); + elseif typ == "def" then + a=sparse(ij,grand(nel1,1,typ),[m,n]); + end + + +endfunction + diff --git a/modules/sparse/macros/spzeros.bin b/modules/sparse/macros/spzeros.bin Binary files differnew file mode 100755 index 000000000..d8d462971 --- /dev/null +++ b/modules/sparse/macros/spzeros.bin diff --git a/modules/sparse/macros/spzeros.sci b/modules/sparse/macros/spzeros.sci new file mode 100755 index 000000000..a6373f4a5 --- /dev/null +++ b/modules/sparse/macros/spzeros.sci @@ -0,0 +1,21 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 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 + +function [sp]=spzeros(m,n) + + [lhs,rhs]=argn(0) + if ( rhs < 1 ) then + error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"), "spzeros" , 1 , 2 )); + end + + if rhs==1 then + [m,n]=size(m) + end + sp=sparse([],[],[m,n]) +endfunction |