1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier: added cgs, bicg and bicgstab
4// Copyright (C) 2008 - INRIA - Michael Baudin
5//
6//  This file is distributed under the same license as the Scilab package.
7// =============================================================================
8//
9// <-- CLI SHELL MODE -->
10
11//------------------------------------------------------------------
12// PCG
13
14// Numerical tests
15// Case where A is sparse
16A=[ 94  0   0   0    0   28  0   0   32  0
170   59  13  5    0   0   0   10  0   0
180   13  72  34   2   0   0   0   0   65
190   5   34  114  0   0   0   0   0   55
200   0   2   0    70  0   28  32  12  0
2128  0   0   0    0   87  20  0   33  0
220   0   0   0    28  20  71  39  0   0
230   10  0   0    32  0   39  46  8   0
2432  0   0   0    12  33  0   8   82  11
250   0   65  55   0   0   0   0   11  100];
26b = [154.
2787.
28186.
29208.
30144.
31168.
32158.
33135.
34178.
35231.];
36Asparse = sparse(A);
37// With the default 10 iterations, the algorithm performs well
38[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"pcg");
39xexpected=ones(10,1);
40if norm(xcomputed-xexpected)>11**3*%eps then pause,end
41if fail<>0 then pause,end
42if iter<>10 then pause,end
43if err > 10**3*%eps then pause,end
44
45//------------------------------------------------------------------
46// CGS
47
48// CGS needs 11 iterations to converge
49[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"cgs",maxIter=11);
50if norm(xcomputed-xexpected)>100**3*%eps then pause,end
51if fail<>0 then pause,end
52if iter<>11 then pause,end
53if err > 100**3*%eps then pause,end
54
55//------------------------------------------------------------------
56// BICG
57
58// With the default 10 iterations, the algorithm performs well
59[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"bicg");
60if norm(xcomputed-xexpected)>11**3*%eps then pause,end
61if fail<>0 then pause,end
62if iter<>10 then pause,end
63if err > 10**3*%eps then pause,end
64
65//------------------------------------------------------------------
66// BICGSTAB
67
68// BICGSTAB only needs 8 iterations to converge to the required tol, but is less accurate on arrival.
69[xcomputed, fail, err, iter, res]=conjgrad(Asparse,b,"bicgstab");
70if norm(xcomputed-xexpected)>10000**3*%eps then pause,end
71if fail<>0 then pause,end
72if iter<>8 then pause,end
73if err > 1000**3*%eps then pause,end
74