← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/linsolve/SuperLU/SRC/spivotgrowth.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:
scipy/linsolve/SuperLU/SRC/spivotgrowth.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
#include <math.h>
 
11
#include "ssp_defs.h"
 
12
#include "util.h"
 
13
 
 
14
float
 
15
sPivotGrowth(int ncols, SuperMatrix *A, int *perm_c, 
 
16
             SuperMatrix *L, SuperMatrix *U)
 
17
{
 
18
/*
 
19
 * Purpose
 
20
 * =======
 
21
 *
 
22
 * Compute the reciprocal pivot growth factor of the leading ncols columns
 
23
 * of the matrix, using the formula:
 
24
 *     min_j ( max_i(abs(A_ij)) / max_i(abs(U_ij)) )
 
25
 *
 
26
 * Arguments
 
27
 * =========
 
28
 *
 
29
 * ncols    (input) int
 
30
 *          The number of columns of matrices A, L and U.
 
31
 *
 
32
 * A        (input) SuperMatrix*
 
33
 *          Original matrix A, permuted by columns, of dimension
 
34
 *          (A->nrow, A->ncol). The type of A can be:
 
35
 *          Stype = NC; Dtype = SLU_S; Mtype = GE.
 
36
 *
 
37
 * L        (output) SuperMatrix*
 
38
 *          The factor L from the factorization Pr*A=L*U; use compressed row 
 
39
 *          subscripts storage for supernodes, i.e., L has type: 
 
40
 *          Stype = SC; Dtype = SLU_S; Mtype = TRLU.
 
41
 *
 
42
 * U        (output) SuperMatrix*
 
43
 *          The factor U from the factorization Pr*A*Pc=L*U. Use column-wise
 
44
 *          storage scheme, i.e., U has types: Stype = NC;
 
45
 *          Dtype = SLU_S; Mtype = TRU.
 
46
 *
 
47
 */
 
48
    NCformat *Astore;
 
49
    SCformat *Lstore;
 
50
    NCformat *Ustore;
 
51
    float  *Aval, *Lval, *Uval;
 
52
    int      fsupc, nsupr, luptr, nz_in_U;
 
53
    int      i, j, k, oldcol;
 
54
    int      *inv_perm_c;
 
55
    float   rpg, maxaj, maxuj;
 
56
    extern   double slamch_(char *);
 
57
    float   smlnum;
 
58
    float   *luval;
 
59
   
 
60
    /* Get machine constants. */
 
61
    smlnum = slamch_("S");
 
62
    rpg = 1. / smlnum;
 
63
 
 
64
    Astore = A->Store;
 
65
    Lstore = L->Store;
 
66
    Ustore = U->Store;
 
67
    Aval = Astore->nzval;
 
68
    Lval = Lstore->nzval;
 
69
    Uval = Ustore->nzval;
 
70
    
 
71
    inv_perm_c = (int *) SUPERLU_MALLOC(A->ncol*sizeof(int));
 
72
    for (j = 0; j < A->ncol; ++j) inv_perm_c[perm_c[j]] = j;
 
73
 
 
74
    for (k = 0; k <= Lstore->nsuper; ++k) {
 
75
        fsupc = L_FST_SUPC(k);
 
76
        nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);
 
77
        luptr = L_NZ_START(fsupc);
 
78
        luval = &Lval[luptr];
 
79
        nz_in_U = 1;
 
80
        
 
81
        for (j = fsupc; j < L_FST_SUPC(k+1) && j < ncols; ++j) {
 
82
            maxaj = 0.;
 
83
            oldcol = inv_perm_c[j];
 
84
            for (i = Astore->colptr[oldcol]; i < Astore->colptr[oldcol+1]; ++i)
 
85
                maxaj = SUPERLU_MAX( maxaj, fabs(Aval[i]) );
 
86
        
 
87
            maxuj = 0.;
 
88
            for (i = Ustore->colptr[j]; i < Ustore->colptr[j+1]; i++)
 
89
                maxuj = SUPERLU_MAX( maxuj, fabs(Uval[i]) );
 
90
            
 
91
            /* Supernode */
 
92
            for (i = 0; i < nz_in_U; ++i)
 
93
                maxuj = SUPERLU_MAX( maxuj, fabs(luval[i]) );
 
94
 
 
95
            ++nz_in_U;
 
96
            luval += nsupr;
 
97
 
 
98
            if ( maxuj == 0. )
 
99
                rpg = SUPERLU_MIN( rpg, 1.);
 
100
            else
 
101
                rpg = SUPERLU_MIN( rpg, maxaj / maxuj );
 
102
        }
 
103
        
 
104
        if ( j >= ncols ) break;
 
105
    }
 
106
 
 
107
    SUPERLU_FREE(inv_perm_c);
 
108
    return (rpg);
 
109
}