← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/SuperLU/SRC/zpivotgrowth.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 "zsp_defs.h"
12
 
#include "util.h"
13
 
 
14
 
double
15
 
zPivotGrowth(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_Z; 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_Z; 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_Z; Mtype = TRU.
46
 
 *
47
 
 */
48
 
    NCformat *Astore;
49
 
    SCformat *Lstore;
50
 
    NCformat *Ustore;
51
 
    doublecomplex  *Aval, *Lval, *Uval;
52
 
    int      fsupc, nsupr, luptr, nz_in_U;
53
 
    int      i, j, k, oldcol;
54
 
    int      *inv_perm_c;
55
 
    double   rpg, maxaj, maxuj;
56
 
    extern   double dlamch_(char *);
57
 
    double   smlnum;
58
 
    doublecomplex   *luval;
59
 
    doublecomplex   temp_comp;
60
 
   
61
 
    /* Get machine constants. */
62
 
    smlnum = dlamch_("S");
63
 
    rpg = 1. / smlnum;
64
 
 
65
 
    Astore = A->Store;
66
 
    Lstore = L->Store;
67
 
    Ustore = U->Store;
68
 
    Aval = Astore->nzval;
69
 
    Lval = Lstore->nzval;
70
 
    Uval = Ustore->nzval;
71
 
    
72
 
    inv_perm_c = (int *) SUPERLU_MALLOC(A->ncol*sizeof(int));
73
 
    for (j = 0; j < A->ncol; ++j) inv_perm_c[perm_c[j]] = j;
74
 
 
75
 
    for (k = 0; k <= Lstore->nsuper; ++k) {
76
 
        fsupc = L_FST_SUPC(k);
77
 
        nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);
78
 
        luptr = L_NZ_START(fsupc);
79
 
        luval = &Lval[luptr];
80
 
        nz_in_U = 1;
81
 
        
82
 
        for (j = fsupc; j < L_FST_SUPC(k+1) && j < ncols; ++j) {
83
 
            maxaj = 0.;
84
 
            oldcol = inv_perm_c[j];
85
 
            for (i = Astore->colptr[oldcol]; i < Astore->colptr[oldcol+1]; ++i)
86
 
                maxaj = SUPERLU_MAX( maxaj, z_abs1( &Aval[i]) );
87
 
        
88
 
            maxuj = 0.;
89
 
            for (i = Ustore->colptr[j]; i < Ustore->colptr[j+1]; i++)
90
 
                maxuj = SUPERLU_MAX( maxuj, z_abs1( &Uval[i]) );
91
 
            
92
 
            /* Supernode */
93
 
            for (i = 0; i < nz_in_U; ++i)
94
 
                maxuj = SUPERLU_MAX( maxuj, z_abs1( &luval[i]) );
95
 
 
96
 
            ++nz_in_U;
97
 
            luval += nsupr;
98
 
 
99
 
            if ( maxuj == 0. )
100
 
                rpg = SUPERLU_MIN( rpg, 1.);
101
 
            else
102
 
                rpg = SUPERLU_MIN( rpg, maxaj / maxuj );
103
 
        }
104
 
        
105
 
        if ( j >= ncols ) break;
106
 
    }
107
 
 
108
 
    SUPERLU_FREE(inv_perm_c);
109
 
    return (rpg);
110
 
}