Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_scale.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === KLU_scale ============================================================ */
00003 /* ========================================================================== */
00004 
00005 /* Scale a matrix and check to see if it is valid.  Can be called by the user.
00006  * This is called by KLU_factor and KLU_refactor.  Returns TRUE if the input
00007  * matrix is valid, FALSE otherwise.  If the W input argument is non-NULL,
00008  * then the input matrix is checked for duplicate entries.
00009  *
00010  * scaling methods:
00011  *  <0: no scaling, do not compute Rs, and do not check input matrix.
00012  *  0: no scaling
00013  *  1: the scale factor for row i is sum (abs (A (i,:)))
00014  *  2 or more: the scale factor for row i is max (abs (A (i,:)))
00015  */
00016 
00017 #include "amesos_klu_internal.h"
00018 
00019 Int KLU_scale   /* return TRUE if successful, FALSE otherwise */
00020 (
00021     /* inputs, not modified */
00022     Int scale,    /* 0: none, 1: sum, 2: max */
00023     Int n,
00024     Int Ap [ ],   /* size n+1, column pointers */
00025     Int Ai [ ],   /* size nz, row indices */
00026     double Ax [ ],
00027     /* outputs, not defined on input */
00028     double Rs [ ],  /* size n, can be NULL if scale <= 0 */
00029     /* workspace, not defined on input or output */
00030     Int W [ ],    /* size n, can be NULL */
00031     /* --------------- */
00032     KLU_common *Common
00033 )
00034 {
00035     double a ;
00036     Entry *Az ;
00037     Int row, col, p, pend, check_duplicates ;
00038 
00039     /* ---------------------------------------------------------------------- */
00040     /* check inputs */
00041     /* ---------------------------------------------------------------------- */
00042 
00043     if (Common == NULL)
00044     {
00045   return (FALSE) ;
00046     }
00047     Common->status = KLU_OK ;
00048 
00049     if (scale < 0)
00050     {
00051   /* return without checking anything and without computing the
00052    * scale factors */
00053   return (TRUE) ;
00054     }
00055 
00056     Az = (Entry *) Ax ;
00057 
00058     if (n <= 0 || Ap == NULL || Ai == NULL || Az == NULL ||
00059   (scale > 0 && Rs == NULL))
00060     {
00061   /* Ap, Ai, Ax and Rs must be present, and n must be > 0 */
00062   Common->status = KLU_INVALID ;
00063   return (FALSE) ;
00064     }
00065     if (Ap [0] != 0 || Ap [n] < 0)
00066     {
00067   /* nz = Ap [n] must be >= 0 and Ap [0] must equal zero */
00068   Common->status = KLU_INVALID ;
00069   return (FALSE) ;
00070     }
00071     for (col = 0 ; col < n ; col++)
00072     {
00073   if (Ap [col] > Ap [col+1])
00074   {
00075       /* column pointers must be non-decreasing */
00076       Common->status = KLU_INVALID ;
00077       return (FALSE) ;
00078   }
00079     }
00080 
00081     /* ---------------------------------------------------------------------- */
00082     /* scale */
00083     /* ---------------------------------------------------------------------- */
00084 
00085     if (scale > 0)
00086     {
00087   /* initialize row sum or row max */
00088   for (row = 0 ; row < n ; row++)
00089   {
00090       Rs [row] = 0 ;
00091   }
00092     }
00093 
00094     /* check for duplicates only if W is present */
00095     check_duplicates = (W != (Int *) NULL) ;
00096     if (check_duplicates)
00097     {
00098   for (row = 0 ; row < n ; row++)
00099   {
00100       W [row] = EMPTY ;
00101   }
00102     }
00103 
00104     for (col = 0 ; col < n ; col++)
00105     {
00106   pend = Ap [col+1] ;
00107   for (p = Ap [col] ; p < pend ; p++)
00108   {
00109       row = Ai [p] ;
00110       if (row < 0 || row >= n)
00111       {
00112     /* row index out of range, or duplicate entry */
00113     Common->status = KLU_INVALID ;
00114     return (FALSE) ;
00115       }
00116       if (check_duplicates)
00117       {
00118     if (W [row] == col)
00119     {
00120         /* duplicate entry */
00121         Common->status = KLU_INVALID ;
00122         return (FALSE) ;
00123     }
00124     /* flag row i as appearing in column col */
00125     W [row] = col ;
00126       }
00127       /* a = ABS (Az [p]) ;*/
00128       ABS (a, Az [p]) ;
00129       if (scale == 1)
00130       {
00131     /* accumulate the abs. row sum */
00132     Rs [row] += a ;
00133       }
00134       else if (scale > 1)
00135       {
00136     /* find the max abs. value in the row */
00137     Rs [row] = MAX (Rs [row], a) ;
00138       }
00139   }
00140     }
00141 
00142     if (scale > 0)
00143     {
00144   /* do not scale empty rows */
00145   for (row = 0 ; row < n ; row++)
00146   {
00147       /* matrix is singular */
00148       PRINTF (("Rs [%d] = %g\n", row, Rs [row])) ;
00149 
00150       if (Rs [row] == 0.0)
00151       {
00152     PRINTF (("Row %d of A is all zero\n", row)) ;
00153     Rs [row] = 1.0 ;
00154       }
00155   }
00156     }
00157 
00158     return (TRUE) ;
00159 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines