← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/sparse/SuperLU/SRC/sgsequ.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/SuperLU/SRC/sgsequ.c
 
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:   sgsequ.c
 
12
 * History:     Modified from LAPACK routine SGEEQU
 
13
 */
 
14
#include <math.h>
 
15
#include "ssp_defs.h"
 
16
#include "util.h"
 
17
 
 
18
void
 
19
sgsequ(SuperMatrix *A, float *r, float *c, float *rowcnd,
 
20
        float *colcnd, float *amax, int *info)
 
21
{
 
22
/*    
 
23
    Purpose   
 
24
    =======   
 
25
 
 
26
    SGSEQU computes row and column scalings intended to equilibrate an   
 
27
    M-by-N sparse matrix A and reduce its condition number. R returns the row
 
28
    scale factors and C the column scale factors, chosen to try to make   
 
29
    the largest element in each row and column of the matrix B with   
 
30
    elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.   
 
31
 
 
32
    R(i) and C(j) are restricted to be between SMLNUM = smallest safe   
 
33
    number and BIGNUM = largest safe number.  Use of these scaling   
 
34
    factors is not guaranteed to reduce the condition number of A but   
 
35
    works well in practice.   
 
36
 
 
37
    See supermatrix.h for the definition of 'SuperMatrix' structure.
 
38
 
 
39
    Arguments   
 
40
    =========   
 
41
 
 
42
    A       (input) SuperMatrix*
 
43
            The matrix of dimension (A->nrow, A->ncol) whose equilibration
 
44
            factors are to be computed. The type of A can be:
 
45
            Stype = SLU_NC; Dtype = SLU_S; Mtype = SLU_GE.
 
46
            
 
47
    R       (output) float*, size A->nrow
 
48
            If INFO = 0 or INFO > M, R contains the row scale factors   
 
49
            for A.
 
50
            
 
51
    C       (output) float*, size A->ncol
 
52
            If INFO = 0,  C contains the column scale factors for A.
 
53
            
 
54
    ROWCND  (output) float*
 
55
            If INFO = 0 or INFO > M, ROWCND contains the ratio of the   
 
56
            smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and   
 
57
            AMAX is neither too large nor too small, it is not worth   
 
58
            scaling by R.
 
59
            
 
60
    COLCND  (output) float*
 
61
            If INFO = 0, COLCND contains the ratio of the smallest   
 
62
            C(i) to the largest C(i).  If COLCND >= 0.1, it is not   
 
63
            worth scaling by C.
 
64
            
 
65
    AMAX    (output) float*
 
66
            Absolute value of largest matrix element.  If AMAX is very   
 
67
            close to overflow or very close to underflow, the matrix   
 
68
            should be scaled.
 
69
            
 
70
    INFO    (output) int*
 
71
            = 0:  successful exit   
 
72
            < 0:  if INFO = -i, the i-th argument had an illegal value   
 
73
            > 0:  if INFO = i,  and i is   
 
74
                  <= A->nrow:  the i-th row of A is exactly zero   
 
75
                  >  A->ncol:  the (i-M)-th column of A is exactly zero   
 
76
 
 
77
    ===================================================================== 
 
78
*/
 
79
 
 
80
    /* Local variables */
 
81
    NCformat *Astore;
 
82
    float   *Aval;
 
83
    int i, j, irow;
 
84
    float rcmin, rcmax;
 
85
    float bignum, smlnum;
 
86
    extern double slamch_(char *);
 
87
    
 
88
    /* Test the input parameters. */
 
89
    *info = 0;
 
90
    if ( A->nrow < 0 || A->ncol < 0 ||
 
91
         A->Stype != SLU_NC || A->Dtype != SLU_S || A->Mtype != SLU_GE )
 
92
        *info = -1;
 
93
    if (*info != 0) {
 
94
        i = -(*info);
 
95
        xerbla_("sgsequ", &i);
 
96
        return;
 
97
    }
 
98
 
 
99
    /* Quick return if possible */
 
100
    if ( A->nrow == 0 || A->ncol == 0 ) {
 
101
        *rowcnd = 1.;
 
102
        *colcnd = 1.;
 
103
        *amax = 0.;
 
104
        return;
 
105
    }
 
106
 
 
107
    Astore = A->Store;
 
108
    Aval = Astore->nzval;
 
109
    
 
110
    /* Get machine constants. */
 
111
    smlnum = slamch_("S");
 
112
    bignum = 1. / smlnum;
 
113
 
 
114
    /* Compute row scale factors. */
 
115
    for (i = 0; i < A->nrow; ++i) r[i] = 0.;
 
116
 
 
117
    /* Find the maximum element in each row. */
 
118
    for (j = 0; j < A->ncol; ++j)
 
119
        for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) {
 
120
            irow = Astore->rowind[i];
 
121
            r[irow] = SUPERLU_MAX( r[irow], fabs(Aval[i]) );
 
122
        }
 
123
 
 
124
    /* Find the maximum and minimum scale factors. */
 
125
    rcmin = bignum;
 
126
    rcmax = 0.;
 
127
    for (i = 0; i < A->nrow; ++i) {
 
128
        rcmax = SUPERLU_MAX(rcmax, r[i]);
 
129
        rcmin = SUPERLU_MIN(rcmin, r[i]);
 
130
    }
 
131
    *amax = rcmax;
 
132
 
 
133
    if (rcmin == 0.) {
 
134
        /* Find the first zero scale factor and return an error code. */
 
135
        for (i = 0; i < A->nrow; ++i)
 
136
            if (r[i] == 0.) {
 
137
                *info = i + 1;
 
138
                return;
 
139
            }
 
140
    } else {
 
141
        /* Invert the scale factors. */
 
142
        for (i = 0; i < A->nrow; ++i)
 
143
            r[i] = 1. / SUPERLU_MIN( SUPERLU_MAX( r[i], smlnum ), bignum );
 
144
        /* Compute ROWCND = min(R(I)) / max(R(I)) */
 
145
        *rowcnd = SUPERLU_MAX( rcmin, smlnum ) / SUPERLU_MIN( rcmax, bignum );
 
146
    }
 
147
 
 
148
    /* Compute column scale factors */
 
149
    for (j = 0; j < A->ncol; ++j) c[j] = 0.;
 
150
 
 
151
    /* Find the maximum element in each column, assuming the row
 
152
       scalings computed above. */
 
153
    for (j = 0; j < A->ncol; ++j)
 
154
        for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) {
 
155
            irow = Astore->rowind[i];
 
156
            c[j] = SUPERLU_MAX( c[j], fabs(Aval[i]) * r[irow] );
 
157
        }
 
158
 
 
159
    /* Find the maximum and minimum scale factors. */
 
160
    rcmin = bignum;
 
161
    rcmax = 0.;
 
162
    for (j = 0; j < A->ncol; ++j) {
 
163
        rcmax = SUPERLU_MAX(rcmax, c[j]);
 
164
        rcmin = SUPERLU_MIN(rcmin, c[j]);
 
165
    }
 
166
 
 
167
    if (rcmin == 0.) {
 
168
        /* Find the first zero scale factor and return an error code. */
 
169
        for (j = 0; j < A->ncol; ++j)
 
170
            if ( c[j] == 0. ) {
 
171
                *info = A->nrow + j + 1;
 
172
                return;
 
173
            }
 
174
    } else {
 
175
        /* Invert the scale factors. */
 
176
        for (j = 0; j < A->ncol; ++j)
 
177
            c[j] = 1. / SUPERLU_MIN( SUPERLU_MAX( c[j], smlnum ), bignum);
 
178
        /* Compute COLCND = min(C(J)) / max(C(J)) */
 
179
        *colcnd = SUPERLU_MAX( rcmin, smlnum ) / SUPERLU_MIN( rcmax, bignum );
 
180
    }
 
181
 
 
182
    return;
 
183
 
 
184
} /* sgsequ */
 
185
 
 
186