← 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/cpivotL.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/cpivotL.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
 
  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
12
 
 
13
 
  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
14
 
  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
15
 
 
16
 
  Permission is hereby granted to use or copy this program for any
17
 
  purpose, provided the above notices are retained on all copies.
18
 
  Permission to modify the code and to distribute modified code is
19
 
  granted, provided the above notices are retained, and a notice that
20
 
  the code was modified is included with the above copyright notice.
21
 
*/
22
 
 
23
 
#include <math.h>
24
 
#include <stdlib.h>
25
 
#include "csp_defs.h"
26
 
#include "util.h"
27
 
 
28
 
#undef DEBUG
29
 
 
30
 
int
31
 
cpivotL(
32
 
        const int  jcol,     /* in */
33
 
        const float u,      /* in - diagonal pivoting threshold */
34
 
        int        *usepr,   /* re-use the pivot sequence given by perm_r/iperm_r */
35
 
        int        *perm_r,  /* may be modified */
36
 
        int        *iperm_r, /* in - inverse of perm_r */
37
 
        int        *iperm_c, /* in - used to find diagonal of Pc*A*Pc' */
38
 
        int        *pivrow,  /* out */
39
 
        GlobalLU_t *Glu      /* modified - global LU data structures */
40
 
       )
41
 
{
42
 
/*
43
 
 * Purpose
44
 
 * =======
45
 
 *   Performs the numerical pivoting on the current column of L,
46
 
 *   and the CDIV operation.
47
 
 *
48
 
 *   Pivot policy:
49
 
 *   (1) Compute thresh = u * max_(i>=j) abs(A_ij);
50
 
 *   (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
51
 
 *           pivot row = k;
52
 
 *       ELSE IF abs(A_jj) >= thresh THEN
53
 
 *           pivot row = j;
54
 
 *       ELSE
55
 
 *           pivot row = m;
56
 
 * 
57
 
 *   Note: If you absolutely want to use a given pivot order, then set u=0.0.
58
 
 *
59
 
 *   Return value: 0      success;
60
 
 *                 i > 0  U(i,i) is exactly zero.
61
 
 *
62
 
 */
63
 
    complex one = {1.0, 0.0};
64
 
    int          fsupc;     /* first column in the supernode */
65
 
    int          nsupc;     /* no of columns in the supernode */
66
 
    int          nsupr;     /* no of rows in the supernode */
67
 
    int          lptr;      /* points to the starting subscript of the supernode */
68
 
    int          pivptr, old_pivptr, diag, diagind;
69
 
    float       pivmax, rtemp, thresh;
70
 
    complex       temp;
71
 
    complex       *lu_sup_ptr; 
72
 
    complex       *lu_col_ptr;
73
 
    int          *lsub_ptr;
74
 
    int          isub, icol, k, itemp;
75
 
    int          *lsub, *xlsub;
76
 
    complex       *lusup;
77
 
    int          *xlusup;
78
 
    extern SuperLUStat_t SuperLUStat;
79
 
    flops_t  *ops = SuperLUStat.ops;
80
 
 
81
 
    /* Initialize pointers */
82
 
    lsub       = Glu->lsub;
83
 
    xlsub      = Glu->xlsub;
84
 
    lusup      = Glu->lusup;
85
 
    xlusup     = Glu->xlusup;
86
 
    fsupc      = (Glu->xsup)[(Glu->supno)[jcol]];
87
 
    nsupc      = jcol - fsupc;          /* excluding jcol; nsupc >= 0 */
88
 
    lptr       = xlsub[fsupc];
89
 
    nsupr      = xlsub[fsupc+1] - lptr;
90
 
    lu_sup_ptr = &lusup[xlusup[fsupc]]; /* start of the current supernode */
91
 
    lu_col_ptr = &lusup[xlusup[jcol]];  /* start of jcol in the supernode */
92
 
    lsub_ptr   = &lsub[lptr];   /* start of row indices of the supernode */
93
 
 
94
 
#ifdef DEBUG
95
 
if ( jcol == MIN_COL ) {
96
 
    printf("Before cdiv: col %d\n", jcol);
97
 
    for (k = nsupc; k < nsupr; k++) 
98
 
        printf("  lu[%d] %f\n", lsub_ptr[k], lu_col_ptr[k]);
99
 
}
100
 
#endif
101
 
    
102
 
    /* Determine the largest abs numerical value for partial pivoting;
103
 
       Also search for user-specified pivot, and diagonal element. */
104
 
    if ( *usepr ) *pivrow = iperm_r[jcol];
105
 
    diagind = iperm_c[jcol];
106
 
    pivmax = 0.0;
107
 
    pivptr = nsupc;
108
 
    diag = EMPTY;
109
 
    old_pivptr = nsupc;
110
 
    for (isub = nsupc; isub < nsupr; ++isub) {
111
 
        rtemp = c_abs1 (&lu_col_ptr[isub]);
112
 
        if ( rtemp > pivmax ) {
113
 
            pivmax = rtemp;
114
 
            pivptr = isub;
115
 
        }
116
 
        if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
117
 
        if ( lsub_ptr[isub] == diagind ) diag = isub;
118
 
    }
119
 
 
120
 
    /* Test for singularity */
121
 
    if ( pivmax == 0.0 ) {
122
 
        *pivrow = lsub_ptr[pivptr];
123
 
        perm_r[*pivrow] = jcol;
124
 
        *usepr = 0;
125
 
        return (jcol+1);
126
 
    }
127
 
 
128
 
    thresh = u * pivmax;
129
 
    
130
 
    /* Choose appropriate pivotal element by our policy. */
131
 
    if ( *usepr ) {
132
 
        rtemp = c_abs1 (&lu_col_ptr[old_pivptr]);
133
 
        if ( rtemp != 0.0 && rtemp >= thresh )
134
 
            pivptr = old_pivptr;
135
 
        else
136
 
            *usepr = 0;
137
 
    }
138
 
    if ( *usepr == 0 ) {
139
 
        /* Use diagonal pivot? */
140
 
        if ( diag >= 0 ) { /* diagonal exists */
141
 
            rtemp = c_abs1 (&lu_col_ptr[diag]);
142
 
            if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = diag;
143
 
        }
144
 
        *pivrow = lsub_ptr[pivptr];
145
 
    }
146
 
    
147
 
    /* Record pivot row */
148
 
    perm_r[*pivrow] = jcol;
149
 
    
150
 
    /* Interchange row subscripts */
151
 
    if ( pivptr != nsupc ) {
152
 
        itemp = lsub_ptr[pivptr];
153
 
        lsub_ptr[pivptr] = lsub_ptr[nsupc];
154
 
        lsub_ptr[nsupc] = itemp;
155
 
 
156
 
        /* Interchange numerical values as well, for the whole snode, such 
157
 
         * that L is indexed the same way as A.
158
 
         */
159
 
        for (icol = 0; icol <= nsupc; icol++) {
160
 
            itemp = pivptr + icol * nsupr;
161
 
            temp = lu_sup_ptr[itemp];
162
 
            lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];
163
 
            lu_sup_ptr[nsupc + icol*nsupr] = temp;
164
 
        }
165
 
    } /* if */
166
 
 
167
 
    /* cdiv operation */
168
 
    ops[FACT] += 10 * (nsupr - nsupc);
169
 
 
170
 
    c_div(&temp, &one, &lu_col_ptr[nsupc]);
171
 
    for (k = nsupc+1; k < nsupr; k++) 
172
 
        cc_mult(&lu_col_ptr[k], &lu_col_ptr[k], &temp);
173
 
 
174
 
    return 0;
175
 
}
176