/* ========================================================================== */
/* === klu_scale ============================================================ */
/* ========================================================================== */
/* Scale a matrix and check to see if it is valid. Can be called by the user.
* This is called by klu_factor and klu_refactor. Returns TRUE if the input
* matrix is valid, FALSE otherwise. If the W input argument is non-NULL,
* then the input matrix is checked for duplicate entries.
*
* scaling methods:
* <0: no scaling, do not compute Rs, and do not check input matrix.
* 0: no scaling
* 1: the scale factor for row i is sum (abs (A (i,:)))
* 2 or more: the scale factor for row i is max (abs (A (i,:)))
*/
#include "klu_internal.h"
int KLU_scale /* return TRUE if successful, FALSE otherwise */
(
/* inputs, not modified */
int scale, /* 0: none, 1: sum, 2: max */
int n,
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
double Ax [ ],
/* outputs, not defined on input */
double Rs [ ], /* size n, can be NULL if scale <= 0 */
/* workspace, not defined on input or output */
int W [ ], /* size n, can be NULL */
/* --------------- */
klu_common *Common
)
{
double a ;
Entry *Az ;
int row, col, p, pend, check_duplicates ;
if (Common == NULL)
{
return (FALSE) ;
}
Common->status = KLU_OK ;
if (scale < 0)
{
/* return without checking anything and without computing the
* scale factors */
return (TRUE) ;
}
Az = (Entry *) Ax ;
if (n <= 0 || (Ap == NULL) || (Ai == NULL) || (Az == NULL))
{
/* Ap, Ai, and Ax must be present, and n must be > 0 */
Common->status = KLU_INVALID ;
return (FALSE) ;
}
if (Ap [0] != 0 || Ap [n] < 0)
{
/* nz = Ap [n] must be >= 0 and Ap [0] must equal zero */
Common->status = KLU_INVALID ;
return (FALSE) ;
}
for (col = 0 ; col < n ; col++)
{
if (Ap [col] > Ap [col+1])
{
/* column pointers must be non-decreasing */
Common->status = KLU_INVALID ;
return (FALSE) ;
}
}
if (scale > 0)
{
/* initialize row sum or row max */
for (row = 0 ; row < n ; row++)
{
Rs [row] = 0 ;
}
}
/* check for duplicates only if W is present */
check_duplicates = (W != (int *) NULL) ;
if (check_duplicates)
{
for (row = 0 ; row < n ; row++)
{
W [row] = EMPTY ;
}
}
for (col = 0 ; col < n ; col++)
{
pend = Ap [col+1] ;
for (p = Ap [col] ; p < pend ; p++)
{
row = Ai [p] ;
if (row < 0 || row >= n)
{
/* row index out of range, or duplicate entry */
Common->status = KLU_INVALID ;
return (FALSE) ;
}
if (check_duplicates)
{
if (W [row] == col)
{
/* duplicate entry */
Common->status = KLU_INVALID ;
return (FALSE) ;
}
/* flag row i as appearing in column col */
W [row] = col ;
}
/* a = ABS (Az [p]) ;*/
ABS (a, Az [p]) ;
if (scale == 1)
{
/* accumulate the abs. row sum */
Rs [row] += a ;
}
else if (scale > 1)
{
/* find the max abs. value in the row */
Rs [row] = MAX (Rs [row], a) ;
}
}
}
if (scale > 0)
{
/* do not scale empty rows */
for (row = 0 ; row < n ; row++)
{
/* matrix is singular */
PRINTF (("Rs [%d] = %g\n", row, Rs [row])) ;
if (Rs [row] == 0.0)
{
PRINTF (("Row %d of A is all zero\n", row)) ;
Rs [row] = 1.0 ;
}
}
}
return (TRUE) ;
}
syntax highlighted by Code2HTML, v. 0.9.1