← 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/zpivotgrowth.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/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
}