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