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