Mercurial > forge
diff main/sparse/SuperLU/SRC/dlaqgs.c @ 0:6b33357c7561 octave-forge
Initial revision
author | pkienzle |
---|---|
date | Wed, 10 Oct 2001 19:54:49 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/sparse/SuperLU/SRC/dlaqgs.c Wed Oct 10 19:54:49 2001 +0000 @@ -0,0 +1,138 @@ + + +/* + * -- SuperLU routine (version 2.0) -- + * Univ. of California Berkeley, Xerox Palo Alto Research Center, + * and Lawrence Berkeley National Lab. + * November 15, 1997 + * + */ +/* + * File name: dlaqgs.c + * History: Modified from LAPACK routine DLAQGE + */ +#include <math.h> +#include "dsp_defs.h" +#include "util.h" + +void +dlaqgs(SuperMatrix *A, double *r, double *c, + double rowcnd, double colcnd, double amax, char *equed) +{ +/* + Purpose + ======= + + DLAQGS equilibrates a general sparse M by N matrix A using the row and + scaling factors in the vectors R and C. + + See supermatrix.h for the definition of 'SuperMatrix' structure. + + Arguments + ========= + + A (input/output) SuperMatrix* + On exit, the equilibrated matrix. See EQUED for the form of + the equilibrated matrix. The type of A can be: + Stype = NC; Dtype = _D; Mtype = GE. + + R (input) double*, dimension (A->nrow) + The row scale factors for A. + + C (input) double*, dimension (A->ncol) + The column scale factors for A. + + ROWCND (input) double + Ratio of the smallest R(i) to the largest R(i). + + COLCND (input) double + Ratio of the smallest C(i) to the largest C(i). + + AMAX (input) double + Absolute value of largest matrix entry. + + EQUED (output) char* + Specifies the form of equilibration that was done. + = 'N': No equilibration + = 'R': Row equilibration, i.e., A has been premultiplied by + diag(R). + = 'C': Column equilibration, i.e., A has been postmultiplied + by diag(C). + = 'B': Both row and column equilibration, i.e., A has been + replaced by diag(R) * A * diag(C). + + Internal Parameters + =================== + + THRESH is a threshold value used to decide if row or column scaling + should be done based on the ratio of the row or column scaling + factors. If ROWCND < THRESH, row scaling is done, and if + COLCND < THRESH, column scaling is done. + + LARGE and SMALL are threshold values used to decide if row scaling + should be done based on the absolute size of the largest matrix + element. If AMAX > LARGE or AMAX < SMALL, row scaling is done. + + ===================================================================== +*/ + +#define THRESH (0.1) + + /* Local variables */ + NCformat *Astore; + double *Aval; + int i, j, irow; + double large, small, cj; + extern double dlamch_(char *); + + + /* Quick return if possible */ + if (A->nrow <= 0 || A->ncol <= 0) { + *(unsigned char *)equed = 'N'; + return; + } + + Astore = A->Store; + Aval = Astore->nzval; + + /* Initialize LARGE and SMALL. */ + small = dlamch_("Safe minimum") / dlamch_("Precision"); + large = 1. / small; + + if (rowcnd >= THRESH && amax >= small && amax <= large) { + if (colcnd >= THRESH) + *(unsigned char *)equed = 'N'; + else { + /* Column scaling */ + for (j = 0; j < A->ncol; ++j) { + cj = c[j]; + for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) { + Aval[i] *= cj; + } + } + *(unsigned char *)equed = 'C'; + } + } else if (colcnd >= THRESH) { + /* Row scaling, no column scaling */ + for (j = 0; j < A->ncol; ++j) + for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) { + irow = Astore->rowind[i]; + Aval[i] *= r[irow]; + } + *(unsigned char *)equed = 'R'; + } else { + /* Row and column scaling */ + for (j = 0; j < A->ncol; ++j) { + cj = c[j]; + for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) { + irow = Astore->rowind[i]; + Aval[i] *= cj * r[irow]; + } + } + *(unsigned char *)equed = 'B'; + } + + return; + +} /* dlaqgs */ +