Add the <-- JVM NOT MANDATORY --> TAG
[scilab.git] / scilab / modules / umfpack / tests / unit_tests / umf_lufact.tst
1 // ============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2007-2008 - Bruno PINCON
4 //
5 //  This file is distributed under the same license as the Scilab package.
6 // ============================================================================
7
8 // <-- JVM NOT MANDATORY -->
9
10 // this is the small linear test system from UMFPACK
11 // whom solution must be [1;2;3;4;5]
12 A = sparse( [ 2  3  0  0  0;
13               3  0  4  0  6;
14               0 -1 -3  2  0;
15               0  0  1  0  0;
16               0  4  2  0  1] );
17 b = [8 ; 45; -3; 3; 19];
18 Lup = umf_lufact(A);
19 x = umf_lusolve(Lup,b);
20
21 if or( (x-[ 1; 2; 3; 4; 5 ]) > %eps) then pause, end
22
23 // solve now A'x=b
24 x = umf_lusolve(Lup,b,"A''x=b");
25 if norm(A'*x - b) <> 0, then pause, end
26
27 // don't forget to clear memory with
28 umf_ludel(Lup)
29
30 // a real (but small)  example
31 // first load a sparse matrix
32 [A] = ReadHBSparse(SCI+"/modules/umfpack/examples/arc130.rua");
33 // compute the factorisation
34 Lup = umf_lufact(A);
35 b = rand(size(A,1),1); // a random rhs
36 // use umf_lusolve for solving Ax=b
37 x = umf_lusolve(Lup,b);
38 firstNorm=norm(A*x - b);
39
40 // now the same thing with iterative refiment
41 x = umf_lusolve(Lup,b,"Ax=b",A);
42 secondNorm=norm(A*x - b);
43
44 if firstNorm <> secondNorm then pause, end
45
46 // don't forget to clear memory
47 umf_ludel(Lup)