src/test/java/edu/ufl/cise/klu/tdouble/demo/Dklu_demo.java

   1 package edu.ufl.cise.klu.tdouble.demo;
   2
   3 import java.io.FileNotFoundException;
   4 import java.io.FileReader;
   5 import java.io.IOException;
   6
   7 import edu.ufl.cise.klu.common.KLU_common;
   8 import edu.ufl.cise.klu.common.KLU_numeric;
   9 import edu.ufl.cise.klu.common.KLU_symbolic;
  10 import edu.ufl.cise.klu.common.KLU_version;
  11 import static edu.ufl.cise.klu.tdouble.Dklu_analyze.klu_analyze;
  12 import static edu.ufl.cise.klu.tdouble.Dklu_defaults.klu_defaults;
  13 import static edu.ufl.cise.klu.tdouble.Dklu_diagnostics.klu_rgrowth;
  14 import static edu.ufl.cise.klu.tdouble.Dklu_diagnostics.klu_condest;
  15 import static edu.ufl.cise.klu.tdouble.Dklu_diagnostics.klu_rcond;
  16 import static edu.ufl.cise.klu.tdouble.Dklu_diagnostics.klu_flops;
  17 import static edu.ufl.cise.klu.tdouble.Dklu_factor.klu_factor;
  18 import static edu.ufl.cise.klu.tdouble.Dklu_solve.klu_solve;
  19 import edu.ufl.cise.klu.tdouble.io.MatrixInfo;
  20 import edu.ufl.cise.klu.tdouble.io.MatrixSize;
  21 import edu.ufl.cise.klu.tdouble.io.MatrixVectorReader;
  22
  23 /**
  24  * Read in a Matrix Market matrix(using CHOLMOD) and solve a linear system.
  25  */
  26 public class Dklu_demo {
  27
  28         private static void REAL (double[] X, int i, double v)
  29         {
  30                 X[2*i] = v;
  31         }
  32
  33         private static void IMAG (double[] X, int i, double v)
  34         {
  35                 X[2*i + 1] = v;
  36         }
  37
  38 //      private static double REAL (double[] X, int i)
  39 //      {
  40 //              return X[2*i] ;
  41 //      }
  42 //
  43 //      private static double IMAG (double[] X, int i)
  44 //      {
  45 //              return X[2*i + 1] ;
  46 //      }
  47 //
  48 //      private static double CABS (double[] X, int i)
  49 //      {
  50 //              return Math.sqrt(REAL(X, i) * REAL(X, i) + IMAG(X, i) * IMAG(X, i)) ;
  51 //      }
  52
  53         private static double MAX (double a, double b)
  54         {
  55                 return (a) > (b) ? (a) : (b) ;
  56         }
  57
  58         /**
  59          *
  60          * @param n A is n-by-n
  61          * @param Ap size n+1, column pointers
  62          * @param Ai size nz = Ap [n], row indices
  63          * @param Ax size nz, numerical values
  64          * @param isreal nonzero if A is real, 0 otherwise
  65          * @param B size n, right-hand-side
  66          * @param X size n, solution to Ax=b
  67          * @param R size n, residual r = b-A*x
  68          * @param lunz size 1, nnz(L+U+F)
  69          * @param rnorm size 1, norm(b-A*x,1) / norm(A,1)
  70          * @param Common default parameters and statistics
  71          * @return 1 if successful, 0 otherwise
  72          */
  73         public static int klu_backslash(int n, int[] Ap, int[] Ai, double[] Ax,
  74                         int isreal, double[] B, double[] X, double[] R, int[] lunz,
  75                         double[] rnorm, KLU_common Common)
  76         {
  77                 double anorm = 0, asum;
  78                 KLU_symbolic Symbolic;
  79                 KLU_numeric Numeric;
  80                 int i, j, p;
  81
  82                 if (Ap == null || Ai == null || Ax == null || B == null || X == null || B == null)
  83                         return(0);
  84
  85                 /* ---------------------------------------------------------------------- */
  86                 /* symbolic ordering and analysis */
  87                 /* ---------------------------------------------------------------------- */
  88
  89                 Symbolic = klu_analyze(n, Ap, Ai, Common);
  90                 if (Symbolic == null) return(0);
  91
  92                 if (isreal != 0)
  93                 {
  94
  95                         /* ------------------------------------------------------------------ */
  96                         /* factorization */
  97                         /* ------------------------------------------------------------------ */
  98
  99                         Numeric = klu_factor(Ap, Ai, Ax, Symbolic, Common);
 100                         if (Numeric == null)
 101                         {
 102                                 //Dklu_free_symbolic.klu_free_symbolic(Symbolic, Common);
 103                                 return(0);
 104                         }
 105
 106                         /* ------------------------------------------------------------------ */
 107                         /* statistics(not required to solve Ax=b) */
 108                         /* ------------------------------------------------------------------ */
 109
 110                         klu_rgrowth(Ap, Ai, Ax, Symbolic, Numeric, Common);
 111                         klu_condest(Ap, Ax, Symbolic, Numeric, Common);
 112                         klu_rcond(Symbolic, Numeric, Common);
 113                         klu_flops(Symbolic, Numeric, Common);
 114                         lunz[0] = Numeric.lnz + Numeric.unz - n +
 115                                 (Numeric.Offp != null ? Numeric.Offp [n] : 0);
 116
 117                         /* ------------------------------------------------------------------ */
 118                         /* solve Ax=b */
 119                         /* ------------------------------------------------------------------ */
 120
 121                         for(i = 0; i < n; i++)
 122                         {
 123                                 X [i] = B [i];
 124                         }
 125                         klu_solve(Symbolic, Numeric, n, 1, X, 0, Common);
 126
 127                         /* ------------------------------------------------------------------ */
 128                         /* compute residual, rnorm = norm(b-Ax,1) / norm(A,1) */
 129                         /* ------------------------------------------------------------------ */
 130
 131                         for(i = 0; i < n; i++)
 132                         {
 133                                 R [i] = B [i];
 134                         }
 135                         for(j = 0; j < n; j++)
 136                         {
 137                                 asum = 0;
 138                                 for(p = Ap [j]; p < Ap [j+1]; p++)
 139                                 {
 140                                         /* R(i) -= A(i,j) * X(j) */
 141                                         R [Ai [p]] -= Ax [p] * X [j];
 142                                         asum += Math.abs(Ax [p]);
 143                                 }
 144                                 anorm = MAX (anorm, asum);
 145                         }
 146                         rnorm[0] = 0;
 147                         for(i = 0; i < n; i++)
 148                         {
 149                                 rnorm[0] = MAX (rnorm[0], Math.abs(R [i]));
 150                         }
 151
 152                         /* ------------------------------------------------------------------ */
 153                         /* free numeric factorization */
 154                         /* ------------------------------------------------------------------ */
 155
 156                         //klu_free_numeric(Numeric, Common);
 157                         Numeric = null;
 158
 159                 }
 160                 else
 161                 {
 162                         throw new UnsupportedOperationException();
 163
 164                         /* ------------------------------------------------------------------ */
 165                         /* statistics(not required to solve Ax=b) */
 166                         /* ------------------------------------------------------------------ */
 167
 168 //                      Numeric = klu_z_factor(Ap, Ai, Ax, Symbolic, Common);
 169 //                      if (Numeric == null)
 170 //                      {
 171 //                              klu_free_symbolic(Symbolic, Common);
 172 //                              return(0);
 173 //                      }
 174 //
 175 //                      /* ------------------------------------------------------------------ */
 176 //                      /* statistics */
 177 //                      /* ------------------------------------------------------------------ */
 178 //
 179 //                      klu_z_rgrowth(Ap, Ai, Ax, Symbolic, Numeric, Common);
 180 //                      klu_z_condest(Ap, Ax, Symbolic, Numeric, Common);
 181 //                      klu_z_rcond(Symbolic, Numeric, Common);
 182 //                      klu_z_flops(Symbolic, Numeric, Common);
 183 //                      lunz = Numeric.lnz + Numeric.unz - n +
 184 //                              (Numeric.Offp != null ? Numeric.Offp [n] : 0);
 185 //
 186 //                      /* ------------------------------------------------------------------ */
 187 //                      /* solve Ax=b */
 188 //                      /* ------------------------------------------------------------------ */
 189 //
 190 //                      for(i = 0; i < 2*n; i++)
 191 //                      {
 192 //                              X [i] = B [i];
 193 //                      }
 194 //                      klu_z_solve(Symbolic, Numeric, n, 1, X, Common);
 195 //
 196 //                      /* ------------------------------------------------------------------ */
 197 //                      /* compute residual, rnorm = norm(b-Ax,1) / norm(A,1) */
 198 //                      /* ------------------------------------------------------------------ */
 199 //
 200 //                      for(i = 0; i < 2*n; i++)
 201 //                      {
 202 //                              R [i] = B [i];
 203 //                      }
 204 //                      for(j = 0; j < n; j++)
 205 //                      {
 206 //                              asum = 0;
 207 //                              for(p = Ap [j]; p < Ap [j+1]; p++)
 208 //                              {
 209 //                                      /* R(i) -= A(i,j) * X(j) */
 210 //                                      i = Ai [p];
 211 //                                      REAL(R,i) -= REAL(Ax,p) * REAL(X,j) - IMAG(Ax,p) * IMAG(X,j);
 212 //                                      IMAG(R,i) -= IMAG(Ax,p) * REAL(X,j) + REAL(Ax,p) * IMAG(X,j);
 213 //                                      asum += CABS(Ax, p);
 214 //                              }
 215 //                              anorm = MAX(anorm, asum);
 216 //                      }
 217 //                      rnorm = 0;
 218 //                      for(i = 0; i < n; i++)
 219 //                      {
 220 //                              rnorm = MAX (rnorm, CABS(R, i));
 221 //                      }
 222 //
 223 //                      /* ------------------------------------------------------------------ */
 224 //                      /* free numeric factorization */
 225 //                      /* ------------------------------------------------------------------ */
 226 //
 227 //                      klu_z_free_numeric(&Numeric, Common);
 228                 }
 229
 230                 /* ---------------------------------------------------------------------- */
 231                 /* free symbolic analysis, and residual */
 232                 /* ---------------------------------------------------------------------- */
 233
 234                 //klu_free_symbolic(Symbolic, Common);
 235                 Symbolic = null;
 236
 237                 return (1);
 238         }
 239
 240         /**
 241          * Given a sparse matrix A, set up a right-hand-side and solve X = A\b.
 242          */
 243         public static void klu_demo(int n, int[] Ap, int[] Ai, double[] Ax,
 244                         int isreal)
 245         {
 246                 KLU_common Common = new KLU_common();
 247                 double[] B, X, R;
 248                 int i;
 249                 int[] lunz = new int[1];
 250                 double[] rnorm = new double[1];
 251
 252                 System.out.printf("KLU: %s, version: %d.%d.%d\n",
 253                                 KLU_version.KLU_DATE, KLU_version.KLU_MAIN_VERSION,
 254                                 KLU_version.KLU_SUB_VERSION, KLU_version.KLU_SUBSUB_VERSION);
 255
 256                 /* ---------------------------------------------------------------------- */
 257                 /* set defaults */
 258                 /* ---------------------------------------------------------------------- */
 259
 260                 klu_defaults(Common);
 261
 262                 /* ---------------------------------------------------------------------- */
 263                 /* create a right-hand-side */
 264                 /* ---------------------------------------------------------------------- */
 265
 266                 if (isreal != 0)
 267                 {
 268                         /* B = 1 +(1:n)/n */
 269                         B = new double[n];
 270                         X = new double[n];
 271                         R = new double[n];
 272                         if (B != null)
 273                         {
 274                                 for(i = 0; i < n; i++)
 275                                 {
 276                                         B [i] = 1 + ((double) i+1) / ((double) n);
 277                                 }
 278                         }
 279                 }
 280                 else
 281                 {
 282                         /* real(B) = 1 +(1:n)/n, imag(B) = (n:-1:1)/n */
 283                         B = new double[2 * n];
 284                         X = new double[2 * n];
 285                         R = new double[2 * n];
 286                         if (B != null)
 287                         {
 288                                 for(i = 0; i < n; i++)
 289                                 {
 290                                         REAL(B, i, 1 + ((double) i+1) / ((double) n));
 291                                         IMAG(B, i,     ((double) n-i) / ((double) n));
 292                                 }
 293                         }
 294                 }
 295
 296                 /* ---------------------------------------------------------------------- */
 297                 /* X = A\b using KLU and print statistics */
 298                 /* ---------------------------------------------------------------------- */
 299
 300                 if (klu_backslash(n, Ap, Ai, Ax, isreal, B, X, R, lunz, rnorm, Common) == 0)
 301                 {
 302                         System.out.printf("KLU failed\n");
 303                 }
 304                 else
 305                 {
 306                         System.out.printf("n %d nnz(A) %d nnz(L+U+F) %d resid %g\n" +
 307                                 "recip growth %g condest %g rcond %g flops %g\n",
 308                                 n, Ap [n], lunz[0], rnorm, Common.rgrowth, Common.condest,
 309                                 Common.rcond, Common.flops);
 310                 }
 311
 312                 /* ---------------------------------------------------------------------- */
 313                 /* free the problem */
 314                 /* ---------------------------------------------------------------------- */
 315
 316                 if (isreal != 0)
 317                 {
 318                         B = null;
 319                         X = null;
 320                         R = null;
 321                 }
 322                 else
 323                 {
 324                         B = null;
 325                         X = null;
 326                         R = null;
 327                 }
 328 //              System.out.printf("peak memory usage: %g bytes\n\n",(double)(Common.mempeak));
 329         }
 330
 331         /**
 332          * Read in a sparse matrix in Matrix Market format using CHOLMOD, and then
 333          * solve Ax=b with KLU.  Note that CHOLMOD is only used to read the matrix.
 334          */
 335         public static void main(String[] args) {
 336
 337                 FileReader fileReader;
 338                 MatrixVectorReader reader;
 339                 MatrixInfo info;
 340                 MatrixSize size;
 341
 342                 try {
 343                         fileReader = new FileReader(args[0]);
 344                         reader = new MatrixVectorReader(fileReader);
 345
 346                         info = reader.readMatrixInfo();
 347                         size = reader.readMatrixSize(info);
 348
 349                         if (size.numRows() != size.numColumns() || !info.isGeneral()
 350                                 ||(!(info.isReal() || info.isComplex()))
 351                                 || !info.isCoordinate())
 352                         {
 353                                 System.out.printf("invalid matrix\n");
 354                         }
 355                         else
 356                         {
 357                                 int nnz = size.numEntries();
 358                                 int[] i = new int[nnz];
 359                                 int[] p = new int[nnz];
 360                                 double[] x;
 361
 362                                 if (info.isComplex()) {
 363                                         double[] xR = new double[nnz];
 364                                         double[] xI = new double[nnz];
 365
 366                                         reader.readCoordinate(i, p, xR, xI);
 367
 368                                         x = new double[2 * nnz];
 369                                         for (int j = 0; j < nnz; j++) {
 370                                                 REAL(x, j, xR[j]);
 371                                                 IMAG(x, j, xI[j]);
 372                                         }
 373                                 }
 374                                 else
 375                                 {
 376                                         x = new double[nnz];
 377                                         reader.readCoordinate(i, p, x);
 378                                 }
 379
 380
 381                                 klu_demo(size.numRows(), p, i, x, info.isReal() ? 1 : 0);
 382                         }
 383
 384                 } catch (FileNotFoundException e) {
 385                         e.printStackTrace();
 386                 } catch (IOException e) {
 387                         e.printStackTrace();
 388                 }
 389         }
 390
 391 }