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