fix Demo UMPFACK
[scilab.git] / scilab / modules / umfpack / demos / sparse_matrices.sce
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright Bruno Pinçon, ESIAL-IECN, Inria CORIDA project
3 //   <bruno.pincon@iecn.u-nancy.fr>
4 //
5 // This file is released under the 3-clause BSD license. See COPYING-BSD.
6
7 // The simple demo for the umfpack interface & others
8 // sparse utility (scispt toolbox)
9
10 function demo_sparse_matrices()
11
12     format(15);
13     cd "SCI/modules/umfpack/demos/"
14
15     disp(_("To load a sparse matrix stored in the Harwell-Boeing format in Scilab, you may use the function ReadHBSparse."));
16     disp(_("4 matrices will be loaded for the demo: they come from the University of Florida Sparse Matrix Collection:"));
17     printf("     www.cise.ufl.edu/research/sparse/matrices\n"+..
18     "     http://math.nist.gov/MatrixMarket/matrices.html");
19     disp(_("    maintained by Tim Davis, author of UMFPACK"));
20
21     mode(1)
22     [A1,descr1,ref1,mtype1] = ReadHBSparse("arc130.rua");
23     [A2,descr2,ref2,mtype2] = ReadHBSparse("bcsstk24.rsa");
24     [A3,descr3,ref3,mtype3] = ReadHBSparse("ex14.rua");
25     [A4,descr4,ref4,mtype4] = ReadHBSparse("young1c.csa");
26     mode(-1)
27
28     disp(_("To see the pattern of non zeros elements, you may use the function PlotSparse."));
29     disp(_("Here we will display the 4 matrices"));
30     halt(_("\nPress Return to continue...\n"));
31
32     mode(1)
33     g = figure("visible", "off", "backgroundcolor", [1,1,1]);
34     old_winsid = winsid();
35     num = max(old_winsid);
36     g.figure_size = [800 800];
37     subplot(221)
38     PlotSparse(A1,"y+")
39     xtitle(ref1+"."+mtype1+": "+descr1)
40
41     subplot(222)
42     PlotSparse(A2)
43     xtitle(ref2+"."+mtype2+": "+descr2)
44
45     subplot(223)
46     PlotSparse(A3,"c.")
47     xtitle(ref3+"."+mtype3+": "+descr3)
48
49     subplot(224)
50     PlotSparse(A4,"r.")
51     xtitle(ref4+"."+mtype4+": "+descr4)
52     set(g, "visible", "on");
53     mode(-1)
54
55     disp(_("Now, using umfpack to solve some linear systems:"));
56     disp("     b1 = rand(size(A1,1),1);    "+_("-> to create a rhs"));
57     disp("     x1 = umfpack(A1,''\'',b1);    "+_("-> to solve A1*x1 = b1"));
58     disp("     norm_res = norm(A1*x1-b1);  "+_("-> norm of the residual"));
59     disp(_("this is done for the 4 matrices A1, A2, A3 and A4."));
60     halt(_("\nPress Return to continue...\n"));
61
62     mode(1)
63     n1 = size(A1,1);
64     b1 = rand(n1,1);
65     tic(); x1 = umfpack(A1,"\",b1); t1=toc();
66     norm_res1 = norm(A1*x1-b1);
67
68     n2 = size(A2,1);
69     b2 = rand(n2,1);
70     tic(); x2 = umfpack(A2,"\",b2); t2=toc();
71     norm_res2 = norm(A2*x2-b2);
72
73     n3 = size(A3,1);
74     b3 = rand(n3,1);
75     tic(); x3 = umfpack(A3,"\",b3); t3=toc();
76     norm_res3 =  norm(A3*x3-b3);
77
78     n4 = size(A4,1);
79     b4 = rand(n4,1);
80     tic(); x4 = umfpack(A4,"\",b4); t4=toc();
81     norm_res4 = norm(A4*x4-b4);
82     mode(-1)
83
84     disp(_("A1 ("+descr1+"): order = "+ string(n1) + " nnz = " + string(nnz(A1))));
85     msgNorm = _("     norm of the residual = ")
86     msgCT   = _("     computing time       = ")
87     printf(msgNorm + string(norm_res1) + "\n" + msgCT+string(t1)+"\n")
88     disp(_("A2 ("+descr2+"): order = "+ string(n2) + " nnz = " + string(nnz(A2))))
89     printf(msgNorm + string(norm_res2) + "\n" + msgCT + string(t2) + "\n")
90     disp(_("A3 ("+descr3+"): order = "+ string(n3) + " nnz = " + string(nnz(A3))))
91     printf(msgNorm + string(norm_res3) + "\n" + msgCT + string(t3)+"\n")
92     disp(_("A4 ("+descr4+"): order = "+ string(n4) + " nnz = " + string(nnz(A4))));
93     printf(msgNorm + string(norm_res4)+"\n" + msgCT + string(t4)+"\n\n")
94
95     disp(_("Now we will see how to use the lu factors:"));
96     disp(_("  1/ lu factors of a sparse matrix A are obtained through:"));
97     printf(_("     lup = umf_lufact(A)\n     lup is a pointer to the lu factors (the memory is outside scilab space)\n"))
98     disp(_("  2/ for solving a linear system A*x = b, use:"));
99     printf(_("     x = umf_lusolve(lup,b)\n     or  x = umf_lusolve(lup,b,""Ax=b"",A)  to add an iterative refinement step\n"))
100     disp(_("  3/ to solve A''*x=b you may use:"));
101     printf(_("     x = umf_lusolve(lup,b,""Ax''''=b"")\n     or  x = umf_lusolve(lup,b,""Ax''''=b"",A)  to add an iterative refinement step\n"))
102     disp(_("  4/ you may also compute the 1-norm condition number quicky with:"));
103     printf(_("     K1 = condestsp(A,lup)\n     K1 = condestsp(A) also works but in this case the lu factors are re-computed inside\n"))
104     disp(_("  5/ if you don''t need the lu factors anymore, it is recommended to free them with:"));
105     printf(_("     umf_ludel(lup)\n     if you have lost your pointer you may use umf_ludel() which frees all the current umf lu factors\n"))
106     halt(_("\nPress Return to continue...\n"));
107
108     mode(1)
109     lup1 = umf_lufact(A1);
110     x1 = umf_lusolve(lup1,b1);
111     norm_res1 =  norm(A1*x1-b1);
112     x1 = umf_lusolve(lup1,b1,"Ax=b",A1);
113     norm_res1r = norm(A1*x1-b1);
114     K1 = condestsp(A1,lup1);
115     umf_ludel(lup1);
116
117     lup2 = umf_lufact(A2);
118     x2 = umf_lusolve(lup2,b2);
119     norm_res2 = norm(A2*x2-b2);
120     x2 = umf_lusolve(lup2,b2,"Ax=b",A2);
121     norm_res2r = norm(A2*x2-b2);
122     K2 = condestsp(A2,lup2);
123     umf_ludel(lup2);
124
125     lup3 = umf_lufact(A3);
126     x3 = umf_lusolve(lup3,b3);
127     norm_res3 = norm(A3*x3-b3);
128     x3 = umf_lusolve(lup3,b3,"Ax=b",A3);
129     norm_res3r = norm(A3*x3-b3);
130     K3 = condestsp(A3,lup3);
131     umf_ludel(lup3);
132
133     lup4 = umf_lufact(A4);
134     x4 = umf_lusolve(lup4,b4);
135     norm_res4 = norm(A4*x4-b4);
136     x4 = umf_lusolve(lup4,b4,"Ax=b",A4);
137     norm_res4r = norm(A4*x4-b4);
138     K4 = condestsp(A4,lup4,5);
139     umf_ludel(lup4);
140     mode(-1)
141
142     disp("A1 ("+descr1+"): order = "+ string(n1) + " nnz = " + string(nnz(A1)))
143     printf("     K1 = " + string(K1)+"\n" + ..
144          _("     norm of the residual    = ") + string(norm_res1) + "\n" + ..
145          _("     same but with refinment = ") + string(norm_res1r) + "\n")
146     disp("A2 ("+descr2+"): order = "+ string(n2) + " nnz = " + string(nnz(A2)))
147     printf("     K2 = " + string(K2)+"\n"+..
148          _("     norm of the residual    = ") + string(norm_res2) + "\n" + ..
149          _("     same but with refinment = ") + string(norm_res2r) + "\n")
150     disp("A3 ("+descr3+"): order = "+ string(n3) + " nnz = " + string(nnz(A3)))
151     printf("     K3 = " + string(K3)+"\n"+..
152          _("     norm of the residual    = ") + string(norm_res3) + "\n" + ..
153          _("     same but with refinment = ") + string(norm_res3r) + "\n")
154     disp("A4 ("+descr4+"): order = "+ string(n4) + " nnz = " + string(nnz(A4)))
155     printf("     K4 = " + string(K4) + "\n" + ..
156          _("     norm of the residual     = ") + string(norm_res4) + "\n" + ..
157          _("     same but with refinement = ") + string(norm_res4r) + "\n\n")
158
159     clear A1 A3 A4 x1 x3 x4 b1 b3 b4
160
161     disp(_("Now we will see how to use the taucs snmf stuff on the matrix A2."));
162     disp(_("This is useful and recommended when your matrix is symmetric positive definite (s.p.d.)."));
163     disp(_("  1/ the Cholesky factorization of a s.p.d. matrix A is obtained with:"));
164     printf("     Cp = taucs_chfact(A)\n" + ..
165          _("     Cp is a pointer to the Cholesky fact. (the memory is outside scilab space)\n"));
166     disp(_("  2/ for solving a linear system A*x = b then use:"));
167     printf(_("     x = taucs_chsolve(Cp,b)\n"));
168     disp(_("  3/ for the same thing with one refinement step, use:"));
169     printf(_("     xr = taucs_chsolve(Cp,b,A)\n"));
170     disp(_("  4/ you may also compute the 2-norm condition number with:"));
171     printf(_("     [K2, lm, vm, lM, vM] = cond2sp(A, Cp [, itermax, eps, verb])\n"));
172     disp(_("  5/ if you don''t need the Cholesky factorization anymore, it is recommended to free it with taucs_chdel(Cp)"));
173     printf(_("     if you have lost your pointer you may use taucs_chdel() which frees memory for all the current Cholesky factorizations.\n"));
174     halt(_("\nPress Return to continue...\n"));
175
176     mode(1)
177     tic();
178     Cp2 = taucs_chfact(A2);
179     x2 = taucs_chsolve(Cp2,b2);
180     t2_chol = toc();
181     norm_res_chol_2 = norm(A2*x2-b2);
182     x2r = taucs_chsolve(Cp2,b2,A2);
183     norm_res_chol_2r = norm(A2*x2r-b2);
184     K2_norm2 = cond2sp(A2, Cp2);
185     taucs_chdel(Cp2);
186     mode(-1)
187
188     disp(_("A2 ("+descr2+"): order = "+ string(n2) + " nnz = " + string(nnz(A2))));
189     printf("     K2 (1-norm) = " + string(K2) + "\n" + ..
190            "     K2 (2-norm) = " + string(K2_norm2) + "\n" )
191     disp(_(" with umfpack:"));
192     printf(_("     norm of the residual = ") + string(norm_res2) + "\n" + ..
193            _("     same but with refinment = ") + string(norm_res2r) + "\n" + ..
194            _("     computing time       = ") + string(t2) + "\n")
195     disp(_(" with the taucs snmf Cholesky solver:"))
196     printf(_("     norm of the residual = ") + string(norm_res_chol_2) + "\n" + ..
197            _("     same but with refinment = ") + string(norm_res_chol_2r) + "\n" + ..
198            _("     computing time       = ") + string(t2_chol) + "\n")
199
200     if num == max(winsid()) then
201         xdel(num);
202     end
203
204 endfunction
205
206 demo_sparse_matrices()
207 clear demo_sparse_matrices;