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