← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/superlu/slaqgs.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/sandbox/pysparse/superlu/slaqgs.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:   slaqgs.c
12
 
 * History:     Modified from LAPACK routine SLAQGE
13
 
 */
14
 
#include <math.h>
15
 
#include "ssp_defs.h"
16
 
#include "util.h"
17
 
 
18
 
void
19
 
slaqgs(SuperMatrix *A, float *r, float *c, 
20
 
        float rowcnd, float colcnd, float amax, char *equed)
21
 
{
22
 
/*
23
 
    Purpose   
24
 
    =======   
25
 
 
26
 
    SLAQGS 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 = S_; Mtype = GE.
38
 
            
39
 
    R       (input) float*, dimension (A->nrow)
40
 
            The row scale factors for A.
41
 
            
42
 
    C       (input) float*, dimension (A->ncol)
43
 
            The column scale factors for A.
44
 
            
45
 
    ROWCND  (input) float
46
 
            Ratio of the smallest R(i) to the largest R(i).
47
 
            
48
 
    COLCND  (input) float
49
 
            Ratio of the smallest C(i) to the largest C(i).
50
 
            
51
 
    AMAX    (input) float
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
 
    float   *Aval;
84
 
    int i, j, irow;
85
 
    float large, small, cj;
86
 
    extern double slamch_(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 = slamch_("Safe minimum") / slamch_("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
 
} /* slaqgs */
138