src/main/java/edu/ufl/cise/klu/tdouble/Dklu_scale.java

   1 /**
   2  * KLU: a sparse LU factorization algorithm.
   3  * Copyright (C) 2004-2009, Timothy A. Davis.
   4  * Copyright (C) 2011, Richard W. Lincoln.
   5  * http://www.cise.ufl.edu/research/sparse/klu
   6  *
   7  * -------------------------------------------------------------------------
   8  *
   9  * KLU is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public
  11  * License as published by the Free Software Foundation; either
  12  * version 2.1 of the License, or (at your option) any later version.
  13  *
  14  * KLU is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with this Module; if not, write to the Free Software
  21  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  22  *
  23  */
  24
  25 package edu.ufl.cise.klu.tdouble;
  26
  27 /**
  28  * Scale a matrix and check to see if it is valid.  Can be called by the user.
  29  * This is called by KLU_factor and KLU_refactor.  Returns true if the input
  30  * matrix is valid, false otherwise.  If the W input argument is non-null,
  31  * then the input matrix is checked for duplicate entries.
  32  *
  33  * scaling methods:
  34  *      <0: no scaling, do not compute Rs, and do not check input matrix.
  35  *      0: no scaling
  36  *      1: the scale factor for row i is sum (abs (A (i,:)))
  37  *      2 or more: the scale factor for row i is max(abs (A (i,:)))
  38  */
  39 public class Dklu_scale {
  40
  41         /**
  42          *
  43          * @param scale 0: none, 1: sum, 2: max
  44          * @param n
  45          * @param Ap size n+1, column pointers
  46          * @param Ai size nz, row indices
  47          * @param Ax outputs, not defined on input
  48          * @param Rs size n, can be null if scale <= 0
  49          * @param W size n, can be null
  50          * @param Common
  51          * @return true if successful, false otherwise
  52          */
  53         public static boolean klu_scale(int scale, int n, int[] Ap, int[] Ai,
  54                     double[] Ax, double[] Rs, int[] W, KLU_common Common)
  55         {
  56                 double a;
  57             Entry[] Az;
  58             int row, col, p, pend;
  59             boolean check_duplicates;
  60
  61             /* ------------------------------------------------------------------ */
  62             /* check inputs */
  63             /* ------------------------------------------------------------------ */
  64
  65             if (Common == null)
  66             {
  67                 return false;
  68             }
  69             Common.status = KLU_OK;
  70
  71             if (scale < 0)
  72             {
  73                 /* return without checking anything and without computing the
  74                  * scale factors */
  75                 return true;
  76             }
  77
  78             Az = (Entry[]) Ax;
  79
  80             if (n <= 0 || Ap == null || Ai == null || Az == null ||
  81                 (scale > 0 && Rs == null))
  82             {
  83                 /* Ap, Ai, Ax and Rs must be present, and n must be > 0 */
  84                 Common.status = KLU_INVALID;
  85                 return false;
  86             }
  87             if (Ap[0] != 0 || Ap[n] < 0)
  88             {
  89                 /* nz = Ap[n] must be >= 0 and Ap[0] must equal zero */
  90                 Common.status = KLU_INVALID;
  91                 return false;
  92             }
  93             for (col = 0; col < n; col++)
  94             {
  95                 if (Ap[col] > Ap[col+1])
  96                 {
  97                     /* column pointers must be non-decreasing */
  98                     Common.status = KLU_INVALID;
  99                     return false;
 100                 }
 101             }
 102
 103             /* ------------------------------------------------------------------ */
 104             /* scale */
 105             /* ------------------------------------------------------------------ */
 106
 107             if (scale > 0)
 108             {
 109                 /* initialize row sum or row max*/
 110                 for (row = 0; row < n; row++)
 111                 {
 112                     Rs[row] = 0;
 113                 }
 114             }
 115
 116             /* check for duplicates only if W is present */
 117             check_duplicates = (W != null);
 118             if (check_duplicates)
 119             {
 120                 for (row = 0; row < n; row++)
 121                 {
 122                     W[row] = EMPTY;
 123                 }
 124             }
 125
 126             for (col = 0; col < n; col++)
 127             {
 128                 pend = Ap[col+1];
 129                 for (p = Ap[col]; p < pend; p++)
 130                 {
 131                     row = Ai[p];
 132                     if (row < 0 || row >= n)
 133                     {
 134                         /* row index out of range, or duplicate entry */
 135                         Common.status = KLU_INVALID;
 136                         return false;
 137                     }
 138                     if (check_duplicates)
 139                     {
 140                         if (W[row] == col)
 141                         {
 142                             /* duplicate entry */
 143                             Common.status = KLU_INVALID;
 144                             return false;
 145                         }
 146                         /* flag row i as appearing in column col */
 147                         W[row] = col;
 148                     }
 149                     /* a = ABS (Az[p]);*/
 150                     ABS (a, Az[p]);
 151                     if (scale == 1)
 152                     {
 153                         /* accumulate the abs. row sum */
 154                         Rs[row] += a;
 155                     }
 156                     else if (scale > 1)
 157                     {
 158                         /* find the maxabs. value in the row */
 159                         Rs[row] = max(Rs[row], a);
 160                     }
 161                 }
 162             }
 163
 164             if (scale > 0)
 165             {
 166                 /* do not scale empty rows */
 167                 for (row = 0; row < n; row++)
 168                 {
 169                     /* matrix is singular */
 170                     printf("Rs[%d] = %g\n", row, Rs[row]);
 171
 172                     if (Rs[row] == 0.0)
 173                     {
 174                         printf("Row %d of A is all zero\n", row);
 175                         Rs[row] = 1.0;
 176                     }
 177                 }
 178             }
 179
 180             return true;
 181         }
 182
 183 }