Bug #14271 fixed - conjgrad() displayed an incorrect error message about number of...
[scilab.git] / scilab / modules / sparse / tests / unit_tests / conjgrad.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
4 //
5 //  This file is distributed under the same license as the Scilab package.
6 // =============================================================================
7 //
8 // <-- CLI SHELL MODE -->
9
10 //-------------------------------------------------------
11 // PCG
12
13 // Test with 2 input arguments and all output arguments
14 A = [10 1; 1 10];
15 b = [11; 11];
16 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg");
17 xexpected = [1; 1];
18 assert_checkalmostequal ( xcomputed , xexpected , %eps);
19 assert_checkequal ( flag , 0);
20 assert_checkequal ( err , 0);
21 assert_checkequal ( iter, 1);
22 // Test with 3 input arguments and all output arguments
23 tol = 100*%eps;
24 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol);
25 assert_checkalmostequal ( xcomputed , xexpected , %eps);
26 assert_checktrue ( err <= tol);
27 // Test with 4 input arguments and all output arguments
28 maxit = 10;
29 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit);
30 assert_checkalmostequal ( xcomputed,xexpected,%eps);
31 // Test with 5 input arguments and all output arguments
32 M1 = [1 0; 0 1];
33 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1);
34 assert_checkalmostequal ( xcomputed , xexpected , %eps);
35 // Test with 6 input arguments and all output arguments
36 M2 = M1;
37 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1, M2);
38 assert_checkalmostequal ( xcomputed , xexpected , %eps);
39 // Test with 7 input arguments and all output arguments
40 x0 = [1; 1];
41 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "pcg", tol, maxit, M1, M2, x0);
42 assert_checkalmostequal ( xcomputed , xexpected , %eps);
43 // Test with non-positional input parameters so that 0 iteration can happen
44 A2 = [100 1; 1 10];
45 b2 = [101; 11];
46 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="pcg", maxIter=0);
47 assert_checkequal ( iter , 0 );
48 // Test with non-positional input parameters so that 1 iteration is sufficient
49 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="pcg", tol=0.1);
50 assert_checkequal ( iter , 1 );
51 // Test with non-positional input parameters so that pre-conditionning is necessary
52 A3 = [100 1; 1 0.0101];
53 M = A3;
54 [xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="pcg", tol, maxit,%M=M, maxIter=5, tol=%eps);
55 assert_checkequal ( flag , 0 );
56 // Test with non-positional input parameters so that good initialization generates 0 iterations
57 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "pcg", x0=[1;1]);
58 assert_checkequal ( iter , 0 );
59 // Test the special case where b = 0
60 b3 = [0; 0];
61 [xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "pcg");
62 xexpected2 = [0; 0];
63 assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
64 assert_checkequal ( flag , 0 );
65 assert_checktrue ( err <= %eps );
66 assert_checkequal ( iter , 0 );
67 // Try a hard case where preconditionning is necessary
68 // This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
69 A4 = [
70 1.                          0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134
71 0.4999999999999982236432    0.3333333333333319825620    0.2499999999999988897770    0.1999999999999990951682    0.1666666666666659080143
72 0.3333333333333320380731    0.2499999999999990563104    0.1999999999999992617017    0.1666666666666660745477    0.1428571428571423496123
73 0.2499999999999990285549    0.1999999999999993449684    0.1666666666666661855700    0.1428571428571424883902    0.1249999999999996808109
74 0.1999999999999991506794    0.1666666666666660745477    0.1428571428571424328791    0.1249999999999996669331    0.11111111111111082739
75 ];
76 b4 = ones(5, 1);
77 M = A4;
78 [xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "pcg", %M=M, tol=%eps, maxIter=3);
79 expected = [
80 5
81 -120
82 630
83 -1120
84 630 ];
85 assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
86 assert_checkequal ( flag , 0 );
87 assert_checkequal ( iter , 2 );
88 // Try a difficult case where preconditionning is necessary
89 // Use two pre-conditionning matrices.
90 // This is the Cholesky factorization of the matrix A : C = chol(A)
91 C = [
92 1.    0.5000000000000001110223    0.3333333333333334258519    0.2500000000000000555112    0.2000000000000000666134
93 0.    0.288675134594810367528     0.2886751345948112557060    0.2598076211353305131624    0.2309401076758494653074
94 0.    0.                          0.0745355992499937836104    0.1118033988749897594817    0.1277753129999876502421
95 0.    0.                          0.                          0.0188982236504644136865    0.0377964473009222076683
96 0.    0.                          0.                          0.                          0.0047619047619250291087
97 ];
98 M1 = C';
99 M2 = C;
100 [xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "pcg", %M=M1, %M2=M2, tol=%eps, maxIter=10);
101 assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
102 assert_checkequal ( flag , 0 );
103 assert_checkequal ( iter , 2 );
104 // Well-conditionned problem
105 A5 = [ 94  0   0   0    0   28  0   0   32  0
106 0   59  13  5    0   0   0   10  0   0
107 0   13  72  34   2   0   0   0   0   65
108 0   5   34  114  0   0   0   0   0   55
109 0   0   2   0    70  0   28  32  12  0
110 28  0   0   0    0   87  20  0   33  0
111 0   0   0   0    28  20  71  39  0   0
112 0   10  0   0    32  0   39  46  8   0
113 32  0   0   0    12  33  0   8   82  11
114 0   0   65  55   0   0   0   0   11  100];
115 b5 = ones(10, 1);
116 [x, fail, err, iter, res] = conjgrad(A5, b5, "pcg", 1d-12, 10);
117 assert_checkequal ( flag ,  0 );
118 assert_checkequal ( iter , 10 );
119
120 //-------------------------------------------------------
121 // CGS
122
123 // Test with 2 input arguments and all output arguments
124 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs");
125 assert_checkalmostequal ( xcomputed , xexpected , %eps);
126 assert_checkequal ( flag , 0);
127 assert_checkequal ( err , 0);
128 assert_checkequal ( iter, 1);
129 // Test with 3 input arguments and all output arguments
130 tol = 100*%eps;
131 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol);
132 assert_checkalmostequal ( xcomputed , xexpected , %eps);
133 assert_checktrue ( err <= tol);
134 // Test with 4 input arguments and all output arguments
135 maxit = 10;
136 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit);
137 assert_checkalmostequal ( xcomputed,xexpected,%eps);
138 // Test with 5 input arguments and all output arguments
139 M1 = [1 0; 0 1];
140 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1);
141 assert_checkalmostequal ( xcomputed , xexpected , %eps);
142 // Test with 6 input arguments and all output arguments
143 M2 = M1;
144 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1, M2);
145 assert_checkalmostequal ( xcomputed , xexpected , %eps);
146 // Test with 7 input arguments and all output arguments
147 x0 = [1; 1];
148 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "cgs", tol, maxit, M1, M2, x0);
149 assert_checkalmostequal ( xcomputed , xexpected , %eps);
150 // Test with non-positional input parameters so that 0 iteration can happen
151 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="cgs", maxIter=0);
152 assert_checkequal ( iter , 0 );
153 // Test with non-positional input parameters so that 1 iteration is sufficient
154 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="cgs", tol=0.1);
155 assert_checkequal ( iter , 1 );
156 // Test with non-positional input parameters so that pre-conditionning is necessary
157 M = A3;
158 [xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="cgs", tol, maxit,%M=M, maxIter=5, tol=%eps);
159 assert_checkequal ( flag , 0 );
160 // Test with non-positional input parameters so that good initialization generates 0 iterations
161 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "cgs", x0=[1;1]);
162 assert_checkequal ( iter , 0 );
163 // Test the special case where b = 0
164 [xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "cgs");
165 xexpected2 = [0; 0];
166 assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
167 assert_checkequal ( flag , 0 );
168 assert_checktrue ( err <= %eps );
169 assert_checkequal ( iter , 0 );
170 // Try a hard case where preconditionning is necessary
171 // This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
172 M = A4;
173 [xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "cgs", %M=M, tol=%eps, maxIter=3);
174 assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
175 assert_checkequal ( flag , 0 );
176 assert_checkequal ( iter , 2 );
177 // Try a difficult case where preconditionning is necessary
178 // Use two pre-conditionning matrices.
179 // This is the Cholesky factorization of the matrix A : C = chol(A)
180 M1 = C';
181 M2 = C;
182 [xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "cgs", %M=M1, %M2=M2, tol=%eps, maxIter=10);
183 assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
184 assert_checkequal ( flag , 0 );
185 assert_checkequal ( iter , 10 );
186 // Well-conditionned problem with a nonsymmetric matrix (top-right element)
187 A5 = [ 94  0   0   0    0   28  0   0   32  20
188 0   59  13  5    0   0   0   10  0   0
189 0   13  72  34   2   0   0   0   0   65
190 0   5   34  114  0   0   0   0   0   55
191 0   0   2   0    70  0   28  32  12  0
192 28  0   0   0    0   87  20  0   33  0
193 0   0   0   0    28  20  71  39  0   0
194 0   10  0   0    32  0   39  46  8   0
195 32  0   0   0    12  33  0   8   82  11
196 0   0   65  55   0   0   0   0   11  100];
197 [x, fail, err, iter, res] = conjgrad(A5, b5, "cgs", 1d-12, 15);
198 assert_checkequal ( flag ,  0 );
199 assert_checkequal ( iter , 11 );
200
201 //-------------------------------------------------------
202 // BICG
203
204 // Test with 2 input arguments and all output arguments
205 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg");
206 assert_checkalmostequal ( xcomputed , xexpected , %eps);
207 assert_checkequal ( flag , 0);
208 assert_checkequal ( err , 0);
209 assert_checkequal ( iter, 1);
210 // Test with 3 input arguments and all output arguments
211 tol = 100*%eps;
212 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol);
213 assert_checkalmostequal ( xcomputed , xexpected , %eps);
214 assert_checktrue ( err <= tol);
215 // Test with 4 input arguments and all output arguments
216 maxit = 10;
217 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit);
218 assert_checkalmostequal ( xcomputed,xexpected,%eps);
219 // Test with 5 input arguments and all output arguments
220 M1 = [1 0; 0 1];
221 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1);
222 assert_checkalmostequal ( xcomputed , xexpected , %eps);
223 // Test with 6 input arguments and all output arguments
224 M2 = M1;
225 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1, M2);
226 assert_checkalmostequal ( xcomputed , xexpected , %eps);
227 // Test with 7 input arguments and all output arguments
228 x0 = [1; 1];
229 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicg", tol, maxit, M1, M2, x0);
230 assert_checkalmostequal ( xcomputed , xexpected , %eps);
231 // Test with non-positional input parameters so that 0 iteration can happen
232 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicg", maxIter=0);
233 assert_checkequal ( iter , 0 );
234 // Test with non-positional input parameters so that 1 iteration is sufficient
235 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicg", tol=0.1);
236 assert_checkequal ( iter , 1 );
237 // Test with non-positional input parameters so that pre-conditionning is necessary
238 M = A3;
239 [xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="bicg", tol, maxit,%M=M, maxIter=5, tol=%eps);
240 assert_checkequal ( flag , 0 );
241 // Test with non-positional input parameters so that good initialization generates 0 iterations
242 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "bicg", x0=[1;1]);
243 assert_checkequal ( iter , 0 );
244 // Test the special case where b = 0
245 [xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "bicg");
246 assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
247 assert_checkequal ( flag , 0 );
248 assert_checktrue ( err <= %eps );
249 assert_checkequal ( iter , 0 );
250 // Try a hard case where preconditionning is necessary
251 // This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
252 M = A4;
253 [xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicg", %M=M, tol=%eps, maxIter=3);
254 assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
255 assert_checkequal ( flag , 0 );
256 assert_checkequal ( iter , 2 );
257 // Try a difficult case where preconditionning is necessary
258 // Use two pre-conditionning matrices.
259 // This is the Cholesky factorization of the matrix A : C = chol(A)
260 M1 = C';
261 M2 = C;
262 [xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicg", %M=M1, %M2=M2, tol=%eps, maxIter=10);
263 assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
264 assert_checkequal ( flag , 0 );
265 assert_checkequal ( iter , 8 );
266 // Well-conditionned problem with a nonsymmetric matrix (top-right element)
267 [x, fail, err, iter, res] = conjgrad(A5, b5, "bicg", 1d-12, 15);
268 assert_checkequal ( flag ,  0 );
269 assert_checkequal ( iter , 10 );
270
271 //-------------------------------------------------------
272 // BICGSTAB
273
274 // Test with 2 input arguments and all output arguments
275 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab");
276 assert_checkalmostequal ( xcomputed , xexpected , %eps);
277 assert_checkequal ( flag , 0);
278 assert_checkequal ( err , 0);
279 assert_checkequal ( iter, 1);
280 // Test with 3 input arguments and all output arguments
281 tol = 100*%eps;
282 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol);
283 assert_checkalmostequal ( xcomputed , xexpected , %eps);
284 assert_checktrue ( err <= tol);
285 // Test with 4 input arguments and all output arguments
286 maxit = 10;
287 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit);
288 assert_checkalmostequal ( xcomputed,xexpected,%eps);
289 // Test with 5 input arguments and all output arguments
290 M1 = [1 0; 0 1];
291 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1);
292 assert_checkalmostequal ( xcomputed , xexpected , %eps);
293 // Test with 6 input arguments and all output arguments
294 M2 = M1;
295 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1, M2);
296 assert_checkalmostequal ( xcomputed , xexpected , %eps);
297 // Test with 7 input arguments and all output arguments
298 x0 = [1; 1];
299 [xcomputed, flag, err, iter, res] = conjgrad(A, b, "bicgstab", tol, maxit, M1, M2, x0);
300 assert_checkalmostequal ( xcomputed , xexpected , %eps);
301 // Test with non-positional input parameters so that 0 iteration can happen
302 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicgstab", maxIter=0);
303 assert_checkequal ( iter , 0 );
304 // Test with non-positional input parameters so that 1 iteration is sufficient
305 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2,  method="bicgstab", tol=0.1);
306 assert_checkequal ( iter , 1 );
307 // Test with non-positional input parameters so that pre-conditionning is necessary
308 M = A3;
309 [xcomputed, flag, err, iter, res] = conjgrad(A3, b2, method="bicgstab", tol, maxit,%M=M, maxIter=5, tol=%eps);
310 assert_checkequal ( flag , 0 );
311 // Test with non-positional input parameters so that good initialization generates 0 iterations
312 [xcomputed, flag, err, iter, res] = conjgrad(A2, b2, "bicgstab", x0=[1;1]);
313 assert_checkequal ( iter , 0 );
314 // Test the special case where b = 0
315 [xcomputed, flag, err, iter, res] = conjgrad(A2, b3, "bicgstab");
316 assert_checkalmostequal ( xcomputed , xexpected2 , %eps);
317 assert_checkequal ( flag , 0 );
318 assert_checktrue ( err <= %eps );
319 assert_checkequal ( iter , 0 );
320 // Try a hard case where preconditionning is necessary
321 // This is the Hilbert 5x5 matrix : A = 1/(testmatrix("hilb",5))
322 M = A4;
323 [xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicgstab", %M=M, tol=%eps, maxIter=3);
324 assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
325 assert_checkequal ( flag , 0 );
326 assert_checkequal ( iter , 1 );
327 // Try a difficult case where preconditionning is necessary
328 // Use two pre-conditionning matrices.
329 // This is the Cholesky factorization of the matrix A : C = chol(A)
330 M1 = C';
331 M2 = C;
332 [xcomputed, flag, err, iter, res] = conjgrad(A4, b4, "bicgstab", %M=M1, %M2=M2, tol=%eps, maxIter=10);
333 assert_checkalmostequal ( xcomputed , expected, 1.e8*%eps );
334 assert_checkequal ( flag , 0 );
335 assert_checkequal ( iter , 9 );
336 // Well-conditionned problem with a nonsymmetric matrix (top-right element)
337 [x, fail, err, iter, res] = conjgrad(A5, b5, "bicgstab", 1d-12, 15);
338 assert_checkequal ( flag ,  0 );
339 assert_checkequal ( iter , 10 );
340
341
342
343 //-------------------------------------------------------
344 // Error checks
345
346 refMsg = msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),"conjgrad",2,9);
347 assert_checkerror("conjgrad(A);", refMsg);
348
349 refMsg = msprintf(_("%s: Wrong type for input argument #%d.\n"),"conjgrad",1);
350 assert_checkerror("conjgrad(""string"", b);", refMsg);
351
352 refMsg = msprintf(_("%s: Wrong type for input argument #%d: Square matrix expected.\n"),"conjgrad",1);
353 assert_checkerror("conjgrad(ones(5, 4), b);", refMsg);
354
355 refMsg = msprintf(_("%s: Wrong type for input argument #%d: A matrix expected.\n"),"conjgrad",2);
356 assert_checkerror("conjgrad(A, ""string"");", refMsg);
357
358 refMsg = msprintf(_("%s: Wrong type for input argument #%d: Column vector expected.\n"),"conjgrad",2);
359 assert_checkerror("conjgrad(A, A);", refMsg);
360
361 refMsg = msprintf(_("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"conjgrad",2,1);
362 assert_checkerror("conjgrad(A, b4);", refMsg);
363
364 refMsg = msprintf(_("%s: Wrong type for input argument #%d: Single String expected.\n"),"conjgrad",3);
365 assert_checkerror("conjgrad(A, b, 1d-12, 15);", refMsg);
366
367 refMsg = msprintf(_("%s: Wrong value for input argument #%d: %s, %s, %s or %s expected.\n"),"conjgrad",3,"pcg","cgs","bicg","bicgstab");
368 assert_checkerror("conjgrad(A, b, ""gmres"", 1d-12, 15);", refMsg);