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