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