* Bugs #11007 11008 11009 fixed - Sparse: cgs, bicg, bicgstab 07/13007/7
Paul BIGNIER [Thu, 24 Oct 2013 15:31:54 +0000 (17:31 +0200)]
New function conjgrad to regroup the Conjugate Gradient solvers pcg, cgs, bicg and bicgstab.

Change-Id: Iebd3943a92f7c3a61b8e53e6c7415338ed52eb99

33 files changed:
SEP/INDEX
SEP/SEP_116_conjgrad.odt [new file with mode: 0644]
scilab/CHANGES_5.5.X
scilab/modules/sparse/help/en_US/iterativesolvers/conjgrad.xml [new file with mode: 0644]
scilab/modules/sparse/help/en_US/iterativesolvers/gmres.xml
scilab/modules/sparse/help/en_US/iterativesolvers/pcg.xml
scilab/modules/sparse/help/en_US/iterativesolvers/qmr.xml
scilab/modules/sparse/macros/%bicg.sci [new file with mode: 0644]
scilab/modules/sparse/macros/%bicgstab.sci [new file with mode: 0644]
scilab/modules/sparse/macros/%cgs.sci [new file with mode: 0644]
scilab/modules/sparse/macros/%pcg.sci [new file with mode: 0644]
scilab/modules/sparse/macros/conjgrad.sci [new file with mode: 0644]
scilab/modules/sparse/macros/pcg.sci
scilab/modules/sparse/tests/nonreg_tests/bug_3310.dia.ref
scilab/modules/sparse/tests/unit_tests/conjgrad.dia.ref [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/conjgrad.tst [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/conjgrad_function.dia.ref [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/conjgrad_function.tst [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/conjgrad_list.dia.ref [moved from scilab/modules/sparse/tests/unit_tests/pcg_list.dia.ref with 59% similarity]
scilab/modules/sparse/tests/unit_tests/conjgrad_list.tst [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/conjgrad_numerical.dia.ref [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/conjgrad_numerical.tst [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/conjgrad_sparse.dia.ref [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/conjgrad_sparse.tst [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/pcg.dia.ref [deleted file]
scilab/modules/sparse/tests/unit_tests/pcg.tst [deleted file]
scilab/modules/sparse/tests/unit_tests/pcg_function.dia.ref [deleted file]
scilab/modules/sparse/tests/unit_tests/pcg_function.tst [deleted file]
scilab/modules/sparse/tests/unit_tests/pcg_list.tst [deleted file]
scilab/modules/sparse/tests/unit_tests/pcg_numerical.dia.ref [deleted file]
scilab/modules/sparse/tests/unit_tests/pcg_numerical.tst [deleted file]
scilab/modules/sparse/tests/unit_tests/pcg_sparse.dia.ref [deleted file]
scilab/modules/sparse/tests/unit_tests/pcg_sparse.tst [deleted file]

index 44a85b4..cef2c09 100644 (file)
--- a/SEP/INDEX
+++ b/SEP/INDEX
@@ -111,3 +111,4 @@ SEP #112: variance and variancef can now return the mean of the input in a new o
 SEP #113: resize_matrix hypermatrices support.
 SEP #114: New function jcreatejar.
 SEP #115: New function bode_asymp.
+SEP #116: conjgrad function, to regroup the Conjugate Gradient solvers. pcg obsolete as it is.
diff --git a/SEP/SEP_116_conjgrad.odt b/SEP/SEP_116_conjgrad.odt
new file mode 100644 (file)
index 0000000..953a3a3
Binary files /dev/null and b/SEP/SEP_116_conjgrad.odt differ
index c8ad0d1..084e2c0 100644 (file)
@@ -19,6 +19,8 @@ Obsolete & Removed Functions
 * msd and st_deviation tagged as obsolete (See bug #7593). Will be removed in Scilab 5.5.1.
   Please use stdev instead.
 
+* pcg tagged as obsolete. Will be removed in Scilab 5.5.1.
+  Please use conjgrad instead.
 
 Scilab Bug Fixes
 ================
@@ -89,6 +91,8 @@ Scilab Bug Fixes
 
 * Bug #10273 fixed - spchol help page now displays an example showing how to use its output arguments.
 
+* Bug #11007, #11008 & #11009 fixed - New function conjgrad (Conjugate Gradient methods "pcg", "cgs", "bicg" and "bicgstab").
+
 * Bug #11571 fixed - x_mdialog did not let the Look&Feel select the window size.
 
 * Bug #11575 fixed - There was no preview of GIF files in exportUI dialog.
@@ -405,8 +409,6 @@ Obsolete & Removed Functions
 
 * milk_drop tagged as obsolete. Will be removed in Scilab 5.5.1.
 
-* st_deviation and msd tagged as obsolete, please use stdev intead. Will be removed in Scilab 5.5.1.
-
 
 Scilab Bug Fixes
 ================
diff --git a/scilab/modules/sparse/help/en_US/iterativesolvers/conjgrad.xml b/scilab/modules/sparse/help/en_US/iterativesolvers/conjgrad.xml
new file mode 100644 (file)
index 0000000..e011282
--- /dev/null
@@ -0,0 +1,406 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * 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
+ *
+ -->
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="conjgrad" xml:lang="en">
+    <refnamediv>
+        <refname>conjgrad</refname>
+        <refpurpose>conjugate gradient solvers</refpurpose>
+    </refnamediv>
+    <refsynopsisdiv>
+        <title>Calling Sequence</title>
+        <synopsis>
+            [x, flag, err, iter, res] = conjgrad(A, b [, method [, tol [, maxIter [, M [, M2 [, x0 [, verbose]]]]]]])
+            [x, flag, err, iter, res] = conjgrad(A, b [, method [, key=value,...]])
+        </synopsis>
+    </refsynopsisdiv>
+    <refsection>
+        <title>Arguments</title>
+        <variablelist>
+            <varlistentry>
+                <term>A</term>
+                <term/>
+                <listitem>
+                    <para>a matrix, or a function, or a list computing
+                        <literal>A*x</literal> for each given <literal>x</literal>. The
+                        following is a description of the computation of A*x depending on
+                        the type of A.
+                    </para>
+                    <itemizedlist>
+                        <listitem>
+                            <para>
+                                <literal>matrix.</literal>If A is a matrix, it can be
+                                dense or sparse
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>
+                                <literal>function.</literal>If A is a function, it must
+                                have the following header :
+                            </para>
+                            <programlisting role=""><![CDATA[
+function y = A ( x )
+ ]]></programlisting>
+                        </listitem>
+                        <listitem>
+                            <para>
+                                <literal>list.</literal>If A is a list, the first element
+                                of the list is expected to be a function and the other elements
+                                in the list are the arguments of the function, from index 2 to
+                                the end. When the function is called, the current value of x is
+                                passed to the function as the first argument. The other
+                                arguments passed are the one given in the list.
+                            </para>
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>b</term>
+                <listitem>
+                    <para>right hand side vector (size: nx1)</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>mehtod</term>
+                <listitem>
+                    <para>scalar string, "pcg", "cgs", "bicg" or "bicgstab" (default "bicgstab")</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>tol</term>
+                <listitem>
+                    <para>error relative tolerance (default: 1e-8).
+                        The termination criteria is based on the 2-norm of the
+                        residual r=b-Ax, divided by the 2-norm of the right hand side b.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>maxIter</term>
+                <listitem>
+                    <para>maximum number of iterations (default: n)</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>M</term>
+                <listitem>
+                    <para>preconditioner: full or sparse matrix or function returning
+                        <literal>M\x</literal> (default: none)
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>M2</term>
+                <listitem>
+                    <para>preconditioner: full or sparse matrix or function returning
+                        <literal>M2\x</literal> for each <literal>x</literal> (default:
+                        none)
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>x0</term>
+                <listitem>
+                    <para>initial guess vector (default: zeros(n,1))</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>verbose</term>
+                <listitem>
+                    <para>set to 1 to enable verbose logging (default 0)</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>x</term>
+                <listitem>
+                    <para>solution vector</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>flag</term>
+                <listitem>
+                    <para>
+                        0 if <literal>conjgrad</literal> converged to the desired tolerance
+                        within <literal>maxi</literal> iterations, 1 else
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>err</term>
+                <listitem>
+                    <para>final relative norm of the residual (the 2-norm of the right-hand side b is used)</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>iter</term>
+                <listitem>
+                    <para>number of iterations performed</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>res</term>
+                <listitem>
+                    <para>vector of the residual relative norms</para>
+                </listitem>
+            </varlistentry>
+        </variablelist>
+    </refsection>
+    <refsection>
+        <title>Description</title>
+        <para>
+            Solves the linear system <literal>Ax=b</literal> using the conjugate
+            gradient method with or without preconditioning. The preconditionning
+            should be defined by a symmetric positive definite matrix
+            <literal>M</literal>, or two matrices <literal>M1</literal> and
+            <literal>M2</literal> such that <literal>M=M1*M2</literal>. in the case
+            the function solves <literal>inv(M)*A*x = inv(M)*b</literal> for
+            <literal>x</literal>. <literal>M</literal>, <literal>M1</literal> and
+            <literal>M2</literal> can be Scilab functions with calling sequence
+            <literal>y=Milx(x)</literal> which computes the corresponding left
+            division <literal>y=Mi\x</literal>.
+        </para>
+        <para>
+            The <literal>method</literal> input argument selects the solver to be used:
+            <simplelist type="inline">
+                <member>
+                    "pcg" Preconditioned Conjugate Gradient
+                </member>
+                <member>
+                    "cgs" preconditioned Conjugate Gradient Squared
+                </member>
+                <member>
+                    "bicg" preconditioned BiConjugate Gradient
+                </member>
+                <member>
+                    "bicgstab" preconditioned BiConjugate Gradient Stabilized (default)
+                </member>
+            </simplelist>
+        </para>
+        <para>
+            If <literal>method="pcg"</literal>, then the <literal>A</literal> matrix
+            must be a symmetric positive definite matrix (full or sparse)
+            or a function with calling sequence <literal>y=Ax(x)</literal> which computes
+            <literal>y=A*x</literal>.
+            Otherwise (<literal>method="cgs", "bicg" or "bicgstab"</literal>),
+            <literal>A</literal> just needs to be square.
+        </para>
+    </refsection>
+    <refsection>
+        <title>Example with well-conditionned and ill-conditionned
+            problems
+        </title>
+        <para>In the following example, two linear systems are solved. The first
+            maxtrix has a condition number equals to ~0.02, which makes the algorithm
+            converge in exactly 10 iterations. Since this is the size of the matrix,
+            it is an expected behaviour for a gradient conjugate method. The second
+            one has a low condition number equals to 1.d-6, which makes the algorithm
+            converge in a larger 22 iterations. This is why the parameter maxIter is
+            set to 30. See below for other examples of the "key=value" syntax.
+        </para>
+        <programlisting role="example"><![CDATA[
+// Well-conditionned problem
+A = [ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
+     0   0   65  55   0   0   0   0   11  100];
+
+b = ones(10, 1);
+[x, fail, err, iter, res] = conjgrad(A, b, "bicg", 1d-12, 15);
+mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)
+
+// Ill-contionned one
+A = [ 894     0   0     0   0   28  0   0   1000  70000
+      0      5   13    5   0   0   0   0   0     0
+      0      13  72    34  0   0   0   0   0     6500
+      0      5   34    1   0   0   0   0   0     55
+      0      0   0     0   70  0   28  32  12    0
+      28     0   0     0   0   87  20  0   33    0
+      0      0   0     0   28  20  71  39  0     0
+      0      0   0     0   32  0   39  46  8     0
+      1000   0   0     0   12  33  0   8   82    11
+      70000  0   6500  55  0   0   0   0   11    100];
+
+[x, fail, err, iter, res] = conjgrad(A, b, method="pcg", maxIter=30, tol=1d-12);
+mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)
+ ]]></programlisting>
+    </refsection>
+    <refsection>
+        <title>Examples with A given as a sparse matrix, or function, or
+            list
+        </title>
+        <para>The following example shows that the method can handle sparse
+            matrices as well. It also shows the case where a function, computing the
+            right-hand side, is given to the "conjgrad" primitive. The final case shown by
+            this example, is when a list is passed to the primitive.
+        </para>
+        <programlisting role="example"><![CDATA[
+// Well-conditionned problem
+A = [ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
+     0   0   65  55   0   0   0   0   11  100];
+b = ones(10, 1);
+
+// Convert A into a sparse matrix
+Asparse=sparse(A);
+[x, fail, err, iter, res] = conjgrad(Asparse, b, "cgs", maxIter=30, tol=1d-12);
+mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)
+
+// Define a function which computes the right-hand side.
+function y = Atimesx(x)
+  A = [ 94  0   0   0    0   28  0   0   32  0
+       0   59  13  5    0   0   0   10  0   0
+       0   13  72  34   2   0   0   0   0   65
+       0   5   34  114  0   0   0   0   0   55
+       0   0   2   0    70  0   28  32  12  0
+       28  0   0   0    0   87  20  0   33  0
+       0   0   0   0    28  20  71  39  0   0
+       0   10  0   0    32  0   39  46  8   0
+       32  0   0   0    12  33  0   8   82  11
+       0   0   65  55   0   0   0   0   11  100];
+  y = A*x
+endfunction
+
+// Pass the script Atimesx to the primitive
+[x, fail, err, iter, res] = conjgrad(Atimesx, b, maxIter=30, tol=1d-12);
+mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)
+
+// Define a function which computes the right-hand side.
+function y = Atimesxbis(x, A)
+  y = A*x
+endfunction
+
+// Pass a list to the primitive
+Alist = list(Atimesxbis, Asparse);
+[x, fail, err, iter, res] = conjgrad(Alist, b, maxIter=30, tol=1d-12);
+mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)
+ ]]></programlisting>
+    </refsection>
+    <refsection>
+        <title>Examples with key=value syntax</title>
+        <para>The following example shows how to pass arguments with the
+            "key=value" syntax. This allows to set non-positionnal arguments, that is,
+            to set arguments which are not depending on their order in the list of
+            arguments. The available keys are the names of the optional arguments,
+            that is : tol, maxIter, %M, %M2, x0, verbose. Notice that, in the
+            following example, the verbose option is given before the maxIter option.
+            Without the "key=value" syntax, the positionnal arguments would require
+            that maxIter come first and verbose after.
+        </para>
+        <programlisting role="example"><![CDATA[
+// Example of an argument passed with key=value syntax
+A = [100 1; 1 10];
+b = [101; 11];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, verbose=1);
+
+// With key=value syntax, the order does not matter
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, verbose=1, maxIter=0);
+ ]]></programlisting>
+    </refsection>
+    <refsection role="see also">
+        <title>See Also</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="backslash">backslash</link>
+            </member>
+            <member>
+                <link linkend="qmr">qmr</link>
+            </member>
+            <member>
+                <link linkend="gmres">gmres</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection>
+        <title>References</title>
+        <para>
+            <emphasis role="bold">PCG</emphasis>
+        </para>
+        <para>"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/linalg/templates.ps
+        </para>
+        <para>"Iterative Methods for Sparse Linear Systems, Second Edition", Saad,
+            SIAM Publications, 2003, ftp
+            ftp.cs.umn.edu/dept/users/saad/PS/all_ps.zip
+        </para>
+        <para>
+            <emphasis role="bold">CGS</emphasis>
+        </para>
+        <para>
+            "CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems" by Peter Sonneveld.
+        </para>
+        <para>
+            <ulink url="http://epubs.siam.org/doi/abs/10.1137/0910004">Original article</ulink>
+        </para>
+        <para>
+            <ulink url="http://dl.acm.org/citation.cfm?id=64888&amp;preflayout=flat">Article on ACM</ulink>
+        </para>
+        <para>
+            <ulink url="http://mathworld.wolfram.com/ConjugateGradientSquaredMethod.html">Some theory around CGS</ulink>
+        </para>
+        <para>
+            <emphasis role="bold">BICG</emphasis>
+        </para>
+        <para>
+            "Numerical Recipes: The Art of Scientific Computing." (third ed.) by William Press, Saul Teukolsky, William Vetterling, Brian Flannery.
+        </para>
+        <para>
+            <ulink url="http://apps.nrbook.com/empanel/index.html?pg=87">http://apps.nrbook.com/empanel/index.html?pg=87</ulink>
+        </para>
+        <para>
+            <ulink url="http://dl.acm.org/citation.cfm?doid=1874391.187410">Article on ACM</ulink>
+        </para>
+        <para>
+            <ulink url="http://mathworld.wolfram.com/BiconjugateGradientMethod.html">Some theory around BICG</ulink>
+        </para>
+        <para>
+            <emphasis role="bold">BICGSTAB</emphasis>
+        </para>
+        <para>
+            "Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems" by Henk van der Vorst.        339
+        </para>
+        <para>
+            <ulink url="http://epubs.siam.org/doi/abs/10.1137/0913035">Original article</ulink>
+        </para>
+        <para>
+            <ulink url="http://dl.acm.org/citation.cfm?id=131916.131930&amp;coll=DL&amp;dl=GUIDE&amp;CFID=372773884&amp;CFTOKEN=56630250">Article on ACM</ulink>
+        </para>
+        <para>
+            <ulink url="http://mathworld.wolfram.com/BiconjugateGradientStabilizedMethod.html">Some theory around BICG</ulink>
+        </para>
+    </refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.5.0</revnumber>
+                <revdescription>
+                    Introduction
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
index 46e0ffa..f54361e 100644 (file)
@@ -3,11 +3,11 @@
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) XXXX-2008 - INRIA
  * Copyright (C) 2012 - Scilab Enterprises - Adeline Carnis
- * 
+ *
  * 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    
+ * are also available at
  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
  *
  -->
@@ -29,7 +29,7 @@
                     <para>
                         n-by-n matrix or function returning <literal>A*x</literal>. If <literal>A</literal> is a function, it must have the following header :
                     </para>
-                    <programlisting role=""><![CDATA[ 
+                    <programlisting role=""><![CDATA[
 function y = A ( x )
  ]]></programlisting>
                 </listitem>
@@ -52,7 +52,7 @@ function y = A ( x )
                     <para>
                         preconditioner: matrix of size n-by-n or function returning <literal>M*x</literal> (In the first case, default: eye(n,n)). If M is a function, it must have the following header :
                     </para>
-                    <programlisting role=""><![CDATA[ 
+                    <programlisting role=""><![CDATA[
 function y = M ( x )
  ]]></programlisting>
                 </listitem>
@@ -154,17 +154,17 @@ function y = M ( x )
     </refsection>
     <refsection>
         <title>Examples</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
        // If A and M are matrices
-A=[ 94  0   0   0    0   28  0   0   32  0  
-     0   59  13  5    0   0   0   10  0   0  
-     0   13  72  34   2   0   0   0   0   65 
-     0   5   34  114  0   0   0   0   0   55 
-     0   0   2   0    70  0   28  32  12  0  
-     28  0   0   0    0   87  20  0   33  0  
-     0   0   0   0    28  20  71  39  0   0  
-     0   10  0   0    32  0   39  46  8   0  
-     32  0   0   0    12  33  0   8   82  11 
+A=[ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
      0   0   65  55   0   0   0   0   11  100];
 b=ones(10,1);
 [x,flag,err,iter,res] = gmres(A, b)
@@ -174,15 +174,15 @@ M = eye(10, 10);
 [x,flag,err,iter,res] = gmres(A, b, 10, 1d-12, 20, M, zeros(10, 1))
 
        // If A is a matrix and M is a function
-       A=[ 94  0   0   0    0   28  0   0   32  0  
-     0   59  13  5    0   0   0   10  0   0  
-     0   13  72  34   2   0   0   0   0   65 
-     0   5   34  114  0   0   0   0   0   55 
-     0   0   2   0    70  0   28  32  12  0  
-     28  0   0   0    0   87  20  0   33  0  
-     0   0   0   0    28  20  71  39  0   0  
-     0   10  0   0    32  0   39  46  8   0  
-     32  0   0   0    12  33  0   8   82  11 
+       A=[ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
      0   0   65  55   0   0   0   0   11  100];
 b=ones(10,1);
 
@@ -195,46 +195,46 @@ endfunction
 
        // If A is a function and M is a matrix
        function y = Atimesx(x)
-       A=[ 94  0   0   0    0   28  0   0   32  0  
-     0   59  13  5    0   0   0   10  0   0  
-     0   13  72  34   2   0   0   0   0   65 
-     0   5   34  114  0   0   0   0   0   55 
-     0   0   2   0    70  0   28  32  12  0  
-     28  0   0   0    0   87  20  0   33  0  
-     0   0   0   0    28  20  71  39  0   0  
-     0   10  0   0    32  0   39  46  8   0  
-     32  0   0   0    12  33  0   8   82  11 
+       A=[ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
      0   0   65  55   0   0   0   0   11  100];
         y = A * x;
         endfunction
-        
+
         b = ones(10,1);
         M = eye(10, 10);
-        
+
         [x,flag,err,iter,res] = gmres(Atimesx, b)
-        
+
         [x,flag,err,iter,res] = gmres(Atimesx, b, 10, 1d-12, 20, M, zeros(10,1))
-        
+
         // If A and M are functions
         function y = Atimesx(x)
-       A=[ 94  0   0   0    0   28  0   0   32  0  
-     0   59  13  5    0   0   0   10  0   0  
-     0   13  72  34   2   0   0   0   0   65 
-     0   5   34  114  0   0   0   0   0   55 
-     0   0   2   0    70  0   28  32  12  0  
-     28  0   0   0    0   87  20  0   33  0  
-     0   0   0   0    28  20  71  39  0   0  
-     0   10  0   0    32  0   39  46  8   0  
-     32  0   0   0    12  33  0   8   82  11 
+       A=[ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
      0   0   65  55   0   0   0   0   11  100];
         y = A * x;
         endfunction
-        
+
         function y = Mtimesx(x)
 M = eye(10,10);
 y = M*x;
 endfunction
-        
+
         [x,flag,err,iter,res] = gmres(Atimesx, b, 10, 1d-12, 20, Mtimesx, zeros(10,1))
  ]]></programlisting>
     </refsection>
@@ -242,7 +242,7 @@ endfunction
         <title>See Also</title>
         <simplelist type="inline">
             <member>
-                <link linkend="pcg">pcg</link>
+                <link linkend="conjgrad">conjgrad</link>
             </member>
             <member>
                 <link linkend="qmr">qmr</link>
index caa0e09..2c7bfb9 100644 (file)
@@ -1,6 +1,7 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: made obsolete, merged in conjgrad
  * Copyright (C) XXXX-2008 - INRIA
  *
  * This file must be used under the terms of the CeCILL.
@@ -149,6 +150,11 @@ function y = A ( x )
     </refsection>
     <refsection>
         <title>Description</title>
+        <warning>
+            This function is obsolete.
+            It has been merged into <link linkend="conjgrad">conjgrad</link>.
+            <literal>pcg(A, b, ...) => conjgrad(A, b, "pcg", ...)</literal>.
+        </warning>
         <para>
             Solves the linear system <literal>Ax=b</literal> using the conjugate
             gradient method with or without preconditioning. The preconditionning
@@ -295,6 +301,9 @@ b=[101;11];
         <title>See Also</title>
         <simplelist type="inline">
             <member>
+                <link linkend="conjgrad">conjgrad</link>
+            </member>
+            <member>
                 <link linkend="backslash">backslash</link>
             </member>
             <member>
@@ -317,4 +326,15 @@ b=[101;11];
             ftp.cs.umn.edu/dept/users/saad/PS/all_ps.zip
         </para>
     </refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.5.0</revnumber>
+                <revdescription>
+                    Calling pcg(A, b, ...) is deprecated. conjgrad(A, b, "pcg", ...) should be used instead.
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
index 395b28f..e89934c 100644 (file)
@@ -2,11 +2,11 @@
 <!--
  * 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    
+ * are also available at
  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
  *
  -->
@@ -28,7 +28,7 @@
                 <term>A</term>
                 <listitem>
                     <para>
-                        matrix of size n-by-n or function.                
+                        matrix of size n-by-n or function.
                     </para>
                     <itemizedlist>
                         <listitem>
                                 <literal>function.</literal>If A is a function which returns <literal>A*x</literal>, it must
                                 have the following header :
                             </para>
-                            <programlisting role=""><![CDATA[ 
+                            <programlisting role=""><![CDATA[
 function y = A ( x )
  ]]></programlisting>
                             <para>
-                                If A is a function which returns <literal>A*x</literal> or <literal>A'*x</literal> depending t. 
-                                If <literal>t = "notransp"</literal>, the function returns <literal>A*x</literal>. 
+                                If A is a function which returns <literal>A*x</literal> or <literal>A'*x</literal> depending t.
+                                If <literal>t = "notransp"</literal>, the function returns <literal>A*x</literal>.
                                 If <literal>t = "transp"</literal>, the function returns <literal>A'*x</literal>. It must
                                 have the following header :
                             </para>
-                            <programlisting role=""><![CDATA[ 
+                            <programlisting role=""><![CDATA[
 function y = A ( x, t )
  ]]></programlisting>
                         </listitem>
@@ -64,7 +64,7 @@ function y = A ( x, t )
                     <para>
                         function returning <literal>A'*x</literal>. It must have the followinf header :
                     </para>
-                    <programlisting role=""><![CDATA[ 
+                    <programlisting role=""><![CDATA[
 function y = Ap ( x )
  ]]></programlisting>
                 </listitem>
@@ -107,8 +107,8 @@ function y = Ap ( x )
                 <term>M1p</term>
                 <listitem>
                     <para>
-                        must only be provided when <literal>M1</literal> is a function returning  <literal>M1*x</literal>. 
-                        In this case <literal>M1p</literal> is the function which returns <literal>M1'*x</literal>. 
+                        must only be provided when <literal>M1</literal> is a function returning  <literal>M1*x</literal>.
+                        In this case <literal>M1p</literal> is the function which returns <literal>M1'*x</literal>.
                     </para>
                 </listitem>
             </varlistentry>
@@ -138,7 +138,7 @@ function y = Ap ( x )
                 <term>M2p</term>
                 <listitem>
                     <para>
-                        must only be provided when <literal>M2</literal> is a function returning  <literal>M2*x</literal>. 
+                        must only be provided when <literal>M2</literal> is a function returning  <literal>M2*x</literal>.
                         In this case <literal>M2p</literal> is the function which returns <literal>M2'*x</literal>
                     </para>
                 </listitem>
@@ -214,17 +214,17 @@ function y = Ap ( x )
     </refsection>
     <refsection>
         <title>Examples</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
        // If A is a matrix
-A=[ 94  0   0   0    0   28  0   0   32  0  
-     0   59  13  5    0   0   0   10  0   0  
-     0   13  72  34   2   0   0   0   0   65 
-     0   5   34  114  0   0   0   0   0   55 
-     0   0   2   0    70  0   28  32  12  0  
-     28  0   0   0    0   87  20  0   33  0  
-     0   0   0   0    28  20  71  39  0   0  
-     0   10  0   0    32  0   39  46  8   0  
-     32  0   0   0    12  33  0   8   82  11 
+A=[ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
      0   0   65  55   0   0   0   0   11  100];
 b=ones(10,1);
 [x,flag,err,iter,res] = qmr(A, b)
@@ -233,15 +233,15 @@ b=ones(10,1);
 
        // If A is a function
        function y = Atimesx(x,t)
-       A=[ 94  0   0   0    0   28  0   0   32  0  
-     0   59  13  5    0   0   0   10  0   0  
-     0   13  72  34   2   0   0   0   0   65 
-     0   5   34  114  0   0   0   0   0   55 
-     0   0   2   0    70  0   28  32  12  0  
-     28  0   0   0    0   87  20  0   33  0  
-     0   0   0   0    28  20  71  39  0   0  
-     0   10  0   0    32  0   39  46  8   0  
-     32  0   0   0    12  33  0   8   82  11 
+       A=[ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
      0   0   65  55   0   0   0   0   11  100];
         if (t == 'notransp') then
         y = A*x;
@@ -249,45 +249,45 @@ b=ones(10,1);
         y = A'*x;
     end
        endfunction
-        
+
         [x,flag,err,iter,res] = qmr(Atimesx, b)
-        
+
         [x,flag,err,iter,res] = qmr(Atimesx, b, zeros(10,1), eye(10,10), eye(10,10), 10, 1d-12)
-        
+
         // OR
-        
+
         function y = funA(x)
-       A = [ 94  0   0   0    0   28  0   0   32  0  
-     0   59  13  5    0   0   0   10  0   0  
-     0   13  72  34   2   0   0   0   0   65 
-     0   5   34  114  0   0   0   0   0   55 
-     0   0   2   0    70  0   28  32  12  0  
-     28  0   0   0    0   87  20  0   33  0  
-     0   0   0   0    28  20  71  39  0   0  
-     0   10  0   0    32  0   39  46  8   0  
-     32  0   0   0    12  33  0   8   82  11 
+       A = [ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
      0   0   65  55   0   0   0   0   11  100];
         y = A*x
        endfunction
-       
+
         function y = funAp(x)
-       A = [ 94  0   0   0    0   28  0   0   32  0  
-     0   59  13  5    0   0   0   10  0   0  
-     0   13  72  34   2   0   0   0   0   65 
-     0   5   34  114  0   0   0   0   0   55 
-     0   0   2   0    70  0   28  32  12  0  
-     28  0   0   0    0   87  20  0   33  0  
-     0   0   0   0    28  20  71  39  0   0  
-     0   10  0   0    32  0   39  46  8   0  
-     32  0   0   0    12  33  0   8   82  11 
+       A = [ 94  0   0   0    0   28  0   0   32  0
+     0   59  13  5    0   0   0   10  0   0
+     0   13  72  34   2   0   0   0   0   65
+     0   5   34  114  0   0   0   0   0   55
+     0   0   2   0    70  0   28  32  12  0
+     28  0   0   0    0   87  20  0   33  0
+     0   0   0   0    28  20  71  39  0   0
+     0   10  0   0    32  0   39  46  8   0
+     32  0   0   0    12  33  0   8   82  11
      0   0   65  55   0   0   0   0   11  100];
         y = A'*x
        endfunction
-        
+
         [x,flag,err,iter,res] = qmr(funA, funAp, b)
-        
+
         [x,flag,err,iter,res] = qmr(funA, funAp, b, zeros(10,1), eye(10,10), eye(10,10), 10, 1d-12)
-        
+
         // If A is a matrix, M1 and M2 are functions
         function y = M1timesx(x,t)
         M1 = eye(10,10);
@@ -297,7 +297,7 @@ b=ones(10,1);
         y = M1'*x;
     end
        endfunction
-       
+
        function y = M2timesx(x,t)
         M2 = eye(10,10);
     if(t=="notransp") then
@@ -306,33 +306,33 @@ b=ones(10,1);
         y = M2'*x;
     end
        endfunction
-       
+
        [x,flag,err,iter,res] = qmr(A, b, zeros(10,1), M1timesx, M2timesx, 10, 1d-12)
-       
+
        // OR
-       
+
        function y = funM1(x)
        M1 = eye(10,10);
        y = M1*x;
        endfunction
-       
+
        function y = funM1p(x)
        M1 = eye(10,10);
        y = M1'*x;
        endfunction
-       
+
        function y = funM2(x)
        M2 = eye(10,10);
        y = M2*x;
        endfunction
-       
+
        function y = funM2p(x)
        M2 = eye(10,10);
        y = M2'*x;
        endfunction
-       
+
        [x,flag,err,iter,res] = qmr(A, b, zeros(10,1), funM1, funM1p, funM2, funM2p, 10, 1d-12)
-       
+
        // If A, M1, M2 are functions
        [x,flag,err,iter,res] = qmr(funA, funAp, b, zeros(10,1), funM1, funM1p, funM2, funM2p, 10, 1d-12)
        [x,flag,err,iter,res] = qmr(Atimesx, b, zeros(10,1), M1timesx, M2timesx, 10, 1d-12)
@@ -346,7 +346,7 @@ b=ones(10,1);
                 <link linkend="gmres">gmres</link>
             </member>
             <member>
-                <link linkend="pcg">pcg</link>
+                <link linkend="conjgrad">conjgrad</link>
             </member>
         </simplelist>
     </refsection>
@@ -357,7 +357,7 @@ b=ones(10,1);
                 <revnumber>5.4.0</revnumber>
                 <revdescription>
                     Calling qmr(A, Ap) is deprecated. qmr(A) should be used instead. The following function is an example :
-                    <programlisting role=""><![CDATA[ 
+                    <programlisting role=""><![CDATA[
 function y = A ( x, t )
 Amat = eye(2,2);
 if ( t== "notransp") then
diff --git a/scilab/modules/sparse/macros/%bicg.sci b/scilab/modules/sparse/macros/%bicg.sci
new file mode 100644 (file)
index 0000000..800fd31
--- /dev/null
@@ -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/scilab/modules/sparse/macros/%bicgstab.sci b/scilab/modules/sparse/macros/%bicgstab.sci
new file mode 100644 (file)
index 0000000..fe90693
--- /dev/null
@@ -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/scilab/modules/sparse/macros/%cgs.sci b/scilab/modules/sparse/macros/%cgs.sci
new file mode 100644 (file)
index 0000000..bec33e0
--- /dev/null
@@ -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/scilab/modules/sparse/macros/%pcg.sci b/scilab/modules/sparse/macros/%pcg.sci
new file mode 100644 (file)
index 0000000..be3162f
--- /dev/null
@@ -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/scilab/modules/sparse/macros/conjgrad.sci b/scilab/modules/sparse/macros/conjgrad.sci
new file mode 100644 (file)
index 0000000..07a7da9
--- /dev/null
@@ -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
index 98f0765..e08d691 100644 (file)
@@ -1,4 +1,6 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: tagged as obsolete,
+//                                                         use conjgrad(A,b,"pcg",...) instead
 // Copyright (C) 2009 - INRIA - Michael Baudin
 // Copyright (C) 2008 - INRIA - Michael Baudin
 // Copyright (C) 2006 - INRIA - Serge Steer
@@ -57,6 +59,7 @@
 //     Matlab's pcg command.
 //
 function [x, flag, resNorm, iter, resVec] = pcg(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
+    warnobsolete("conjgrad", "5.5.1");
     [lhs,rhs] = argn(0);
     if (rhs < 2),
         error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"pcg",2,7));
@@ -306,4 +309,3 @@ function [x, flag, resNorm, iter, resVec] = pcg(%A, %b, tol, maxIter, %M, %M2, x
         end
     end
 endfunction
-
index 195f347..24d1d31 100644 (file)
@@ -4,6 +4,7 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
+// <-- CLI SHELL MODE -->
 // <-- Non-regression test for bug 3310 -->
 //
 // <-- Bugzilla URL -->
@@ -15,5 +16,9 @@ A=[10,1;1,10];
 // Suppose that the rhs is in myrhs instead of b
 myrhs=[11;11];
 [xcomputed, flag, err, iter, res]=pcg(A,myrhs);
+WARNING: Feature pcg is obsolete.
+WARNING: Please use conjgrad instead.
+WARNING: This feature will be permanently removed in Scilab 5.5.1
+
 xexpected=[1;1];
 if norm(xcomputed-xexpected)>%eps then bugmes();quit;end
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad.dia.ref b/scilab/modules/sparse/tests/unit_tests/conjgrad.dia.ref
new file mode 100644 (file)
index 0000000..6708187
--- /dev/null
@@ -0,0 +1,349 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//-------------------------------------------------------
+// PCG
+// Test with 2 input arguments and all output arguments
+A = [10 1; 1 10];
+b = [11; 11];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg");
+xexpected = [1; 1];
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checkequal ( flag , 0);
+assert_checkequal ( err , 0);
+assert_checkequal ( iter, 1);
+// Test with 3 input arguments and all output arguments
+tol = 100*%eps;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checktrue ( err <= tol);
+// Test with 4 input arguments and all output arguments
+maxit = 10;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit);
+assert_checkalmostequal ( xcomputed,xexpected,%eps);
+// Test with 5 input arguments and all output arguments
+M1 = [1 0; 0 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 6 input arguments and all output arguments
+M2 = M1;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1, M2);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 7 input arguments and all output arguments
+x0 = [1; 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1, M2, x0);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with non-positional input parameters so that 0 iteration can happen
+A2 = [100 1; 1 10];
+b2 = [101; 11];
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="pcg", maxIter=0);
+assert_checkequal ( iter , 0 );
+// Test with non-positional input parameters so that 1 iteration is sufficient
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="pcg", tol=0.1);
+assert_checkequal ( iter , 1 );
+// Test with non-positional input parameters so that pre-conditionning is necessary
+A3 = [100 1; 1 0.0101];
+M = A3;
+[xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="pcg", tol, maxit,%M=M, maxIter=5, tol=%eps);
+assert_checkequal ( flag , 0 );
+// Test with non-positional input parameters so that good initialization generates 0 iterations
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "pcg", x0=[1;1]);
+assert_checkequal ( iter , 0 );
+// Test the special case where b = 0
+b3 = [0; 0];
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "pcg");
+xexpected2 = [0; 0];
+assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
+assert_checkequal ( flag , 0 );
+assert_checktrue ( err <= %eps );
+assert_checkequal ( iter , 0 );
+// Try a hard case where preconditionning is necessary
+// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
+A4 = [
+1.                          0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134
+0.4999999999999982236432    0.3333333333333319825620    0.2499999999999988897770    0.1999999999999990951682    0.1666666666666659080143
+0.3333333333333320380731    0.2499999999999990563104    0.1999999999999992617017    0.1666666666666660745477    0.1428571428571423496123
+0.2499999999999990285549    0.1999999999999993449684    0.1666666666666661855700    0.1428571428571424883902    0.1249999999999996808109
+0.1999999999999991506794    0.1666666666666660745477    0.1428571428571424328791    0.1249999999999996669331    0.11111111111111082739
+];
+b4 = ones(5, 1);
+M = A4;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "pcg", %M=M, tol=%eps, maxIter=3);
+expected = [
+5
+-120
+630
+-1120
+630 ];
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 2 );
+// Try a difficult case where preconditionning is necessary
+// Use two pre-conditionning matrices.
+// This is the Cholesky factorization of the matrix A : C = chol(A)
+C = [
+1.    0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134
+0.    0.288675134594810367528     0.2886751345948112557060    0.2598076211353305131624    0.2309401076758494653074
+0.    0.                          0.0745355992499937836104    0.1118033988749897594817    0.1277753129999876502421
+0.    0.                          0.                          0.0188982236504644136865    0.0377964473009222076683
+0.    0.                          0.                          0.                          0.0047619047619250291087
+];
+M1 = C';
+M2 = C;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "pcg", %M=M1, %M2=M2, tol=%eps, maxIter=10);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 2 );
+// Well-conditionned problem
+A5 = [ 94  0   0   0    0   28  0   0   32  0
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+b5 = ones(10, 1);
+[x, fail, err, iter, res] = conjgrad(A5, b5, "pcg", 1d-12, 10);
+assert_checkequal ( flag ,  0 );
+assert_checkequal ( iter , 10 );
+//-------------------------------------------------------
+// CGS
+// Test with 2 input arguments and all output arguments
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs");
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checkequal ( flag , 0);
+assert_checkequal ( err , 0);
+assert_checkequal ( iter, 1);
+// Test with 3 input arguments and all output arguments
+tol = 100*%eps;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checktrue ( err <= tol);
+// Test with 4 input arguments and all output arguments
+maxit = 10;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit);
+assert_checkalmostequal ( xcomputed,xexpected,%eps);
+// Test with 5 input arguments and all output arguments
+M1 = [1 0; 0 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 6 input arguments and all output arguments
+M2 = M1;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1, M2);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 7 input arguments and all output arguments
+x0 = [1; 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1, M2, x0);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with non-positional input parameters so that 0 iteration can happen
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="cgs", maxIter=0);
+assert_checkequal ( iter , 0 );
+// Test with non-positional input parameters so that 1 iteration is sufficient
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="cgs", tol=0.1);
+assert_checkequal ( iter , 1 );
+// Test with non-positional input parameters so that pre-conditionning is necessary
+M = A3;
+[xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="cgs", tol, maxit,%M=M, maxIter=5, tol=%eps);
+assert_checkequal ( flag , 0 );
+// Test with non-positional input parameters so that good initialization generates 0 iterations
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "cgs", x0=[1;1]);
+assert_checkequal ( iter , 0 );
+// Test the special case where b = 0
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "cgs");
+xexpected2 = [0; 0];
+assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
+assert_checkequal ( flag , 0 );
+assert_checktrue ( err <= %eps );
+assert_checkequal ( iter , 0 );
+// Try a hard case where preconditionning is necessary
+// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
+M = A4;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "cgs", %M=M, tol=%eps, maxIter=3);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 2 );
+// Try a difficult case where preconditionning is necessary
+// Use two pre-conditionning matrices.
+// This is the Cholesky factorization of the matrix A : C = chol(A)
+M1 = C';
+M2 = C;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "cgs", %M=M1, %M2=M2, tol=%eps, maxIter=10);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 10 );
+// Well-conditionned problem with a nonsymmetric matrix (top-right element)
+A5 = [ 94  0   0   0    0   28  0   0   32  20
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+[x, fail, err, iter, res] = conjgrad(A5, b5, "cgs", 1d-12, 15);
+assert_checkequal ( flag ,  0 );
+assert_checkequal ( iter , 11 );
+//-------------------------------------------------------
+// BICG
+// Test with 2 input arguments and all output arguments
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg");
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checkequal ( flag , 0);
+assert_checkequal ( err , 0);
+assert_checkequal ( iter, 1);
+// Test with 3 input arguments and all output arguments
+tol = 100*%eps;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checktrue ( err <= tol);
+// Test with 4 input arguments and all output arguments
+maxit = 10;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit);
+assert_checkalmostequal ( xcomputed,xexpected,%eps);
+// Test with 5 input arguments and all output arguments
+M1 = [1 0; 0 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 6 input arguments and all output arguments
+M2 = M1;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1, M2);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 7 input arguments and all output arguments
+x0 = [1; 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1, M2, x0);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with non-positional input parameters so that 0 iteration can happen
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicg", maxIter=0);
+assert_checkequal ( iter , 0 );
+// Test with non-positional input parameters so that 1 iteration is sufficient
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicg", tol=0.1);
+assert_checkequal ( iter , 1 );
+// Test with non-positional input parameters so that pre-conditionning is necessary
+M = A3;
+[xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="bicg", tol, maxit,%M=M, maxIter=5, tol=%eps);
+assert_checkequal ( flag , 0 );
+// Test with non-positional input parameters so that good initialization generates 0 iterations
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "bicg", x0=[1;1]);
+assert_checkequal ( iter , 0 );
+// Test the special case where b = 0
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "bicg");
+assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
+assert_checkequal ( flag , 0 );
+assert_checktrue ( err <= %eps );
+assert_checkequal ( iter , 0 );
+// Try a hard case where preconditionning is necessary
+// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
+M = A4;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicg", %M=M, tol=%eps, maxIter=3);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 2 );
+// Try a difficult case where preconditionning is necessary
+// Use two pre-conditionning matrices.
+// This is the Cholesky factorization of the matrix A : C = chol(A)
+M1 = C';
+M2 = C;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicg", %M=M1, %M2=M2, tol=%eps, maxIter=10);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 8 );
+// Well-conditionned problem with a nonsymmetric matrix (top-right element)
+[x, fail, err, iter, res] = conjgrad(A5, b5, "bicg", 1d-12, 15);
+assert_checkequal ( flag ,  0 );
+assert_checkequal ( iter , 10 );
+//-------------------------------------------------------
+// BICGSTAB
+// Test with 2 input arguments and all output arguments
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab");
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checkequal ( flag , 0);
+assert_checkequal ( err , 0);
+assert_checkequal ( iter, 1);
+// Test with 3 input arguments and all output arguments
+tol = 100*%eps;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checktrue ( err <= tol);
+// Test with 4 input arguments and all output arguments
+maxit = 10;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit);
+assert_checkalmostequal ( xcomputed,xexpected,%eps);
+// Test with 5 input arguments and all output arguments
+M1 = [1 0; 0 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 6 input arguments and all output arguments
+M2 = M1;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1, M2);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 7 input arguments and all output arguments
+x0 = [1; 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1, M2, x0);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with non-positional input parameters so that 0 iteration can happen
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicgstab", maxIter=0);
+assert_checkequal ( iter , 0 );
+// Test with non-positional input parameters so that 1 iteration is sufficient
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicgstab", tol=0.1);
+assert_checkequal ( iter , 1 );
+// Test with non-positional input parameters so that pre-conditionning is necessary
+M = A3;
+[xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="bicgstab", tol, maxit,%M=M, maxIter=5, tol=%eps);
+assert_checkequal ( flag , 0 );
+// Test with non-positional input parameters so that good initialization generates 0 iterations
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "bicgstab", x0=[1;1]);
+assert_checkequal ( iter , 0 );
+// Test the special case where b = 0
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "bicgstab");
+assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
+assert_checkequal ( flag , 0 );
+assert_checktrue ( err <= %eps );
+assert_checkequal ( iter , 0 );
+// Try a hard case where preconditionning is necessary
+// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
+M = A4;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicgstab", %M=M, tol=%eps, maxIter=3);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 1 );
+// Try a difficult case where preconditionning is necessary
+// Use two pre-conditionning matrices.
+// This is the Cholesky factorization of the matrix A : C = chol(A)
+M1 = C';
+M2 = C;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicgstab", %M=M1, %M2=M2, tol=%eps, maxIter=10);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 9 );
+// Well-conditionned problem with a nonsymmetric matrix (top-right element)
+[x, fail, err, iter, res] = conjgrad(A5, b5, "bicgstab", 1d-12, 15);
+assert_checkequal ( flag ,  0 );
+assert_checkequal ( iter , 10 );
+//-------------------------------------------------------
+// Error checks
+refMsg = msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),"conjgrad",2,7);
+assert_checkerror("conjgrad(A);", refMsg);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d.\n"),"conjgrad",1);
+assert_checkerror("conjgrad(""string"", b);", refMsg);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",1);
+assert_checkerror("conjgrad(ones(5, 4), b);", refMsg);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: A matrix expected.\n"),"conjgrad",2);
+assert_checkerror("conjgrad(A, ""string"");", refMsg);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: Column vector expected.\n"),"conjgrad",2);
+assert_checkerror("conjgrad(A, A);", refMsg);
+refMsg = msprintf(_("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1);
+assert_checkerror("conjgrad(A, b4);", refMsg);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: Single String expected.\n"),"conjgrad",3);
+assert_checkerror("conjgrad(A, b, 1d-12, 15);", refMsg);
+refMsg = msprintf(_("%s: Wrong value for input argument #%d: %s, %s, %s or %s expected.\n"),"conjgrad",3,"pcg","cgs","bicg","bicgstab");
+assert_checkerror("conjgrad(A, b, ""gmres"", 1d-12, 15);", refMsg);
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad.tst b/scilab/modules/sparse/tests/unit_tests/conjgrad.tst
new file mode 100644 (file)
index 0000000..6039e2a
--- /dev/null
@@ -0,0 +1,368 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+
+//-------------------------------------------------------
+// PCG
+
+// Test with 2 input arguments and all output arguments
+A = [10 1; 1 10];
+b = [11; 11];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg");
+xexpected = [1; 1];
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checkequal ( flag , 0);
+assert_checkequal ( err , 0);
+assert_checkequal ( iter, 1);
+// Test with 3 input arguments and all output arguments
+tol = 100*%eps;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checktrue ( err <= tol);
+// Test with 4 input arguments and all output arguments
+maxit = 10;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit);
+assert_checkalmostequal ( xcomputed,xexpected,%eps);
+// Test with 5 input arguments and all output arguments
+M1 = [1 0; 0 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 6 input arguments and all output arguments
+M2 = M1;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1, M2);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 7 input arguments and all output arguments
+x0 = [1; 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1, M2, x0);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with non-positional input parameters so that 0 iteration can happen
+A2 = [100 1; 1 10];
+b2 = [101; 11];
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="pcg", maxIter=0);
+assert_checkequal ( iter , 0 );
+// Test with non-positional input parameters so that 1 iteration is sufficient
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="pcg", tol=0.1);
+assert_checkequal ( iter , 1 );
+// Test with non-positional input parameters so that pre-conditionning is necessary
+A3 = [100 1; 1 0.0101];
+M = A3;
+[xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="pcg", tol, maxit,%M=M, maxIter=5, tol=%eps);
+assert_checkequal ( flag , 0 );
+// Test with non-positional input parameters so that good initialization generates 0 iterations
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "pcg", x0=[1;1]);
+assert_checkequal ( iter , 0 );
+// Test the special case where b = 0
+b3 = [0; 0];
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "pcg");
+xexpected2 = [0; 0];
+assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
+assert_checkequal ( flag , 0 );
+assert_checktrue ( err <= %eps );
+assert_checkequal ( iter , 0 );
+// Try a hard case where preconditionning is necessary
+// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
+A4 = [
+1.                          0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134
+0.4999999999999982236432    0.3333333333333319825620    0.2499999999999988897770    0.1999999999999990951682    0.1666666666666659080143
+0.3333333333333320380731    0.2499999999999990563104    0.1999999999999992617017    0.1666666666666660745477    0.1428571428571423496123
+0.2499999999999990285549    0.1999999999999993449684    0.1666666666666661855700    0.1428571428571424883902    0.1249999999999996808109
+0.1999999999999991506794    0.1666666666666660745477    0.1428571428571424328791    0.1249999999999996669331    0.11111111111111082739
+];
+b4 = ones(5, 1);
+M = A4;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "pcg", %M=M, tol=%eps, maxIter=3);
+expected = [
+5
+-120
+630
+-1120
+630 ];
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 2 );
+// Try a difficult case where preconditionning is necessary
+// Use two pre-conditionning matrices.
+// This is the Cholesky factorization of the matrix A : C = chol(A)
+C = [
+1.    0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134
+0.    0.288675134594810367528     0.2886751345948112557060    0.2598076211353305131624    0.2309401076758494653074
+0.    0.                          0.0745355992499937836104    0.1118033988749897594817    0.1277753129999876502421
+0.    0.                          0.                          0.0188982236504644136865    0.0377964473009222076683
+0.    0.                          0.                          0.                          0.0047619047619250291087
+];
+M1 = C';
+M2 = C;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "pcg", %M=M1, %M2=M2, tol=%eps, maxIter=10);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 2 );
+// Well-conditionned problem
+A5 = [ 94  0   0   0    0   28  0   0   32  0
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+b5 = ones(10, 1);
+[x, fail, err, iter, res] = conjgrad(A5, b5, "pcg", 1d-12, 10);
+assert_checkequal ( flag ,  0 );
+assert_checkequal ( iter , 10 );
+
+//-------------------------------------------------------
+// CGS
+
+// Test with 2 input arguments and all output arguments
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs");
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checkequal ( flag , 0);
+assert_checkequal ( err , 0);
+assert_checkequal ( iter, 1);
+// Test with 3 input arguments and all output arguments
+tol = 100*%eps;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checktrue ( err <= tol);
+// Test with 4 input arguments and all output arguments
+maxit = 10;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit);
+assert_checkalmostequal ( xcomputed,xexpected,%eps);
+// Test with 5 input arguments and all output arguments
+M1 = [1 0; 0 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 6 input arguments and all output arguments
+M2 = M1;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1, M2);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 7 input arguments and all output arguments
+x0 = [1; 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1, M2, x0);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with non-positional input parameters so that 0 iteration can happen
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="cgs", maxIter=0);
+assert_checkequal ( iter , 0 );
+// Test with non-positional input parameters so that 1 iteration is sufficient
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="cgs", tol=0.1);
+assert_checkequal ( iter , 1 );
+// Test with non-positional input parameters so that pre-conditionning is necessary
+M = A3;
+[xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="cgs", tol, maxit,%M=M, maxIter=5, tol=%eps);
+assert_checkequal ( flag , 0 );
+// Test with non-positional input parameters so that good initialization generates 0 iterations
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "cgs", x0=[1;1]);
+assert_checkequal ( iter , 0 );
+// Test the special case where b = 0
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "cgs");
+xexpected2 = [0; 0];
+assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
+assert_checkequal ( flag , 0 );
+assert_checktrue ( err <= %eps );
+assert_checkequal ( iter , 0 );
+// Try a hard case where preconditionning is necessary
+// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
+M = A4;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "cgs", %M=M, tol=%eps, maxIter=3);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 2 );
+// Try a difficult case where preconditionning is necessary
+// Use two pre-conditionning matrices.
+// This is the Cholesky factorization of the matrix A : C = chol(A)
+M1 = C';
+M2 = C;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "cgs", %M=M1, %M2=M2, tol=%eps, maxIter=10);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 10 );
+// Well-conditionned problem with a nonsymmetric matrix (top-right element)
+A5 = [ 94  0   0   0    0   28  0   0   32  20
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+[x, fail, err, iter, res] = conjgrad(A5, b5, "cgs", 1d-12, 15);
+assert_checkequal ( flag ,  0 );
+assert_checkequal ( iter , 11 );
+
+//-------------------------------------------------------
+// BICG
+
+// Test with 2 input arguments and all output arguments
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg");
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checkequal ( flag , 0);
+assert_checkequal ( err , 0);
+assert_checkequal ( iter, 1);
+// Test with 3 input arguments and all output arguments
+tol = 100*%eps;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checktrue ( err <= tol);
+// Test with 4 input arguments and all output arguments
+maxit = 10;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit);
+assert_checkalmostequal ( xcomputed,xexpected,%eps);
+// Test with 5 input arguments and all output arguments
+M1 = [1 0; 0 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 6 input arguments and all output arguments
+M2 = M1;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1, M2);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 7 input arguments and all output arguments
+x0 = [1; 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1, M2, x0);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with non-positional input parameters so that 0 iteration can happen
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicg", maxIter=0);
+assert_checkequal ( iter , 0 );
+// Test with non-positional input parameters so that 1 iteration is sufficient
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicg", tol=0.1);
+assert_checkequal ( iter , 1 );
+// Test with non-positional input parameters so that pre-conditionning is necessary
+M = A3;
+[xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="bicg", tol, maxit,%M=M, maxIter=5, tol=%eps);
+assert_checkequal ( flag , 0 );
+// Test with non-positional input parameters so that good initialization generates 0 iterations
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "bicg", x0=[1;1]);
+assert_checkequal ( iter , 0 );
+// Test the special case where b = 0
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "bicg");
+assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
+assert_checkequal ( flag , 0 );
+assert_checktrue ( err <= %eps );
+assert_checkequal ( iter , 0 );
+// Try a hard case where preconditionning is necessary
+// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
+M = A4;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicg", %M=M, tol=%eps, maxIter=3);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 2 );
+// Try a difficult case where preconditionning is necessary
+// Use two pre-conditionning matrices.
+// This is the Cholesky factorization of the matrix A : C = chol(A)
+M1 = C';
+M2 = C;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicg", %M=M1, %M2=M2, tol=%eps, maxIter=10);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 8 );
+// Well-conditionned problem with a nonsymmetric matrix (top-right element)
+[x, fail, err, iter, res] = conjgrad(A5, b5, "bicg", 1d-12, 15);
+assert_checkequal ( flag ,  0 );
+assert_checkequal ( iter , 10 );
+
+//-------------------------------------------------------
+// BICGSTAB
+
+// Test with 2 input arguments and all output arguments
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab");
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checkequal ( flag , 0);
+assert_checkequal ( err , 0);
+assert_checkequal ( iter, 1);
+// Test with 3 input arguments and all output arguments
+tol = 100*%eps;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+assert_checktrue ( err <= tol);
+// Test with 4 input arguments and all output arguments
+maxit = 10;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit);
+assert_checkalmostequal ( xcomputed,xexpected,%eps);
+// Test with 5 input arguments and all output arguments
+M1 = [1 0; 0 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 6 input arguments and all output arguments
+M2 = M1;
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1, M2);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with 7 input arguments and all output arguments
+x0 = [1; 1];
+[xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1, M2, x0);
+assert_checkalmostequal ( xcomputed , xexpected , %eps);
+// Test with non-positional input parameters so that 0 iteration can happen
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicgstab", maxIter=0);
+assert_checkequal ( iter , 0 );
+// Test with non-positional input parameters so that 1 iteration is sufficient
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicgstab", tol=0.1);
+assert_checkequal ( iter , 1 );
+// Test with non-positional input parameters so that pre-conditionning is necessary
+M = A3;
+[xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="bicgstab", tol, maxit,%M=M, maxIter=5, tol=%eps);
+assert_checkequal ( flag , 0 );
+// Test with non-positional input parameters so that good initialization generates 0 iterations
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "bicgstab", x0=[1;1]);
+assert_checkequal ( iter , 0 );
+// Test the special case where b = 0
+[xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "bicgstab");
+assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
+assert_checkequal ( flag , 0 );
+assert_checktrue ( err <= %eps );
+assert_checkequal ( iter , 0 );
+// Try a hard case where preconditionning is necessary
+// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
+M = A4;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicgstab", %M=M, tol=%eps, maxIter=3);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 1 );
+// Try a difficult case where preconditionning is necessary
+// Use two pre-conditionning matrices.
+// This is the Cholesky factorization of the matrix A : C = chol(A)
+M1 = C';
+M2 = C;
+[xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicgstab", %M=M1, %M2=M2, tol=%eps, maxIter=10);
+assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
+assert_checkequal ( flag , 0 );
+assert_checkequal ( iter , 9 );
+// Well-conditionned problem with a nonsymmetric matrix (top-right element)
+[x, fail, err, iter, res] = conjgrad(A5, b5, "bicgstab", 1d-12, 15);
+assert_checkequal ( flag ,  0 );
+assert_checkequal ( iter , 10 );
+
+
+
+//-------------------------------------------------------
+// Error checks
+
+refMsg = msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),"conjgrad",2,7);
+assert_checkerror("conjgrad(A);", refMsg);
+
+refMsg = msprintf(_("%s: Wrong type for input argument #%d.\n"),"conjgrad",1);
+assert_checkerror("conjgrad(""string"", b);", refMsg);
+
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",1);
+assert_checkerror("conjgrad(ones(5, 4), b);", refMsg);
+
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: A matrix expected.\n"),"conjgrad",2);
+assert_checkerror("conjgrad(A, ""string"");", refMsg);
+
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: Column vector expected.\n"),"conjgrad",2);
+assert_checkerror("conjgrad(A, A);", refMsg);
+
+refMsg = msprintf(_("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1);
+assert_checkerror("conjgrad(A, b4);", refMsg);
+
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: Single String expected.\n"),"conjgrad",3);
+assert_checkerror("conjgrad(A, b, 1d-12, 15);", refMsg);
+
+refMsg = msprintf(_("%s: Wrong value for input argument #%d: %s, %s, %s or %s expected.\n"),"conjgrad",3,"pcg","cgs","bicg","bicgstab");
+assert_checkerror("conjgrad(A, b, ""gmres"", 1d-12, 15);", refMsg);
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad_function.dia.ref b/scilab/modules/sparse/tests/unit_tests/conjgrad_function.dia.ref
new file mode 100644 (file)
index 0000000..5186917
--- /dev/null
@@ -0,0 +1,68 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
+// Copyright (C) 2008 - INRIA - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//------------------------------------------------------------------
+// PCG
+// Numerical tests
+// Case where A is given as a function computing the right-hand side
+b = [154.
+87.
+186.
+208.
+144.
+168.
+158.
+135.
+178.
+231.];
+function y=Afunction(x)
+    mymatrix=[ 94  0   0   0    0   28  0   0   32  0
+    0   59  13  5    0   0   0   10  0   0
+    0   13  72  34   2   0   0   0   0   65
+    0   5   34  114  0   0   0   0   0   55
+    0   0   2   0    70  0   28  32  12  0
+    28  0   0   0    0   87  20  0   33  0
+    0   0   0   0    28  20  71  39  0   0
+    0   10  0   0    32  0   39  46  8   0
+    32  0   0   0    12  33  0   8   82  11
+    0   0   65  55   0   0   0   0   11  100];
+    y=mymatrix*x
+endfunction
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Afunction,b,"pcg");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>10 then bugmes();quit;end
+if err > 10**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// CGS
+// CGS needs 11 iterations to converge
+[xcomputed, fail, err, iter, res]=conjgrad(Afunction,b,"cgs",maxIter=11);
+if norm(xcomputed-xexpected)>100**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>11 then bugmes();quit;end
+if err > 100**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// BICG
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Afunction,b,"bicg");
+if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>10 then bugmes();quit;end
+if err > 10**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// BICGSTAB
+// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
+[xcomputed, fail, err, iter, res]=conjgrad(Afunction,b,"bicgstab");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10000**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>8 then bugmes();quit;end
+if err > 1000**3*%eps then bugmes();quit;end
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad_function.tst b/scilab/modules/sparse/tests/unit_tests/conjgrad_function.tst
new file mode 100644 (file)
index 0000000..d6dadd4
--- /dev/null
@@ -0,0 +1,76 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
+// Copyright (C) 2008 - INRIA - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+
+//------------------------------------------------------------------
+// PCG
+
+// Numerical tests
+// Case where A is given as a function computing the right-hand side
+b = [154.
+87.
+186.
+208.
+144.
+168.
+158.
+135.
+178.
+231.];
+function y=Afunction(x)
+    mymatrix=[ 94  0   0   0    0   28  0   0   32  0
+    0   59  13  5    0   0   0   10  0   0
+    0   13  72  34   2   0   0   0   0   65
+    0   5   34  114  0   0   0   0   0   55
+    0   0   2   0    70  0   28  32  12  0
+    28  0   0   0    0   87  20  0   33  0
+    0   0   0   0    28  20  71  39  0   0
+    0   10  0   0    32  0   39  46  8   0
+    32  0   0   0    12  33  0   8   82  11
+    0   0   65  55   0   0   0   0   11  100];
+    y=mymatrix*x
+endfunction
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Afunction,b,"pcg");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>10 then pause,end
+if err > 10**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// CGS
+
+// CGS needs 11 iterations to converge
+[xcomputed, fail, err, iter, res]=conjgrad(Afunction,b,"cgs",maxIter=11);
+if norm(xcomputed-xexpected)>100**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>11 then pause,end
+if err > 100**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// BICG
+
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Afunction,b,"bicg");
+if norm(xcomputed-xexpected)>10**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>10 then pause,end
+if err > 10**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// BICGSTAB
+
+// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
+[xcomputed, fail, err, iter, res]=conjgrad(Afunction,b,"bicgstab");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10000**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>8 then pause,end
+if err > 1000**3*%eps then pause,end
@@ -1,9 +1,14 @@
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
 // Copyright (C) 2008 - INRIA - Michael Baudin
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//------------------------------------------------------------------
+// PCG
 // Numerical tests
 // Case where A is given as a function computing the right-hand side
 mymatrix=[ 94  0   0   0    0   28  0   0   32  0
@@ -27,7 +32,7 @@ b = [154.
 178.
 231.];
 function y=Atimesx(x,mymatrix)
-  y=mymatrix*x
+    y=mymatrix*x
 endfunction
 // With the default 10 iterations, the algorithm performs well
 Alist = list(Atimesx,mymatrix)
@@ -66,9 +71,33 @@ Alist = list(Atimesx,mymatrix)
     0.    
     11.   
     100.  
-[xcomputed, fail, err, iter, res]=pcg(Alist,b);
+[xcomputed, fail, err, iter, res]=conjgrad(Alist,b,"pcg");
 xexpected=ones(10,1);
 if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
 if fail<>0 then bugmes();quit;end
 if iter<>10 then bugmes();quit;end
 if err > 10**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// CGS
+// CGS needs 11 iterations to converge
+[xcomputed, fail, err, iter, res]=conjgrad(Alist,b,"cgs",maxIter=11);
+if norm(xcomputed-xexpected)>100**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>11 then bugmes();quit;end
+if err > 100**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// BICG
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Alist,b,"bicg");
+if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>10 then bugmes();quit;end
+if err > 10**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// BICGSTAB
+// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
+[xcomputed, fail, err, iter, res]=conjgrad(Alist,b,"bicgstab");
+if norm(xcomputed-xexpected)>10000**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>8 then bugmes();quit;end
+if err > 1000**3*%eps then bugmes();quit;end
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad_list.tst b/scilab/modules/sparse/tests/unit_tests/conjgrad_list.tst
new file mode 100644 (file)
index 0000000..e1c1af7
--- /dev/null
@@ -0,0 +1,76 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
+// Copyright (C) 2008 - INRIA - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+
+//------------------------------------------------------------------
+// PCG
+
+// Numerical tests
+// Case where A is given as a function computing the right-hand side
+mymatrix=[ 94  0   0   0    0   28  0   0   32  0
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+b = [154.
+87.
+186.
+208.
+144.
+168.
+158.
+135.
+178.
+231.];
+function y=Atimesx(x,mymatrix)
+    y=mymatrix*x
+endfunction
+// With the default 10 iterations, the algorithm performs well
+Alist = list(Atimesx,mymatrix)
+[xcomputed, fail, err, iter, res]=conjgrad(Alist,b,"pcg");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>10 then pause,end
+if err > 10**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// CGS
+
+// CGS needs 11 iterations to converge
+[xcomputed, fail, err, iter, res]=conjgrad(Alist,b,"cgs",maxIter=11);
+if norm(xcomputed-xexpected)>100**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>11 then pause,end
+if err > 100**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// BICG
+
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Alist,b,"bicg");
+if norm(xcomputed-xexpected)>10**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>10 then pause,end
+if err > 10**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// BICGSTAB
+
+// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
+[xcomputed, fail, err, iter, res]=conjgrad(Alist,b,"bicgstab");
+if norm(xcomputed-xexpected)>10000**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>8 then pause,end
+if err > 1000**3*%eps then pause,end
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad_numerical.dia.ref b/scilab/modules/sparse/tests/unit_tests/conjgrad_numerical.dia.ref
new file mode 100644 (file)
index 0000000..bacd6d0
--- /dev/null
@@ -0,0 +1,108 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
+// Copyright (C) 2008 - INRIA - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//------------------------------------------------------------------
+// PCG
+// Numerical tests
+//Well conditionned problem
+A=[ 94  0   0   0    0   28  0   0   32  0
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+b = [154.
+87.
+186.
+208.
+144.
+168.
+158.
+135.
+178.
+231.];
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"pcg");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>10 then bugmes();quit;end
+if err > 10**3*%eps then bugmes();quit;end
+// With a tolerance of 1.e-3, there are 5 iterations and the status is "success"
+tol=1.d-3;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"pcg",tol);
+if fail<>0 then bugmes();quit;end
+if iter>10 then bugmes();quit;end
+// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
+tol=%eps;
+maxIter = 5;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"pcg",tol,maxIter);
+if fail<>1 then bugmes();quit;end
+if iter<>maxIter then bugmes();quit;end
+//------------------------------------------------------------------
+// CGS
+// CGS needs 11 iterations to converge
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"cgs",maxIter=11);
+if norm(xcomputed-xexpected)>100**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>11 then bugmes();quit;end
+if err > 100**3*%eps then bugmes();quit;end
+// With a tolerance of 1.e-3, there are 3 iterations and the status is "success"
+tol=1.d-3;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"cgs",tol);
+if fail<>0 then bugmes();quit;end
+if iter>10 then bugmes();quit;end
+// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
+tol=%eps;
+maxIter = 5;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"cgs",tol,maxIter);
+if fail<>1 then bugmes();quit;end
+if iter<>maxIter then bugmes();quit;end
+//------------------------------------------------------------------
+// BICG
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicg");
+if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>10 then bugmes();quit;end
+if err > 10**3*%eps then bugmes();quit;end
+// With a tolerance of 1.e-3, there are 5 iterations and the status is "success"
+tol=1.d-3;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicg",tol);
+if fail<>0 then bugmes();quit;end
+if iter>10 then bugmes();quit;end
+// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
+tol=%eps;
+maxIter = 5;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicg",tol,maxIter);
+if fail<>1 then bugmes();quit;end
+if iter<>maxIter then bugmes();quit;end
+//------------------------------------------------------------------
+// BICGSTAB
+// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicgstab");
+if norm(xcomputed-xexpected)>10000**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>8 then bugmes();quit;end
+if err > 1000**3*%eps then bugmes();quit;end
+// With a tolerance of 1.e-3, there are 3 iterations and the status is "success"
+tol=1.d-3;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicgstab",tol);
+if fail<>0 then bugmes();quit;end
+if iter>10 then bugmes();quit;end
+// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
+tol=%eps;
+maxIter = 5;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicgstab",tol,maxIter);
+if fail<>1 then bugmes();quit;end
+if iter<>maxIter then bugmes();quit;end
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad_numerical.tst b/scilab/modules/sparse/tests/unit_tests/conjgrad_numerical.tst
new file mode 100644 (file)
index 0000000..a1e9313
--- /dev/null
@@ -0,0 +1,116 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
+// Copyright (C) 2008 - INRIA - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+
+//------------------------------------------------------------------
+// PCG
+
+// Numerical tests
+//Well conditionned problem
+A=[ 94  0   0   0    0   28  0   0   32  0
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+b = [154.
+87.
+186.
+208.
+144.
+168.
+158.
+135.
+178.
+231.];
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"pcg");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>10 then pause,end
+if err > 10**3*%eps then pause,end
+// With a tolerance of 1.e-3, there are 5 iterations and the status is "success"
+tol=1.d-3;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"pcg",tol);
+if fail<>0 then pause,end
+if iter>10 then pause,end
+// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
+tol=%eps;
+maxIter = 5;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"pcg",tol,maxIter);
+if fail<>1 then pause,end
+if iter<>maxIter then pause,end
+
+//------------------------------------------------------------------
+// CGS
+
+// CGS needs 11 iterations to converge
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"cgs",maxIter=11);
+if norm(xcomputed-xexpected)>100**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>11 then pause,end
+if err > 100**3*%eps then pause,end
+// With a tolerance of 1.e-3, there are 3 iterations and the status is "success"
+tol=1.d-3;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"cgs",tol);
+if fail<>0 then pause,end
+if iter>10 then pause,end
+// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
+tol=%eps;
+maxIter = 5;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"cgs",tol,maxIter);
+if fail<>1 then pause,end
+if iter<>maxIter then pause,end
+
+//------------------------------------------------------------------
+// BICG
+
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicg");
+if norm(xcomputed-xexpected)>10**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>10 then pause,end
+if err > 10**3*%eps then pause,end
+// With a tolerance of 1.e-3, there are 5 iterations and the status is "success"
+tol=1.d-3;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicg",tol);
+if fail<>0 then pause,end
+if iter>10 then pause,end
+// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
+tol=%eps;
+maxIter = 5;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicg",tol,maxIter);
+if fail<>1 then pause,end
+if iter<>maxIter then pause,end
+
+//------------------------------------------------------------------
+// BICGSTAB
+
+// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicgstab");
+if norm(xcomputed-xexpected)>10000**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>8 then pause,end
+if err > 1000**3*%eps then pause,end
+// With a tolerance of 1.e-3, there are 3 iterations and the status is "success"
+tol=1.d-3;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicgstab",tol);
+if fail<>0 then pause,end
+if iter>10 then pause,end
+// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
+tol=%eps;
+maxIter = 5;
+[xcomputed, fail, err, iter, res]=conjgrad(A,b,"bicgstab",tol,maxIter);
+if fail<>1 then pause,end
+if iter<>maxIter then pause,end
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad_sparse.dia.ref b/scilab/modules/sparse/tests/unit_tests/conjgrad_sparse.dia.ref
new file mode 100644 (file)
index 0000000..491bc24
--- /dev/null
@@ -0,0 +1,65 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
+// Copyright (C) 2008 - INRIA - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//------------------------------------------------------------------
+// PCG
+// Numerical tests
+// Case where A is sparse
+A=[ 94  0   0   0    0   28  0   0   32  0
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+b = [154.
+87.
+186.
+208.
+144.
+168.
+158.
+135.
+178.
+231.];
+Asparse = sparse(A);
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"pcg");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>10 then bugmes();quit;end
+if err > 10**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// CGS
+// CGS needs 11 iterations to converge
+[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"cgs",maxIter=11);
+if norm(xcomputed-xexpected)>100**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>11 then bugmes();quit;end
+if err > 100**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// BICG
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"bicg");
+if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>10 then bugmes();quit;end
+if err > 10**3*%eps then bugmes();quit;end
+//------------------------------------------------------------------
+// BICGSTAB
+// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
+[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"bicgstab");
+if norm(xcomputed-xexpected)>10000**3*%eps then bugmes();quit;end
+if fail<>0 then bugmes();quit;end
+if iter<>8 then bugmes();quit;end
+if err > 1000**3*%eps then bugmes();quit;end
diff --git a/scilab/modules/sparse/tests/unit_tests/conjgrad_sparse.tst b/scilab/modules/sparse/tests/unit_tests/conjgrad_sparse.tst
new file mode 100644 (file)
index 0000000..dfe585b
--- /dev/null
@@ -0,0 +1,73 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
+// Copyright (C) 2008 - INRIA - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+
+//------------------------------------------------------------------
+// PCG
+
+// Numerical tests
+// Case where A is sparse
+A=[ 94  0   0   0    0   28  0   0   32  0
+0   59  13  5    0   0   0   10  0   0
+0   13  72  34   2   0   0   0   0   65
+0   5   34  114  0   0   0   0   0   55
+0   0   2   0    70  0   28  32  12  0
+28  0   0   0    0   87  20  0   33  0
+0   0   0   0    28  20  71  39  0   0
+0   10  0   0    32  0   39  46  8   0
+32  0   0   0    12  33  0   8   82  11
+0   0   65  55   0   0   0   0   11  100];
+b = [154.
+87.
+186.
+208.
+144.
+168.
+158.
+135.
+178.
+231.];
+Asparse = sparse(A);
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"pcg");
+xexpected=ones(10,1);
+if norm(xcomputed-xexpected)>10**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>10 then pause,end
+if err > 10**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// CGS
+
+// CGS needs 11 iterations to converge
+[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"cgs",maxIter=11);
+if norm(xcomputed-xexpected)>100**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>11 then pause,end
+if err > 100**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// BICG
+
+// With the default 10 iterations, the algorithm performs well
+[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"bicg");
+if norm(xcomputed-xexpected)>10**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>10 then pause,end
+if err > 10**3*%eps then pause,end
+
+//------------------------------------------------------------------
+// BICGSTAB
+
+// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
+[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"bicgstab");
+if norm(xcomputed-xexpected)>10000**3*%eps then pause,end
+if fail<>0 then pause,end
+if iter<>8 then pause,end
+if err > 1000**3*%eps then pause,end
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg.dia.ref b/scilab/modules/sparse/tests/unit_tests/pcg.dia.ref
deleted file mode 100644 (file)
index d06beed..0000000
+++ /dev/null
@@ -1,147 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008-2009 - INRIA - Michael Baudin
-// Copyright (C) 2011 - DIGITEO - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// <-- CLI SHELL MODE -->
-// Test with 2 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-[xcomputed, flag, err, iter, res]=pcg(A,b);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-assert_checkequal ( flag , 0);
-assert_checkequal ( err , 0);
-assert_checkequal ( iter, 1);
-// Test with 3 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-assert_checktrue ( err <= tol);
-// Test with 4 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-maxit = 10;
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol,maxit);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed,xexpected,%eps);
-// Test with 5 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-maxit = 10;
-M1=[1,0;0,1];
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol,maxit,M1);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-// Test with 6 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-maxit = 10;
-M1=[1,0;0,1];
-M2=[1,0;0,1];
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol,maxit,M1,M2);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-// Test with 7 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-maxit = 10;
-M1=[1,0;0,1];
-M2=[1,0;0,1];
-x0=[1;1];
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol,maxit,M1,M2,x0);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-// Test with non-positionnal input parameters so that 0 iteration can happen
-A=[100,1;1,10];
-b=[101;11];
-[xcomputed, flag, err, iter, res]=pcg(A,b,maxIter=0);
-assert_checkequal ( iter , 0 );
-// Test with non-positionnal input parameters so that 1 iteration is sufficient
-A=[100,1;1,10];
-b=[101;11];
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol=0.1);
-assert_checkequal ( iter , 1 );
-// Test with non-positionnal input parameters so that pre-conditionning is necessary
-A=[100,1;1,0.0101];
-b=[101;11];
-M=A;
-[xcomputed, flag, err, iter, res]=pcg(A,b,%M=M,maxIter=3,tol=%eps);
-assert_checkequal ( flag , 0 );
-// Test with non-positionnal input parameters so that good initialization generates 0 iterations
-A=[100,1;1,10.];
-b=[101;11];
-[xcomputed, flag, err, iter, res]=pcg(A,b,x0=[1.;1.]);
-assert_checkequal ( iter , 0 );
-// Test the special case where b=0
-A=[100,1;1,10.];
-b=[0;0];
-[xcomputed, flag, err, iter, res]=pcg(A,b);
-xexpected=[0;0];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-assert_checkequal ( flag , 0 );
-assert_checktrue ( err <= %eps );
-assert_checkequal ( iter , 0 );
-// Try a hard case where preconditionning is necessary
-// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
-A = [
-    1.                          0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134  
-    0.4999999999999982236432    0.3333333333333319825620    0.2499999999999988897770    0.1999999999999990951682    0.1666666666666659080143  
-    0.3333333333333320380731    0.2499999999999990563104    0.1999999999999992617017    0.1666666666666660745477    0.1428571428571423496123  
-    0.2499999999999990285549    0.1999999999999993449684    0.1666666666666661855700    0.1428571428571424883902    0.1249999999999996808109  
-    0.1999999999999991506794    0.1666666666666660745477    0.1428571428571424328791    0.1249999999999996669331    0.11111111111111082739    
-];
-b = ones(5,1);
-M=A;
-[xcomputed, flag, err, iter, res]=pcg(A,b,%M=M,tol=%eps,maxIter = 3 );
-expected = [
-    5.
-  -120.
-    630.
-  -1120.
-    630.
-];
-assert_checkalmostequal ( xcomputed , expected, 1.e8 * %eps );
-assert_checkequal ( flag , 0 );
-assert_checkequal ( iter , 2 );
-// Try a difficult case where preconditionning is necessary
-// Use two pre-conditionning matrices.
-// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
-A = [
-    1.                          0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134  
-    0.4999999999999982236432    0.3333333333333319825620    0.2499999999999988897770    0.1999999999999990951682    0.1666666666666659080143  
-    0.3333333333333320380731    0.2499999999999990563104    0.1999999999999992617017    0.1666666666666660745477    0.1428571428571423496123  
-    0.2499999999999990285549    0.1999999999999993449684    0.1666666666666661855700    0.1428571428571424883902    0.1249999999999996808109  
-    0.1999999999999991506794    0.1666666666666660745477    0.1428571428571424328791    0.1249999999999996669331    0.11111111111111082739    
-];
-b = ones(5,1);
-// This is the cholesky factorization of the matrix A : C = chol(A)
-C = [
-    1.    0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134  
-    0.    0.288675134594810367528     0.2886751345948112557060    0.2598076211353305131624    0.2309401076758494653074  
-    0.    0.                          0.0745355992499937836104    0.1118033988749897594817    0.1277753129999876502421  
-    0.    0.                          0.                          0.0188982236504644136865    0.0377964473009222076683  
-    0.    0.                          0.                          0.                          0.0047619047619250291087  
-];
-M1 = C';
-M2 = C;
-[xcomputed, flag, err, iter, res]=pcg(A,b , %M=M1 , %M2 = M2 , tol=%eps,maxIter = 3 );
-expected = [
-      5.
-   -120.
-    630.
-  -1120.
-    630.
-];
-assert_checkalmostequal ( xcomputed , expected, 1.e8 * %eps );
-assert_checkequal ( flag , 0 );
-assert_checkequal ( iter , 2 );
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg.tst b/scilab/modules/sparse/tests/unit_tests/pcg.tst
deleted file mode 100644 (file)
index 2cbfe0f..0000000
+++ /dev/null
@@ -1,149 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008-2009 - INRIA - Michael Baudin
-// Copyright (C) 2011 - DIGITEO - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// <-- CLI SHELL MODE -->
-
-// Test with 2 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-[xcomputed, flag, err, iter, res]=pcg(A,b);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-assert_checkequal ( flag , 0);
-assert_checkequal ( err , 0);
-assert_checkequal ( iter, 1);
-// Test with 3 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-assert_checktrue ( err <= tol);
-// Test with 4 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-maxit = 10;
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol,maxit);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed,xexpected,%eps);
-// Test with 5 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-maxit = 10;
-M1=[1,0;0,1];
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol,maxit,M1);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-// Test with 6 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-maxit = 10;
-M1=[1,0;0,1];
-M2=[1,0;0,1];
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol,maxit,M1,M2);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-// Test with 7 input arguments and all output arguments
-A=[10,1;1,10];
-b=[11;11];
-tol = 100*%eps;
-maxit = 10;
-M1=[1,0;0,1];
-M2=[1,0;0,1];
-x0=[1;1];
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol,maxit,M1,M2,x0);
-xexpected=[1;1];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-// Test with non-positionnal input parameters so that 0 iteration can happen
-A=[100,1;1,10];
-b=[101;11];
-[xcomputed, flag, err, iter, res]=pcg(A,b,maxIter=0);
-assert_checkequal ( iter , 0 );
-// Test with non-positionnal input parameters so that 1 iteration is sufficient
-A=[100,1;1,10];
-b=[101;11];
-[xcomputed, flag, err, iter, res]=pcg(A,b,tol=0.1);
-assert_checkequal ( iter , 1 );
-// Test with non-positionnal input parameters so that pre-conditionning is necessary
-A=[100,1;1,0.0101];
-b=[101;11];
-M=A;
-[xcomputed, flag, err, iter, res]=pcg(A,b,%M=M,maxIter=3,tol=%eps);
-assert_checkequal ( flag , 0 );
-// Test with non-positionnal input parameters so that good initialization generates 0 iterations
-A=[100,1;1,10.];
-b=[101;11];
-[xcomputed, flag, err, iter, res]=pcg(A,b,x0=[1.;1.]);
-assert_checkequal ( iter , 0 );
-// Test the special case where b=0
-A=[100,1;1,10.];
-b=[0;0];
-[xcomputed, flag, err, iter, res]=pcg(A,b);
-xexpected=[0;0];
-assert_checkalmostequal ( xcomputed , xexpected , %eps);
-assert_checkequal ( flag , 0 );
-assert_checktrue ( err <= %eps );
-assert_checkequal ( iter , 0 );
-// Try a hard case where preconditionning is necessary
-// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
-A = [
-    1.                          0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134  
-    0.4999999999999982236432    0.3333333333333319825620    0.2499999999999988897770    0.1999999999999990951682    0.1666666666666659080143  
-    0.3333333333333320380731    0.2499999999999990563104    0.1999999999999992617017    0.1666666666666660745477    0.1428571428571423496123  
-    0.2499999999999990285549    0.1999999999999993449684    0.1666666666666661855700    0.1428571428571424883902    0.1249999999999996808109  
-    0.1999999999999991506794    0.1666666666666660745477    0.1428571428571424328791    0.1249999999999996669331    0.11111111111111082739    
-];
-b = ones(5,1);
-M=A;
-[xcomputed, flag, err, iter, res]=pcg(A,b,%M=M,tol=%eps,maxIter = 3 );
-expected = [
-    5.
-  -120.
-    630.
-  -1120.
-    630.
-];
-assert_checkalmostequal ( xcomputed , expected, 1.e8 * %eps );
-assert_checkequal ( flag , 0 );
-assert_checkequal ( iter , 2 );
-// Try a difficult case where preconditionning is necessary
-// Use two pre-conditionning matrices.
-// This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
-A = [
-    1.                          0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134  
-    0.4999999999999982236432    0.3333333333333319825620    0.2499999999999988897770    0.1999999999999990951682    0.1666666666666659080143  
-    0.3333333333333320380731    0.2499999999999990563104    0.1999999999999992617017    0.1666666666666660745477    0.1428571428571423496123  
-    0.2499999999999990285549    0.1999999999999993449684    0.1666666666666661855700    0.1428571428571424883902    0.1249999999999996808109  
-    0.1999999999999991506794    0.1666666666666660745477    0.1428571428571424328791    0.1249999999999996669331    0.11111111111111082739    
-];
-b = ones(5,1);
-// This is the cholesky factorization of the matrix A : C = chol(A)
-C = [
-    1.    0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134  
-    0.    0.288675134594810367528     0.2886751345948112557060    0.2598076211353305131624    0.2309401076758494653074  
-    0.    0.                          0.0745355992499937836104    0.1118033988749897594817    0.1277753129999876502421  
-    0.    0.                          0.                          0.0188982236504644136865    0.0377964473009222076683  
-    0.    0.                          0.                          0.                          0.0047619047619250291087  
-];
-M1 = C';
-M2 = C;
-[xcomputed, flag, err, iter, res]=pcg(A,b , %M=M1 , %M2 = M2 , tol=%eps,maxIter = 3 );
-expected = [
-      5.
-   -120.
-    630.
-  -1120.
-    630.
-];
-assert_checkalmostequal ( xcomputed , expected, 1.e8 * %eps );
-assert_checkequal ( flag , 0 );
-assert_checkequal ( iter , 2 );
-
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg_function.dia.ref b/scilab/modules/sparse/tests/unit_tests/pcg_function.dia.ref
deleted file mode 100644 (file)
index 105d9ae..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - INRIA - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// <-- CLI SHELL MODE -->
-// Numerical tests
-// Case where A is given as a function computing the right-hand side
-b = [154.
-87.
-186.
-208.
-144.
-168.
-158.
-135.
-178.
-231.];
-function y=Afunction(x)
-  mymatrix=[ 94  0   0   0    0   28  0   0   32  0
-  0   59  13  5    0   0   0   10  0   0
-  0   13  72  34   2   0   0   0   0   65
-  0   5   34  114  0   0   0   0   0   55
-  0   0   2   0    70  0   28  32  12  0
-  28  0   0   0    0   87  20  0   33  0
-  0   0   0   0    28  20  71  39  0   0
-  0   10  0   0    32  0   39  46  8   0
-  32  0   0   0    12  33  0   8   82  11
-  0   0   65  55   0   0   0   0   11  100];
-  y=mymatrix*x
-endfunction
-// With the default 10 iterations, the algorithm performs well
-[xcomputed, fail, err, iter, res]=pcg(Afunction,b);
-xexpected=ones(10,1);
-if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
-if fail<>0 then bugmes();quit;end
-if iter<>10 then bugmes();quit;end
-if err > 10**3*%eps then bugmes();quit;end
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg_function.tst b/scilab/modules/sparse/tests/unit_tests/pcg_function.tst
deleted file mode 100644 (file)
index 5a54ef8..0000000
+++ /dev/null
@@ -1,43 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - INRIA - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-
-// <-- CLI SHELL MODE -->
-
-// Numerical tests
-// Case where A is given as a function computing the right-hand side
-b = [154.
-87.
-186.
-208.
-144.
-168.
-158.
-135.
-178.
-231.];
-function y=Afunction(x)
-  mymatrix=[ 94  0   0   0    0   28  0   0   32  0
-  0   59  13  5    0   0   0   10  0   0
-  0   13  72  34   2   0   0   0   0   65
-  0   5   34  114  0   0   0   0   0   55
-  0   0   2   0    70  0   28  32  12  0
-  28  0   0   0    0   87  20  0   33  0
-  0   0   0   0    28  20  71  39  0   0
-  0   10  0   0    32  0   39  46  8   0
-  32  0   0   0    12  33  0   8   82  11
-  0   0   65  55   0   0   0   0   11  100];
-  y=mymatrix*x
-endfunction
-// With the default 10 iterations, the algorithm performs well
-[xcomputed, fail, err, iter, res]=pcg(Afunction,b);
-xexpected=ones(10,1);
-if norm(xcomputed-xexpected)>10**3*%eps then pause,end
-if fail<>0 then pause,end
-if iter<>10 then pause,end
-if err > 10**3*%eps then pause,end
-
-
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg_list.tst b/scilab/modules/sparse/tests/unit_tests/pcg_list.tst
deleted file mode 100644 (file)
index 94bd11b..0000000
+++ /dev/null
@@ -1,43 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - INRIA - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-
-// <-- CLI SHELL MODE -->
-
-// Numerical tests
-// Case where A is given as a function computing the right-hand side
-mymatrix=[ 94  0   0   0    0   28  0   0   32  0
-0   59  13  5    0   0   0   10  0   0
-0   13  72  34   2   0   0   0   0   65
-0   5   34  114  0   0   0   0   0   55
-0   0   2   0    70  0   28  32  12  0
-28  0   0   0    0   87  20  0   33  0
-0   0   0   0    28  20  71  39  0   0
-0   10  0   0    32  0   39  46  8   0
-32  0   0   0    12  33  0   8   82  11
-0   0   65  55   0   0   0   0   11  100];
-b = [154.
-87.
-186.
-208.
-144.
-168.
-158.
-135.
-178.
-231.];
-function y=Atimesx(x,mymatrix)
-  y=mymatrix*x
-endfunction
-// With the default 10 iterations, the algorithm performs well
-Alist = list(Atimesx,mymatrix)
-[xcomputed, fail, err, iter, res]=pcg(Alist,b);
-xexpected=ones(10,1);
-if norm(xcomputed-xexpected)>10**3*%eps then pause,end
-if fail<>0 then pause,end
-if iter<>10 then pause,end
-if err > 10**3*%eps then pause,end
-
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg_numerical.dia.ref b/scilab/modules/sparse/tests/unit_tests/pcg_numerical.dia.ref
deleted file mode 100644 (file)
index 3ac93be..0000000
+++ /dev/null
@@ -1,46 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - INRIA - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// Numerical tests
-//Well conditionned problem
-A=[ 94  0   0   0    0   28  0   0   32  0
-0   59  13  5    0   0   0   10  0   0
-0   13  72  34   2   0   0   0   0   65
-0   5   34  114  0   0   0   0   0   55
-0   0   2   0    70  0   28  32  12  0
-28  0   0   0    0   87  20  0   33  0
-0   0   0   0    28  20  71  39  0   0
-0   10  0   0    32  0   39  46  8   0
-32  0   0   0    12  33  0   8   82  11
-0   0   65  55   0   0   0   0   11  100];
-b = [154.
-87.
-186.
-208.
-144.
-168.
-158.
-135.
-178.
-231.];
-// With the default 10 iterations, the algorithm performs well
-[xcomputed, fail, err, iter, res]=pcg(A,b);
-xexpected=ones(10,1);
-if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
-if fail<>0 then bugmes();quit;end
-if iter<>10 then bugmes();quit;end
-if err > 10**3*%eps then bugmes();quit;end
-// With a tolerance of 1.e-3, there are 5 iterations and the status is "sucess"
-tol=1.d-3;
-[xcomputed, fail, err, iter, res]=pcg(A,b,tol);
-if fail<>0 then bugmes();quit;end
-if iter>10 then bugmes();quit;end
-// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
-tol=%eps;
-maxIter = 5;
-[xcomputed, fail, err, iter, res]=pcg(A,b,tol,maxIter);
-if fail<>1 then bugmes();quit;end
-if iter<>maxIter then bugmes();quit;end
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg_numerical.tst b/scilab/modules/sparse/tests/unit_tests/pcg_numerical.tst
deleted file mode 100644 (file)
index a934543..0000000
+++ /dev/null
@@ -1,50 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - INRIA - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-
-// <-- CLI SHELL MODE -->
-
-// Numerical tests
-//Well conditionned problem
-A=[ 94  0   0   0    0   28  0   0   32  0
-0   59  13  5    0   0   0   10  0   0
-0   13  72  34   2   0   0   0   0   65
-0   5   34  114  0   0   0   0   0   55
-0   0   2   0    70  0   28  32  12  0
-28  0   0   0    0   87  20  0   33  0
-0   0   0   0    28  20  71  39  0   0
-0   10  0   0    32  0   39  46  8   0
-32  0   0   0    12  33  0   8   82  11
-0   0   65  55   0   0   0   0   11  100];
-b = [154.
-87.
-186.
-208.
-144.
-168.
-158.
-135.
-178.
-231.];
-// With the default 10 iterations, the algorithm performs well
-[xcomputed, fail, err, iter, res]=pcg(A,b);
-xexpected=ones(10,1);
-if norm(xcomputed-xexpected)>10**3*%eps then pause,end
-if fail<>0 then pause,end
-if iter<>10 then pause,end
-if err > 10**3*%eps then pause,end
-// With a tolerance of 1.e-3, there are 5 iterations and the status is "sucess"
-tol=1.d-3;
-[xcomputed, fail, err, iter, res]=pcg(A,b,tol);
-if fail<>0 then pause,end
-if iter>10 then pause,end
-// With a tolerance of %eps but only 5 iterations allowed, the status is "fail"
-tol=%eps;
-maxIter = 5;
-[xcomputed, fail, err, iter, res]=pcg(A,b,tol,maxIter);
-if fail<>1 then pause,end
-if iter<>maxIter then pause,end
-
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg_sparse.dia.ref b/scilab/modules/sparse/tests/unit_tests/pcg_sparse.dia.ref
deleted file mode 100644 (file)
index 1e7c1ce..0000000
+++ /dev/null
@@ -1,36 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - INRIA - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// Numerical tests
-// Case where A is sparse
-A=[ 94  0   0   0    0   28  0   0   32  0
-0   59  13  5    0   0   0   10  0   0
-0   13  72  34   2   0   0   0   0   65
-0   5   34  114  0   0   0   0   0   55
-0   0   2   0    70  0   28  32  12  0
-28  0   0   0    0   87  20  0   33  0
-0   0   0   0    28  20  71  39  0   0
-0   10  0   0    32  0   39  46  8   0
-32  0   0   0    12  33  0   8   82  11
-0   0   65  55   0   0   0   0   11  100];
-b = [154.
-87.
-186.
-208.
-144.
-168.
-158.
-135.
-178.
-231.];
-Asparse = sparse(A);
-// With the default 10 iterations, the algorithm performs well
-[xcomputed, fail, err, iter, res]=pcg(Asparse,b);
-xexpected=ones(10,1);
-if norm(xcomputed-xexpected)>10**3*%eps then bugmes();quit;end
-if fail<>0 then bugmes();quit;end
-if iter<>10 then bugmes();quit;end
-if err > 10**3*%eps then bugmes();quit;end
diff --git a/scilab/modules/sparse/tests/unit_tests/pcg_sparse.tst b/scilab/modules/sparse/tests/unit_tests/pcg_sparse.tst
deleted file mode 100644 (file)
index 8e95d16..0000000
+++ /dev/null
@@ -1,40 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - INRIA - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-
-// <-- CLI SHELL MODE -->
-
-// Numerical tests
-// Case where A is sparse
-A=[ 94  0   0   0    0   28  0   0   32  0
-0   59  13  5    0   0   0   10  0   0
-0   13  72  34   2   0   0   0   0   65
-0   5   34  114  0   0   0   0   0   55
-0   0   2   0    70  0   28  32  12  0
-28  0   0   0    0   87  20  0   33  0
-0   0   0   0    28  20  71  39  0   0
-0   10  0   0    32  0   39  46  8   0
-32  0   0   0    12  33  0   8   82  11
-0   0   65  55   0   0   0   0   11  100];
-b = [154.
-87.
-186.
-208.
-144.
-168.
-158.
-135.
-178.
-231.];
-Asparse = sparse(A);
-// With the default 10 iterations, the algorithm performs well
-[xcomputed, fail, err, iter, res]=pcg(Asparse,b);
-xexpected=ones(10,1);
-if norm(xcomputed-xexpected)>10**3*%eps then pause,end
-if fail<>0 then pause,end
-if iter<>10 then pause,end
-if err > 10**3*%eps then pause,end
-