← Back to branch summary

~ubuntu-branches/ubuntu/karmic/hypre/karmic

« back to all changes in this revision

Viewing changes to src/FEI_mv/SuperLU/SRC/cpivotgrowth.c

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2009-03-20 11:40:12 UTC
  • mfrom: (4.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20090320114012-132h6ok9w2r6o609
Tags: 2.4.0b-2
Rebuild against new openmpi.

Show diffs side-by-side

added added

removed removed

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