f21695790df9c99fe66ef642894fa454e69655c9
[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 then
95         error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"conjgrad",2,7));
96     end
97     if rhs > 8 then
98         error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"conjgrad",2,7));
99     end
100     if exists("method", "local") == 0 then
101         method = "bicgstab";
102     end
103     if exists("tol", "local") == 0 then
104         tol = 1e-8
105     end
106     if exists("maxIter", "local") == 0 then
107         maxIter = size(%b, 1)
108     end
109     if exists("%M", "local") == 0 then
110         %M = []
111     end
112     if exists("%M2", "local") == 0 then
113         %M2 = []
114     end
115     if exists("x0", "local") == 0 then
116         x0 = zeros(%b);
117     end
118     if exists("verbose", "local") == 0 then
119         verbose = 0;
120     end
121     if verbose == 1 then
122         printf(gettext("Arguments:\n"));
123         printf("  tol = "+string(tol)+"\n");
124         printf("  maxIter = "+string(maxIter)+"\n");
125         printf("  M = \n")
126         disp(%M)
127         printf("  M2 = \n");
128         disp(%M2)
129         printf("  x0 = \n");
130         disp(x0)
131         printf("  verbose = "+string(verbose)+"\n");
132     end
133     // Compute matrixType
134     select type(%A)
135     case 1 then
136         matrixType = 1;
137     case 5 then
138         matrixType = 1;
139     case 13 then
140         matrixType = 0;
141         Aargs = list()
142     case 15 then
143         Aargs = list(%A(2:$))
144         // Caution : modify the input argument %A !
145         %A = %A(1);
146         matrixType = 0;
147     else
148         error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",1));
149     end
150     // If %A is a matrix (dense or sparse)
151     if matrixType == 1 then
152         if size(%A,1) ~= size(%A,2) then
153             error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",1));
154         end
155     end
156     // Check right hand side %b
157     if type(%b) ~= 1
158         error(msprintf(gettext("%s: Wrong type for input argument #%d: A matrix expected.\n"),"conjgrad",2));
159     end
160     if size(%b,2) ~= 1 then
161         error(msprintf(gettext("%s: Wrong type for input argument #%d: Column vector expected.\n"),"conjgrad",2));
162     end
163     if matrixType == 1 then
164         if size(%b,1) ~= size(%A,1) then
165             error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1));
166         end
167     end
168     if matrixType == 1 then
169         if size(%b,1) ~= size(%A,1) then
170             error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1));
171         end
172     end
173     // Check method
174     if type(method) ~= 10 | size(method) ~= [1 1]
175         error(msprintf(gettext("%s: Wrong type for input argument #%d: Single String expected.\n"),"conjgrad",3));
176     elseif and(method ~= ["pcg" "cgs" "bicg" "bicgstab"]),
177         error(msprintf(gettext("%s: Wrong value for input argument #%d: %s, %s, %s or %s expected.\n"),"conjgrad",3,"pcg","cgs","bicg","bicgstab"));
178     end
179     // Check type of the error tolerance tol
180     if or(size(tol) ~= [1 1]) then
181         error(msprintf(gettext("%s: Wrong type for input argument #%d: Scalar expected.\n"),"conjgrad",4));
182     end
183     // Check the type of maximum number of iterations maxIter
184     if or(size(maxIter) ~= [1 1]) then
185         error(msprintf(gettext("%s: Wrong type for input argument #%d: Scalar expected.\n"),"conjgrad",5));
186     end
187     // Compute precondType
188     select type(%M)
189     case 1 then
190         // M is a matrix
191         // precondType = 0 if the M is empty
192         //             = 1 if the M is not empty
193         precondType = bool2s(size(%M,"*")>=1);
194     case 5 then
195         precondType = 1;
196     case 13 then
197         Margs = list()
198         precondType = 2;
199     case 15 then
200         Margs = list(%M(2:$))
201         // Caution : modify the input argument %M !
202         %M = %M(1);
203         precondType = 2;
204     else
205         error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",6));
206     end
207     if precondType == 1 then
208         if size(%M,1) ~= size(%M,2) then
209             error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",6));
210         end
211         if size(%M,1) ~= size(%b,1) then
212             error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",6,2));
213         end
214     end
215     // Compute precondBis
216     select type(%M2)
217     case 1 then
218         // M2 is a matrix
219         // precondBis = 0 if the M2 is empty
220         //            = 1 if the M2 is not empty
221         precondBis = bool2s(size(%M2,"*")>=1);
222     case 5 then
223         precondBis = 1;
224     case 13 then
225         M2args = list()
226         precondBis = 2;
227     case 15 then
228         M2args = list(%M2(2:$))
229         // Caution : modify the input argument %M2 !
230         %M2 = %M2(1);
231         // Caution : modify precondType again !
232         precondType = 2;
233     else
234         error(msprintf(gettext("%s: Wrong type for input argument #%d.\n"),"conjgrad",7));
235     end
236     if precondBis == 1 then
237         if size(%M2,1) ~= size(%M2,2) then
238             error(msprintf(gettext("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",7));
239         end
240         if size(%M2,1) ~= size(%b,1) then
241             error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",7,2));
242         end
243     end
244     // Check size of the initial vector x0
245     if size(x0,2) ~= 1 then
246         error(msprintf(gettext("%s: Wrong value for input argument #%d: Column vector expected.\n"),"conjgrad",8));
247     end
248     if size(x0,1) ~= size(%b,1) then
249         error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",8,2));
250     end
251
252     // ------------
253     // Computations
254     // ------------
255     select method
256     case "pcg"
257         [x, resNorm, iter, resVec] = %pcg(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
258     case "cgs"
259         [x, resNorm, iter, resVec] = %cgs(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
260     case "bicg"
261         [x, resNorm, iter, resVec] = %bicg(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
262     else // "bicgstab"
263         [x, resNorm, iter, resVec] = %bicgstab(%A, %b, tol, maxIter, %M, %M2, x0, verbose )
264     end
265
266     // Test for convergence
267     if resNorm > tol then
268         if verbose == 1 then
269             printf(gettext("Final residual = %s > tol =%s\n"),string(resNorm),string(tol));
270             printf(gettext("Algorithm fails\n"));
271         end
272         flag = 1;
273         if lhs < 2 then
274             warning(msprintf(gettext("%s: Convergence error.\n"),"conjgrad"));
275         end
276     else
277         flag = 0;
278         if verbose == 1 then
279             printf(gettext("Algorithm pass\n"));
280         end
281     end
282
283 endfunction