b85427419b66f128dc4c5e00f4c5557303c8ba55
[scilab.git] / scilab / modules / umfpack / examples / scisptdem1.dem
1 // the simple demo for the umfpack interface & others
2 // sparse utilility (scispt toolbox)
3
4 old_winsid = winsid();
5
6 num1=max(winsid())+1
7 num2=num1+1
8 num3=num1+2
9 num4=num1+3
10
11 x_message(["To load a sparse matrix  stored  in  the  Harwell";
12            "Boeing format in scilab, you may use the function";
13            "                                                 ";
14            "                 ReadHBSparse                    ";
15            "                                                 ";
16            "4 matrices will be loaded for the demo : they co-";
17            "mes from the University of Florida Sparse Matrix "
18            "Collection :                                     ";
19            "                                                 ";
20            "    www.cise.ufl.edu/research/sparse/matrices    ";
21            "                                                 ";    
22            "maintained by Tim Davis the author of UMFPACK    ";
23            "                                                 ";
24            "       CLICK ON OK TO CONTINUE                   "]);
25
26 mode(1)
27 [A1,descr1,ref1,mtype1] = ReadHBSparse(SCI+"/modules/umfpack/examples/arc130.rua");
28 [A2,descr2,ref2,mtype2] = ReadHBSparse(SCI+"/modules/umfpack/examples/bcsstk24.rsa");
29 [A3,descr3,ref3,mtype3] = ReadHBSparse(SCI+"/modules/umfpack/examples/ex14.rua");
30 [A4,descr4,ref4,mtype4] = ReadHBSparse(SCI+"/modules/umfpack/examples/young1c.csa");
31
32 mode(-1)
33 x_message(["To see the pattern of non zeros elements";
34            "you may use the function :              ";
35            "                                        ";
36            "         PlotSparse                     ";
37            "                                        ";
38            "  here we will display the 4 matrices   ";
39            "                                        ";
40            "       CLICK ON OK TO CONTINUE          "]);
41
42 mode(1)
43
44 scf(num1)
45 PlotSparse(A1,"y+")
46 xtitle(ref1+"."+mtype1+" : "+descr1)
47
48 scf(num2)
49 PlotSparse(A2)
50 xtitle(ref2+"."+mtype2+" : "+descr2)
51
52 scf(num3)
53 PlotSparse(A3,"c.")
54 xtitle(ref3+"."+mtype3+" : "+descr3)
55
56 scf(num4)
57 PlotSparse(A4,"r.")
58 xtitle(ref4+"."+mtype4+" : "+descr4)
59
60 mode(-1)
61
62 x_message(["Now use umfpack to solve some linear system :      ";
63            "                                                   ";
64            "b1 = rand(size(A1,1),1);    -> to create a rhs     ";
65            "x1 = umfpack(A1,""\"",b1);    -> to solve A1x1 = b1  ";
66            "norm_res = norm(A1*x1-b1);  -> norm of the residual";
67            "                                                   ";
68            "this is done for the 4 matrices A1, A2, A3 and A4  ";
69            "                                                   ";
70            "      CLICK ON OK TO CONTINUE                      "]);
71 mode(1)
72 n1 = size(A1,1);
73 b1 = rand(n1,1);           
74 timer(); x1 = umfpack(A1,"\",b1); t1=timer();
75 norm_res1 = norm(A1*x1-b1);
76
77 n2 = size(A2,1);
78 b2 = rand(n2,1);           
79 timer(); x2 = umfpack(A2,"\",b2); t2=timer();
80 norm_res2 = norm(A2*x2-b2);
81
82 n3 = size(A3,1);
83 b3 = rand(n3,1);           
84 timer(); x3 = umfpack(A3,"\",b3); t3=timer();
85 norm_res3 =  norm(A3*x3-b3);
86
87 n4 = size(A4,1);
88 b4 = rand(n4,1);           
89 timer(); x4 = umfpack(A4,"\",b4); t4=timer();
90 norm_res4 = norm(A4*x4-b4);
91
92 mode(-1)
93
94 x_message([" A1 ("+descr1+") : order = "+ string(n1) + " nnz = " + string(nnz(A1));
95            " norm of the residual = "+string(norm_res1) ;
96            " computing time       = "+string(t1)  ;
97            "                                     ";
98            " A2 ("+descr2+") : order = "+ string(n2) + " nnz = " + string(nnz(A2));
99            " norm of the residual = "+string(norm_res2) ;
100            " computing time       = "+string(t2)  ;
101            "                                     ";
102            " A3 ("+descr3+") : order = "+ string(n3) + " nnz = " + string(nnz(A3));
103            " norm of the residual = "+string(norm_res3) ;
104            " computing time       = "+string(t3)  ;
105            "                                     ";
106            " A4 ("+descr4+") : order = "+ string(n4) + " nnz = " + string(nnz(A4));
107            " norm of the residual = "+string(norm_res4) ;
108            " computing time       = "+string(t4)  ;
109            "                                     ";
110            "      CLICK ON OK TO CONTINUE        "]);
111
112 x_message(["             Now we will see how to use the lu factors                 ";
113            "                                                                       ";
114            " 1/ lu factors of a sparse matrix A are gotten with :                  ";
115            "                                                                       ";
116            "    lup = umf_lufact(A)                                                ";
117            "                                                                       ";
118            " lup is a pointer to the lu factors (the  memory  is outside scilab space)  ";
119            "                                                                       ";
120            " 2/ for solving a linear system A x = b use :                          ";
121            "                                                                       ";
122            "    x = umf_lusolve(lup,b)                                             ";
123            " or x = umf_lusolve(lup,b,""Ax=b"",A)  to add an iterative refinement step";
124            "                                                                       ";
125            " 3/ to solve A''x=b you may use :                                      ";
126            "                                                                       ";                    
127            "    x = umf_lusolve(lup,b,""Ax''''=b"")                                ";
128            " or x = umf_lusolve(lup,b,""Ax''''=b"",A) to add an iterative refinement step";
129            "                                                                       ";
130            " 4/ you may also compute the 1-norm condition number quicky with :     ";
131            "                                                                       ";
132            "    K1 = condestsp(A,lup)                                              ";
133            "                                                                       ";
134            " K1 = condestsp(A) works also but in this case the lu factors are re-computed";
135            " inside.                                                               ";
136            "                                                                       ";
137            " 5/  WHEN you does''t need anymore the lu factors it is recommended to ";
138            "     free the memory used by them with :                               ";
139            "                                                                       ";
140            "        umf_ludel(lup)                                                 ";
141            "                                                                       ";
142            " If you have lost your pointer you may use  umf_ludel() which free all the ";
143            " current umf lu factors.                                               ";
144            "                                                                       ";
145            "                       CLICK ON OK TO CONTINUE                         "]);
146            
147 mode(1)
148
149 lup1 = umf_lufact(A1);
150 x1 = umf_lusolve(lup1,b1);
151 norm_res1 =  norm(A1*x1-b1);
152 x1 = umf_lusolve(lup1,b1,"Ax=b",A1);
153 norm_res1r = norm(A1*x1-b1);
154 K1 = condestsp(A1,lup1);
155 umf_ludel(lup1);
156
157 lup2 = umf_lufact(A2);
158 x2 = umf_lusolve(lup2,b2);
159 norm_res2 = norm(A2*x2-b2);
160 x2 = umf_lusolve(lup2,b2,"Ax=b",A2);
161 norm_res2r = norm(A2*x2-b2);
162 K2 = condestsp(A2,lup2);
163 umf_ludel(lup2);
164
165 lup3 = umf_lufact(A3);
166 x3 = umf_lusolve(lup3,b3);
167 norm_res3 = norm(A3*x3-b3);
168 x3 = umf_lusolve(lup3,b3,"Ax=b",A3);
169 norm_res3r = norm(A3*x3-b3);
170 K3 = condestsp(A3,lup3);
171 umf_ludel(lup3);
172
173 lup4 = umf_lufact(A4);
174 x4 = umf_lusolve(lup4,b4);
175 norm_res4 = norm(A4*x4-b4);
176 x4 = umf_lusolve(lup4,b4,"Ax=b",A4);
177 norm_res4r = norm(A4*x4-b4);
178 K4 = condestsp(A4,lup4,5);
179 umf_ludel(lup4);
180
181 mode(-1)
182
183 x_message([" A1 ("+descr1+") : order = "+ string(n1) + " nnz = " + string(nnz(A1));
184            " K1 = " + string(K1);
185            " norm of the residual      = "+string(norm_res1) ;
186            " idem but with raffinement = "+string(norm_res1r);
187            "                                     ";
188            " A2 ("+descr2+") : order = "+ string(n2) + " nnz = " + string(nnz(A2));
189            " K2 = " + string(K2);
190            " norm of the residual      = "+string(norm_res2) ;
191            " idem but with raffinement = "+string(norm_res2r);
192            "                                     ";
193            " A3 ("+descr3+") : order = "+ string(n3) + " nnz = " + string(nnz(A3));
194            " K3 = " + string(K3);
195            " norm of the residual      = "+string(norm_res3) ;
196            " idem but with raffinement = "+string(norm_res3r);
197            "                                     ";
198            " A4 ("+descr4+") : order = "+ string(n4) + " nnz = " + string(nnz(A4));
199            " K4 = " + string(K4);
200            " norm of the residual      = "+string(norm_res4) ;
201            " idem but with raffinement = "+string(norm_res4r);
202            "                                     ";
203            "      CLICK ON OK TO CONTINUE        "]);
204
205
206 clear A1 A3 A4 x1 x3 x4 b1 b3 b4
207
208 x_message(["            Now we will see how to use the taucs snmf stuff            ";
209            " This is useful and recommended when your matrix is symetric positive  ";
210            " definite (s.p.d.)                                                     ";
211            "                                                                       ";
212            " 1/ The Cholesky factorization of a s.p.d. matrix A is gotten with :   ";
213            "                                                                       ";
214            "    Cp = taucs_chfact(A)                                               ";
215            "                                                                       ";
216            " Cp is a pointer to the Cholesky fact. (the  memory  is outside scilab ";
217            " space)                                                                ";
218            "                                                                       ";
219            " 2/ for solving a linear system A x = b then use :                     ";
220            "                                                                       ";
221            "    x = taucs_chsolve(Cp,b)                                            ";
222            "                                                                       ";
223            " 3/ for the same thing with one refinement step use :                  ";
224            "                                                                       ";
225            "    xr = taucs_chsolve(Cp,b,A)                                         ";
226            "                                                                       ";
227            " 4/ you may also compute the 2-norm condition number with :            ";
228            "                                                                       ";
229            "    [K2, lm, vm, lM, vM] = cond2sp(A, Cp [, itermax, eps, verb])       ";
230            "                                                                       ";
231            "                                                                       ";
232            " 5/  WHEN you does''t need the Cholesky factorization is recommended to";
233            "     free the memory used by them with :                               ";
234            "                                                                       ";
235            "        taucs_chdel(Cp)                                                ";
236            "                                                                       ";
237            " If you have lost your pointer you may use  taucs_chdel() which free   ";
238            " memory for all the current Cholesky factorizations.                   ";
239            "                                                                       ";
240            "                       CLICK ON OK TO CONTINUE                         "]);
241
242 mode(1)
243 timer();
244 Cp2 = taucs_chfact(A2);
245 x2 = taucs_chsolve(Cp2,b2);
246 t2_chol = timer();
247 norm_res_chol_2 = norm(A2*x2-b2);
248 [x2r, r2] = rafiter(A2, Cp2, b2, x2);
249 norm_res_chol_2r = norm(r2);
250 K2_norm2 = cond2sp(A2, Cp2);
251 taucs_chdel(Cp2);
252 mode(-1)
253
254 x_message([" On the s.p.d. matrix A2 : ";
255            " A2 ("+descr2+") : order = "+ string(n2) + " nnz = " + string(nnz(A2));
256            " K2 (1-norm) = " + string(K2);
257            " K2 (2-norm) = " + string(K2_norm2);
258            "   ";
259            " with umfpack : ";
260            " norm of the residual      = "+string(norm_res2) ;
261            " idem but with raffinement = "+string(norm_res2r);
262            " computing time            = "+string(t2)  ;
263            "                                     ";
264            " with the taucs snmf Cholesky solver ";
265            " norm of the residual      = "+string(norm_res_chol_2) ;
266            " idem but with raffinement = "+string(norm_res_chol_2r);
267            " computing time            = "+string(t2_chol)  ;
268            "                                     ";
269            "      CLICK ON OK TO CONTINUE        "]);
270
271
272 x_message(["   Demo is out : I hope that this stuff will be  ";
273            " useful for you ! Don''t hesitate to mail me about";
274            " bugs (Bruno.Pincon@iecn.u-nancy.fr) or simply to";
275            " tell that this interface is of interest for you ";
276            "                                                 ";
277            "              CLICK ON OK TO FINISH              "]);
278           
279 xdel(num1); xdel(num2); xdel(num3) ; xdel(num4)
280 if old_winsid == [] then, xdel(0), end
281
282
283