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;