1 
2 /*! @file slaqgs.c
3  * \brief Equlibrates a general sprase matrix
4  *
5  * <pre>
6  * -- SuperLU routine (version 2.0) --
7  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
8  * and Lawrence Berkeley National Lab.
9  * November 15, 1997
10  *
11  * Modified from LAPACK routine SLAQGE
12  * </pre>
13  */
14 /*
15  * File name:	slaqgs.c
16  * History:     Modified from LAPACK routine SLAQGE
17  */
18 #include <math.h>
19 #include "slu_sdefs.h"
20 
21 /*! \brief
22  *
23  * <pre>
24  *   Purpose
25  *   =======
26  *
27  *   SLAQGS equilibrates a general sparse M by N matrix A using the row and
28  *   scaling factors in the vectors R and C.
29  *
30  *   See supermatrix.h for the definition of 'SuperMatrix' structure.
31  *
32  *   Arguments
33  *   =========
34  *
35  *   A       (input/output) SuperMatrix*
36  *           On exit, the equilibrated matrix.  See EQUED for the form of
37  *           the equilibrated matrix. The type of A can be:
38  *	    Stype = NC; Dtype = SLU_S; Mtype = GE.
39  *
40  *   R       (input) float*, dimension (A->nrow)
41  *           The row scale factors for A.
42  *
43  *   C       (input) float*, dimension (A->ncol)
44  *           The column scale factors for A.
45  *
46  *   ROWCND  (input) float
47  *           Ratio of the smallest R(i) to the largest R(i).
48  *
49  *   COLCND  (input) float
50  *           Ratio of the smallest C(i) to the largest C(i).
51  *
52  *   AMAX    (input) float
53  *           Absolute value of largest matrix entry.
54  *
55  *   EQUED   (output) char*
56  *           Specifies the form of equilibration that was done.
57  *           = 'N':  No equilibration
58  *           = 'R':  Row equilibration, i.e., A has been premultiplied by
59  *                   diag(R).
60  *           = 'C':  Column equilibration, i.e., A has been postmultiplied
61  *                   by diag(C).
62  *           = 'B':  Both row and column equilibration, i.e., A has been
63  *                   replaced by diag(R) * A * diag(C).
64  *
65  *   Internal Parameters
66  *   ===================
67  *
68  *   THRESH is a threshold value used to decide if row or column scaling
69  *   should be done based on the ratio of the row or column scaling
70  *   factors.  If ROWCND < THRESH, row scaling is done, and if
71  *   COLCND < THRESH, column scaling is done.
72  *
73  *   LARGE and SMALL are threshold values used to decide if row scaling
74  *   should be done based on the absolute size of the largest matrix
75  *   element.  If AMAX > LARGE or AMAX < SMALL, row scaling is done.
76  *
77  *   =====================================================================
78  * </pre>
79  */
80 
81 void
slaqgs(SuperMatrix * A,float * r,float * c,float rowcnd,float colcnd,float amax,char * equed)82 slaqgs(SuperMatrix *A, float *r, float *c,
83 	float rowcnd, float colcnd, float amax, char *equed)
84 {
85 
86 
87 #define THRESH    (0.1)
88 
89     /* Local variables */
90     NCformat *Astore;
91     float   *Aval;
92     int i, j, irow;
93     float large, small, cj;
94 
95 
96     /* Quick return if possible */
97     if (A->nrow <= 0 || A->ncol <= 0) {
98 	*(unsigned char *)equed = 'N';
99 	return;
100     }
101 
102     Astore = A->Store;
103     Aval = Astore->nzval;
104 
105     /* Initialize LARGE and SMALL. */
106     small = slamch_("Safe minimum") / slamch_("Precision");
107     large = 1. / small;
108 
109     if (rowcnd >= THRESH && amax >= small && amax <= large) {
110 	if (colcnd >= THRESH)
111 	    *(unsigned char *)equed = 'N';
112 	else {
113 	    /* Column scaling */
114 	    for (j = 0; j < A->ncol; ++j) {
115 		cj = c[j];
116 		for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) {
117 		    Aval[i] *= cj;
118                 }
119 	    }
120 	    *(unsigned char *)equed = 'C';
121 	}
122     } else if (colcnd >= THRESH) {
123 	/* Row scaling, no column scaling */
124 	for (j = 0; j < A->ncol; ++j)
125 	    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) {
126 		irow = Astore->rowind[i];
127 		Aval[i] *= r[irow];
128 	    }
129 	*(unsigned char *)equed = 'R';
130     } else {
131 	/* Row and column scaling */
132 	for (j = 0; j < A->ncol; ++j) {
133 	    cj = c[j];
134 	    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) {
135 		irow = Astore->rowind[i];
136 		Aval[i] *= cj * r[irow];
137 	    }
138 	}
139 	*(unsigned char *)equed = 'B';
140     }
141 
142     return;
143 
144 } /* slaqgs */
145 
146