← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/superlu/dcolumn_bmod.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/sandbox/pysparse/superlu/dcolumn_bmod.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 <stdio.h>
 
24
#include <stdlib.h>
 
25
#include "dsp_defs.h"
 
26
#include "util.h"
 
27
 
 
28
/* 
 
29
 * Function prototypes 
 
30
 */
 
31
void dusolve(int, int, double*, double*);
 
32
void dlsolve(int, int, double*, double*);
 
33
void dmatvec(int, int, int, double*, double*, double*);
 
34
 
 
35
 
 
36
 
 
37
/* Return value:   0 - successful return
 
38
 *               > 0 - number of bytes allocated when run out of space
 
39
 */
 
40
int
 
41
dcolumn_bmod (
 
42
             const int  jcol,     /* in */
 
43
             const int  nseg,     /* in */
 
44
             double     *dense,   /* in */
 
45
             double     *tempv,   /* working array */
 
46
             int        *segrep,  /* in */
 
47
             int        *repfnz,  /* in */
 
48
             int        fpanelc,  /* in -- first column in the current panel */
 
49
             GlobalLU_t *Glu      /* modified */
 
50
             )
 
51
{
 
52
/*
 
53
 * Purpose:
 
54
 * ========
 
55
 *    Performs numeric block updates (sup-col) in topological order.
 
56
 *    It features: col-col, 2cols-col, 3cols-col, and sup-col updates.
 
57
 *    Special processing on the supernodal portion of L\U[*,j]
 
58
 *
 
59
 */
 
60
#ifdef _CRAY
 
61
    _fcd ftcs1 = _cptofcd("L", strlen("L")),
 
62
         ftcs2 = _cptofcd("N", strlen("N")),
 
63
         ftcs3 = _cptofcd("U", strlen("U"));
 
64
#endif
 
65
    int         incx = 1, incy = 1;
 
66
    double      alpha, beta;
 
67
    
 
68
    /* krep = representative of current k-th supernode
 
69
     * fsupc = first supernodal column
 
70
     * nsupc = no of columns in supernode
 
71
     * nsupr = no of rows in supernode (used as leading dimension)
 
72
     * luptr = location of supernodal LU-block in storage
 
73
     * kfnz = first nonz in the k-th supernodal segment
 
74
     * no_zeros = no of leading zeros in a supernodal U-segment
 
75
     */
 
76
    double       ukj, ukj1, ukj2;
 
77
    int          luptr, luptr1, luptr2;
 
78
    int          fsupc, nsupc, nsupr, segsze;
 
79
    int          nrow;    /* No of rows in the matrix of matrix-vector */
 
80
    int          jcolp1, jsupno, k, ksub, krep, krep_ind, ksupno;
 
81
    register int lptr, kfnz, isub, irow, i;
 
82
    register int no_zeros, new_next; 
 
83
    int          ufirst, nextlu;
 
84
    int          fst_col; /* First column within small LU update */
 
85
    int          d_fsupc; /* Distance between the first column of the current
 
86
                             panel and the first column of the current snode. */
 
87
    int          *xsup, *supno;
 
88
    int          *lsub, *xlsub;
 
89
    double       *lusup;
 
90
    int          *xlusup;
 
91
    int          nzlumax;
 
92
    double       *tempv1;
 
93
    double      zero = 0.0;
 
94
    double      one = 1.0;
 
95
    double      none = -1.0;
 
96
    int          mem_error;
 
97
    extern SuperLUStat_t SuperLUStat;
 
98
    flops_t  *ops = SuperLUStat.ops;
 
99
 
 
100
    xsup    = Glu->xsup;
 
101
    supno   = Glu->supno;
 
102
    lsub    = Glu->lsub;
 
103
    xlsub   = Glu->xlsub;
 
104
    lusup   = Glu->lusup;
 
105
    xlusup  = Glu->xlusup;
 
106
    nzlumax = Glu->nzlumax;
 
107
    jcolp1 = jcol + 1;
 
108
    jsupno = supno[jcol];
 
109
    
 
110
    /* 
 
111
     * For each nonz supernode segment of U[*,j] in topological order 
 
112
     */
 
113
    k = nseg - 1;
 
114
    for (ksub = 0; ksub < nseg; ksub++) {
 
115
 
 
116
        krep = segrep[k];
 
117
        k--;
 
118
        ksupno = supno[krep];
 
119
        if ( jsupno != ksupno ) { /* Outside the rectangular supernode */
 
120
 
 
121
            fsupc = xsup[ksupno];
 
122
            fst_col = SUPERLU_MAX ( fsupc, fpanelc );
 
123
 
 
124
            /* Distance from the current supernode to the current panel; 
 
125
               d_fsupc=0 if fsupc > fpanelc. */
 
126
            d_fsupc = fst_col - fsupc; 
 
127
 
 
128
            luptr = xlusup[fst_col] + d_fsupc;
 
129
            lptr = xlsub[fsupc] + d_fsupc;
 
130
 
 
131
            kfnz = repfnz[krep];
 
132
            kfnz = SUPERLU_MAX ( kfnz, fpanelc );
 
133
 
 
134
            segsze = krep - kfnz + 1;
 
135
            nsupc = krep - fst_col + 1;
 
136
            nsupr = xlsub[fsupc+1] - xlsub[fsupc];      /* Leading dimension */
 
137
            nrow = nsupr - d_fsupc - nsupc;
 
138
            krep_ind = lptr + nsupc - 1;
 
139
 
 
140
            ops[TRSV] += segsze * (segsze - 1);
 
141
            ops[GEMV] += 2 * nrow * segsze;
 
142
 
 
143
 
 
144
            /* 
 
145
             * Case 1: Update U-segment of size 1 -- col-col update 
 
146
             */
 
147
            if ( segsze == 1 ) {
 
148
                ukj = dense[lsub[krep_ind]];
 
149
                luptr += nsupr*(nsupc-1) + nsupc;
 
150
 
 
151
                for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
 
152
                    irow = lsub[i];
 
153
                    dense[irow] -=  ukj*lusup[luptr];
 
154
                    luptr++;
 
155
                }
 
156
 
 
157
            } else if ( segsze <= 3 ) {
 
158
                ukj = dense[lsub[krep_ind]];
 
159
                luptr += nsupr*(nsupc-1) + nsupc-1;
 
160
                ukj1 = dense[lsub[krep_ind - 1]];
 
161
                luptr1 = luptr - nsupr;
 
162
 
 
163
                if ( segsze == 2 ) { /* Case 2: 2cols-col update */
 
164
                    ukj -= ukj1 * lusup[luptr1];
 
165
                    dense[lsub[krep_ind]] = ukj;
 
166
                    for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
 
167
                        irow = lsub[i];
 
168
                        luptr++;
 
169
                        luptr1++;
 
170
                        dense[irow] -= ( ukj*lusup[luptr]
 
171
                                        + ukj1*lusup[luptr1] );
 
172
                    }
 
173
                } else { /* Case 3: 3cols-col update */
 
174
                    ukj2 = dense[lsub[krep_ind - 2]];
 
175
                    luptr2 = luptr1 - nsupr;
 
176
                    ukj1 -= ukj2 * lusup[luptr2-1];
 
177
                    ukj = ukj - ukj1*lusup[luptr1] - ukj2*lusup[luptr2];
 
178
                    dense[lsub[krep_ind]] = ukj;
 
179
                    dense[lsub[krep_ind-1]] = ukj1;
 
180
                    for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
 
181
                        irow = lsub[i];
 
182
                        luptr++;
 
183
                        luptr1++;
 
184
                        luptr2++;
 
185
                        dense[irow] -= ( ukj*lusup[luptr]
 
186
                             + ukj1*lusup[luptr1] + ukj2*lusup[luptr2] );
 
187
                    }
 
188
                }
 
189
 
 
190
 
 
191
 
 
192
            } else {
 
193
                /*
 
194
                 * Case: sup-col update
 
195
                 * Perform a triangular solve and block update,
 
196
                 * then scatter the result of sup-col update to dense
 
197
                 */
 
198
 
 
199
                no_zeros = kfnz - fst_col;
 
200
 
 
201
                /* Copy U[*,j] segment from dense[*] to tempv[*] */
 
202
                isub = lptr + no_zeros;
 
203
                for (i = 0; i < segsze; i++) {
 
204
                    irow = lsub[isub];
 
205
                    tempv[i] = dense[irow];
 
206
                    ++isub; 
 
207
                }
 
208
 
 
209
                /* Dense triangular solve -- start effective triangle */
 
210
                luptr += nsupr * no_zeros + no_zeros; 
 
211
                
 
212
#ifdef USE_VENDOR_BLAS
 
213
#ifdef _CRAY
 
214
                STRSV( ftcs1, ftcs2, ftcs3, &segsze, &lusup[luptr], 
 
215
                       &nsupr, tempv, &incx );
 
216
#else           
 
217
                dtrsv_( "L", "N", "U", &segsze, &lusup[luptr], 
 
218
                       &nsupr, tempv, &incx );
 
219
#endif          
 
220
                luptr += segsze;  /* Dense matrix-vector */
 
221
                tempv1 = &tempv[segsze];
 
222
                alpha = one;
 
223
                beta = zero;
 
224
#ifdef _CRAY
 
225
                SGEMV( ftcs2, &nrow, &segsze, &alpha, &lusup[luptr], 
 
226
                       &nsupr, tempv, &incx, &beta, tempv1, &incy );
 
227
#else
 
228
                dgemv_( "N", &nrow, &segsze, &alpha, &lusup[luptr], 
 
229
                       &nsupr, tempv, &incx, &beta, tempv1, &incy );
 
230
#endif
 
231
#else
 
232
                dlsolve ( nsupr, segsze, &lusup[luptr], tempv );
 
233
 
 
234
                luptr += segsze;  /* Dense matrix-vector */
 
235
                tempv1 = &tempv[segsze];
 
236
                dmatvec (nsupr, nrow , segsze, &lusup[luptr], tempv, tempv1);
 
237
#endif
 
238
                
 
239
                
 
240
                /* Scatter tempv[] into SPA dense[] as a temporary storage */
 
241
                isub = lptr + no_zeros;
 
242
                for (i = 0; i < segsze; i++) {
 
243
                    irow = lsub[isub];
 
244
                    dense[irow] = tempv[i];
 
245
                    tempv[i] = zero;
 
246
                    ++isub;
 
247
                }
 
248
 
 
249
                /* Scatter tempv1[] into SPA dense[] */
 
250
                for (i = 0; i < nrow; i++) {
 
251
                    irow = lsub[isub];
 
252
                    dense[irow] -= tempv1[i];
 
253
                    tempv1[i] = zero;
 
254
                    ++isub;
 
255
                }
 
256
            }
 
257
            
 
258
        } /* if jsupno ... */
 
259
 
 
260
    } /* for each segment... */
 
261
 
 
262
    /*
 
263
     *  Process the supernodal portion of L\U[*,j]
 
264
     */
 
265
    nextlu = xlusup[jcol];
 
266
    fsupc = xsup[jsupno];
 
267
 
 
268
    /* Copy the SPA dense into L\U[*,j] */
 
269
    new_next = nextlu + xlsub[fsupc+1] - xlsub[fsupc];
 
270
    while ( new_next > nzlumax ) {
 
271
        if (mem_error = dLUMemXpand(jcol, nextlu, LUSUP, &nzlumax, Glu))
 
272
            return (mem_error);
 
273
        lusup = Glu->lusup;
 
274
        lsub = Glu->lsub;
 
275
    }
 
276
 
 
277
    for (isub = xlsub[fsupc]; isub < xlsub[fsupc+1]; isub++) {
 
278
        irow = lsub[isub];
 
279
        lusup[nextlu] = dense[irow];
 
280
        dense[irow] = zero;
 
281
        ++nextlu;
 
282
    }
 
283
 
 
284
    xlusup[jcolp1] = nextlu;    /* Close L\U[*,jcol] */
 
285
 
 
286
    /* For more updates within the panel (also within the current supernode), 
 
287
     * should start from the first column of the panel, or the first column 
 
288
     * of the supernode, whichever is bigger. There are 2 cases:
 
289
     *    1) fsupc < fpanelc, then fst_col := fpanelc
 
290
     *    2) fsupc >= fpanelc, then fst_col := fsupc
 
291
     */
 
292
    fst_col = SUPERLU_MAX ( fsupc, fpanelc );
 
293
 
 
294
    if ( fst_col < jcol ) {
 
295
 
 
296
        /* Distance between the current supernode and the current panel.
 
297
           d_fsupc=0 if fsupc >= fpanelc. */
 
298
        d_fsupc = fst_col - fsupc;
 
299
 
 
300
        lptr = xlsub[fsupc] + d_fsupc;
 
301
        luptr = xlusup[fst_col] + d_fsupc;
 
302
        nsupr = xlsub[fsupc+1] - xlsub[fsupc];  /* Leading dimension */
 
303
        nsupc = jcol - fst_col; /* Excluding jcol */
 
304
        nrow = nsupr - d_fsupc - nsupc;
 
305
 
 
306
        /* Points to the beginning of jcol in snode L\U(jsupno) */
 
307
        ufirst = xlusup[jcol] + d_fsupc;        
 
308
 
 
309
        ops[TRSV] += nsupc * (nsupc - 1);
 
310
        ops[GEMV] += 2 * nrow * nsupc;
 
311
        
 
312
#ifdef USE_VENDOR_BLAS
 
313
#ifdef _CRAY
 
314
        STRSV( ftcs1, ftcs2, ftcs3, &nsupc, &lusup[luptr], 
 
315
               &nsupr, &lusup[ufirst], &incx );
 
316
#else
 
317
        dtrsv_( "L", "N", "U", &nsupc, &lusup[luptr], 
 
318
               &nsupr, &lusup[ufirst], &incx );
 
319
#endif
 
320
        
 
321
        alpha = none; beta = one; /* y := beta*y + alpha*A*x */
 
322
 
 
323
#ifdef _CRAY
 
324
        SGEMV( ftcs2, &nrow, &nsupc, &alpha, &lusup[luptr+nsupc], &nsupr,
 
325
               &lusup[ufirst], &incx, &beta, &lusup[ufirst+nsupc], &incy );
 
326
#else
 
327
        dgemv_( "N", &nrow, &nsupc, &alpha, &lusup[luptr+nsupc], &nsupr,
 
328
               &lusup[ufirst], &incx, &beta, &lusup[ufirst+nsupc], &incy );
 
329
#endif
 
330
#else
 
331
        dlsolve ( nsupr, nsupc, &lusup[luptr], &lusup[ufirst] );
 
332
 
 
333
        dmatvec ( nsupr, nrow, nsupc, &lusup[luptr+nsupc],
 
334
                &lusup[ufirst], tempv );
 
335
        
 
336
        /* Copy updates from tempv[*] into lusup[*] */
 
337
        isub = ufirst + nsupc;
 
338
        for (i = 0; i < nrow; i++) {
 
339
            lusup[isub] -= tempv[i];
 
340
            tempv[i] = 0.0;
 
341
            ++isub;
 
342
        }
 
343
 
 
344
#endif
 
345
        
 
346
        
 
347
    } /* if fst_col < jcol ... */ 
 
348
 
 
349
    return 0;
 
350
}