← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/linsolve/SuperLU/SRC/zgsequ.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/linsolve/SuperLU/SRC/zgsequ.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:   zgsequ.c
12
 
 * History:     Modified from LAPACK routine ZGEEQU
13
 
 */
14
 
#include <math.h>
15
 
#include "zsp_defs.h"
16
 
#include "util.h"
17
 
 
18
 
void
19
 
zgsequ(SuperMatrix *A, double *r, double *c, double *rowcnd,
20
 
        double *colcnd, double *amax, int *info)
21
 
{
22
 
/*    
23
 
    Purpose   
24
 
    =======   
25
 
 
26
 
    ZGSEQU 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_Z; Mtype = SLU_GE.
46
 
            
47
 
    R       (output) double*, size A->nrow
48
 
            If INFO = 0 or INFO > M, R contains the row scale factors   
49
 
            for A.
50
 
            
51
 
    C       (output) double*, size A->ncol
52
 
            If INFO = 0,  C contains the column scale factors for A.
53
 
            
54
 
    ROWCND  (output) double*
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) double*
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) double*
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
 
    doublecomplex   *Aval;
83
 
    int i, j, irow;
84
 
    double rcmin, rcmax;
85
 
    double bignum, smlnum;
86
 
    extern double dlamch_(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_Z || A->Mtype != SLU_GE )
92
 
        *info = -1;
93
 
    if (*info != 0) {
94
 
        i = -(*info);
95
 
        xerbla_("zgsequ", &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 = dlamch_("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], z_abs1(&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], z_abs1(&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
 
} /* zgsequ */
185
 
 
186