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