Bug #14271 fixed - conjgrad() displayed an incorrect error message about number of...
[scilab.git] / scilab / modules / sparse / macros / conjgrad.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: transformed into a gateway to
3 //                                                         propose pcg, cgs, bicg and bicgstab.
4 // Copyright (C) 2009 - INRIA - Michael Baudin
5 // Copyright (C) 2008 - INRIA - Michael Baudin
6 // Copyright (C) 2006 - INRIA - Serge Steer
7 // Copyright (C) 2005 - IRISA - Sage Group
8 //
9 // Copyright (C) 2012 - 2016 - Scilab Enterprises
10 //
11 // This file is hereby licensed under the terms of the GNU GPL v2.0,
12 // pursuant to article 5.3.4 of the CeCILL v.2.1.
13 // This file was originally licensed under the terms of the CeCILL v2.1,
14 // and continues to be available under such terms.
15 // For more information, see the COPYING file which you should have received
16 // along with this program.
17
18 //
19 // conjgrad --
20 //   This function regroups four methods from the "Conjugate Gradient family" to solve the linear system %Ax=b:
21 //     - PCG (Preconditioned Conjugate Gradient): A must be symmetric positive definite,
22 //     - CGS (preconditioned Conjugate Gradient Squared): A must be square,
23 //     - BICG (preconditioned BiConjugate Gradient): A must be square,
24 //     - BICGSTAB (preconditioned BiConjugate Gradient Stabilized): A must be square (default method).
25 //   If M is given, it is used as a preconditionning matrix.
26 //   If both M and M2 are given, the matrix M * M2 is used as a preconditionning
27 //   matrix.
28 //
29 // input   %A       REAL matrix or a function y=Ax(x) which computes y=%A*x for a given x
30 //         b        REAL right hand side vector
31 //         tol, optional      REAL error tolerance (default: 1e-8)
32 //         maxIter, optional  INTEGER maximum number of iterations (default: size(%b))
33 //         %M, optional       REAL preconditioner matrix (default: none)
34 //         %M2, optional      REAL preconditioner matrix (default: none)
35 //         x0, optional       REAL initial guess vector (default: the zero vector)
36 //         verbose, optional  INTEGER set to 1 to enable verbose logging (default : 0)
37 //
38 // output  x        REAL solution vector
39 //         flag     INTEGER: 0 = solution found to tolerance
40 //                           1 = no convergence given maxIter
41 //         resNorm  REAL final relative norm of the residual
42 //         iter     INTEGER number of iterations performed
43 //         resVec   REAL residual vector
44 //
45 // References
46 //
47 //   PCG
48 //     "Templates for the Solution of Linear Systems: Building Blocks
49 //     for Iterative Methods",
50 //     Barrett, Berry, Chan, Demmel, Donato, Dongarra, Eijkhout,
51 //     Pozo, Romine, and Van der Vorst, SIAM Publications, 1993
52 //     (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
53 //
54 //     "Iterative Methods for Sparse Linear Systems, Second Edition"
55 //     Saad, SIAM Publications, 2003
56 //     (ftp ftp.cs.umn.edu; cd dept/users/saad/PS; get all_ps.zip).
57 //
58 //     Golub and Van Loan, Matrix Computations
59 //
60 //   CGS
61 //     "CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems"
62 //     by Peter Sonneveld
63 //
64 //     http://epubs.siam.org/doi/abs/10.1137/0910004
65 //     http://dl.acm.org/citation.cfm?id=64888&preflayout=flat
66 //     http://mathworld.wolfram.com/ConjugateGradientSquaredMethod.html
67 //
68 //   BICG
69 //     "Numerical Recipes: The Art of Scientific Computing." (third ed.)
70 //     by William Press, Saul Teukolsky, William Vetterling, Brian Flannery.
71 //
72 //     http://apps.nrbook.com/empanel/index.html?pg=87
73 //     http://dl.acm.org/citation.cfm?doid=1874391.187410
74 //     http://mathworld.wolfram.com/BiconjugateGradientMethod.html
75 //
76 //   BICGSTAB
77 //     "Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems"
78 //     by Henk van der Vorst.
79 //
80 //     http://epubs.siam.org/doi/abs/10.1137/0913035
81 //     http://dl.acm.org/citation.cfm?id=131916.131930&coll=DL&dl=GUIDE&CFID=372773884&CFTOKEN=56630250
82 //     http://mathworld.wolfram.com/BiconjugateGradientStabilizedMethod.html
83 //
84 // Notes
85 //
86 //     The input / output arguments of this command are the same as
87 //     Matlab's pcg, cgs, bicg and bicgstab commands, augmented with the 'method' argument
88 //
89
90 function [x, flag, resNorm, iter, resVec] = conjgrad(%A, %b, method, tol, maxIter, %M, %M2, x0, verbose )
91
92     [lhs, rhs] = argn(0);
93
94     if rhs < 2 | rhs >9 then
95         error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"conjgrad",2,9));
96     end
97
98     if exists("method", "local") == 0 then
99         method = "bicgstab";
100     end
101     if exists("tol", "local") == 0 then
102         tol = 1e-8
103     end
104     if exists("maxIter", "local") == 0 then
105         maxIter = size(%b, 1)
106     end
107     if exists("%M", "local") == 0 then
108         %M = []
109     end
110     if exists("%M2", "local") == 0 then
111         %M2 = []
112     end
113     if exists("x0", "local") == 0 then
114         x0 = zeros(%b);
115     end
116     if exists("verbose", "local") == 0 then
117         verbose = 0;
118     end
119     if verbose == 1 then
120         printf(gettext("Arguments:\n"));
121         printf("  tol = "+string(tol)+"\n");
122         printf("  maxIter = "+string(maxIter)+"\n");
123         printf("  M = \n")
124         disp(%M)
125         printf("  M2 = \n");
126         disp(%M2)
127         printf("  x0 = \n");
128         disp(x0)
129         printf("  verbose = "+string(verbose)+"\n");
130     end
131     // Compute matrixType
132     select type(%A)
133     case 1 then
134         matrixType = 1;
135     case 5 then
136         matrixType = 1;
137     case 13 then
138         matrixType = 0;
139         Aargs = list()
140     case 15 then
141         Aargs = list(%A(2:$))
142         // Caution : modify the input argument %A !
143         %A = %A(1);
144         matrixType = 0;
145     else
146         error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",1));
147     end
148     // If %A is a matrix (dense or sparse)
149     if matrixType == 1 then
150         if size(%A,1) ~= size(%A,2) then
151             error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",1));
152         end
153     end
154     // Check right hand side %b
155     if type(%b) ~= 1
156         error(msprintf(gettext("%s: Wrong type for input argument #%d: A matrix expected.\n"),"conjgrad",2));
157     end
158     if size(%b,2) ~= 1 then
159         error(msprintf(gettext("%s: Wrong type for input argument #%d: Column vector expected.\n"),"conjgrad",2));
160     end
161     if matrixType == 1 then
162         if size(%b,1) ~= size(%A,1) then
163             error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1));
164         end
165     end
166     if matrixType == 1 then
167         if size(%b,1) ~= size(%A,1) then
168             error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1));
169         end
170     end
171     // Check method
172     if type(method) ~= 10 | size(method) ~= [1 1]
173         error(msprintf(gettext("%s: Wrong type for input argument #%d: Single String expected.\n"),"conjgrad",3));
174     elseif and(method ~= ["pcg" "cgs" "bicg" "bicgstab"]),
175         error(msprintf(gettext("%s: Wrong value for input argument #%d: %s, %s, %s or %s expected.\n"),"conjgrad",3,"pcg","cgs","bicg","bicgstab"));
176     end
177     // Check type of the error tolerance tol
178     if or(size(tol) ~= [1 1]) then
179         error(msprintf(gettext("%s: Wrong type for input argument #%d: Scalar expected.\n"),"conjgrad",4));
180     end
181     // Check the type of maximum number of iterations maxIter
182     if or(size(maxIter) ~= [1 1]) then
183         error(msprintf(gettext("%s: Wrong type for input argument #%d: Scalar expected.\n"),"conjgrad",5));
184     end
185     // Compute precondType
186     select type(%M)
187     case 1 then
188         // M is a matrix
189         // precondType = 0 if the M is empty
190         //             = 1 if the M is not empty
191         precondType = bool2s(size(%M,"*")>=1);
192     case 5 then
193         precondType = 1;
194     case 13 then
195         Margs = list()
196         precondType = 2;
197     case 15 then
198         Margs = list(%M(2:$))
199         // Caution : modify the input argument %M !
200         %M = %M(1);
201         precondType = 2;
202     else
203         error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",6));
204     end
205     if precondType == 1 then
206         if size(%M,1) ~= size(%M,2) then
207             error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",6));
208         end
209         if size(%M,1) ~= size(%b,1) then
210             error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",6,2));
211         end
212     end
213     // Compute precondBis
214     select type(%M2)
215     case 1 then
216         // M2 is a matrix
217         // precondBis = 0 if the M2 is empty
218         //            = 1 if the M2 is not empty
219         precondBis = bool2s(size(%M2,"*")>=1);
220     case 5 then
221         precondBis = 1;
222     case 13 then
223         M2args = list()
224         precondBis = 2;
225     case 15 then
226         M2args = list(%M2(2:$))
227         // Caution : modify the input argument %M2 !
228         %M2 = %M2(1);
229         // Caution : modify precondType again !
230         precondType = 2;
231     else
232         error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",7));
233     end
234     if precondBis == 1 then
235         if size(%M2,1) ~= size(%M2,2) then
236             error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",7));
237         end
238         if size(%M2,1) ~= size(%b,1) then
239             error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",7,2));
240         end
241     end
242     // Check size of the initial vector x0
243     if size(x0,2) ~= 1 then
244         error(msprintf(gettext("%s: Wrong value for input argument #%d: Column vector expected.\n"),"conjgrad",8));
245     end
246     if size(x0,1) ~= size(%b,1) then
247         error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",8,2));
248     end
249
250     // ------------
251     // Computations
252     // ------------
253     select method
254     case "pcg"
255         [x, resNorm, iter, resVec] = %pcg(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
256     case "cgs"
257         [x, resNorm, iter, resVec] = %cgs(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
258     case "bicg"
259         [x, resNorm, iter, resVec] = %bicg(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
260     else // "bicgstab"
261         [x, resNorm, iter, resVec] = %bicgstab(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
262     end
263
264     // Test for convergence
265     if resNorm > tol then
266         if verbose == 1 then
267             printf(gettext("Final residual = %s > tol =%s\n"),string(resNorm),string(tol));
268             printf(gettext("Algorithm fails\n"));
269         end
270         flag = 1;
271         if lhs < 2 then
272             warning(msprintf(gettext("%s: Convergence error.\n"),"conjgrad"));
273         end
274     else
275         flag = 0;
276         if verbose == 1 then
277             printf(gettext("Algorithm pass\n"));
278         end
279     end
280
281 endfunction