← 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/dpanel_bmod.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/dpanel_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 dlsolve(int, int, double *, double *);
32
 
void dmatvec(int, int, int, double *, double *, double *);
33
 
extern void dcheck_tempv();
34
 
 
35
 
void
36
 
dpanel_bmod (
37
 
            const int  m,          /* in - number of rows in the matrix */
38
 
            const int  w,          /* in */
39
 
            const int  jcol,       /* in */
40
 
            const int  nseg,       /* in */
41
 
            double     *dense,     /* out, of size n by w */
42
 
            double     *tempv,     /* working array */
43
 
            int        *segrep,    /* in */
44
 
            int        *repfnz,    /* in, of size n by w */
45
 
            GlobalLU_t *Glu        /* modified */
46
 
            )
47
 
{
48
 
/* 
49
 
 * Purpose
50
 
 * =======
51
 
 *
52
 
 *    Performs numeric block updates (sup-panel) in topological order.
53
 
 *    It features: col-col, 2cols-col, 3cols-col, and sup-col updates.
54
 
 *    Special processing on the supernodal portion of L\U[*,j]
55
 
 *
56
 
 *    Before entering this routine, the original nonzeros in the panel 
57
 
 *    were already copied into the spa[m,w].
58
 
 *
59
 
 *    Updated/Output parameters-
60
 
 *      dense[0:m-1,w]: L[*,j:j+w-1] and U[*,j:j+w-1] are returned 
61
 
 *      collectively in the m-by-w vector dense[*]. 
62
 
 *
63
 
 */
64
 
 
65
 
#ifdef USE_VENDOR_BLAS
66
 
#ifdef _CRAY
67
 
    _fcd ftcs1 = _cptofcd("L", strlen("L")),
68
 
         ftcs2 = _cptofcd("N", strlen("N")),
69
 
         ftcs3 = _cptofcd("U", strlen("U"));
70
 
#endif
71
 
    int          incx = 1, incy = 1;
72
 
    double       alpha, beta;
73
 
#endif
74
 
 
75
 
    register int k, ksub;
76
 
    int          fsupc, nsupc, nsupr, nrow;
77
 
    int          krep, krep_ind;
78
 
    double       ukj, ukj1, ukj2;
79
 
    int          luptr, luptr1, luptr2;
80
 
    int          segsze;
81
 
    int          block_nrow;  /* no of rows in a block row */
82
 
    register int lptr;        /* Points to the row subscripts of a supernode */
83
 
    int          kfnz, irow, no_zeros; 
84
 
    register int isub, isub1, i;
85
 
    register int jj;          /* Index through each column in the panel */
86
 
    int          *xsup, *supno;
87
 
    int          *lsub, *xlsub;
88
 
    double       *lusup;
89
 
    int          *xlusup;
90
 
    int          *repfnz_col; /* repfnz[] for a column in the panel */
91
 
    double       *dense_col;  /* dense[] for a column in the panel */
92
 
    double       *tempv1;             /* Used in 1-D update */
93
 
    double       *TriTmp, *MatvecTmp; /* used in 2-D update */
94
 
    double      zero = 0.0;
95
 
    double      one = 1.0;
96
 
    register int ldaTmp;
97
 
    register int r_ind, r_hi;
98
 
    static   int first = 1, maxsuper, rowblk, colblk;
99
 
    extern SuperLUStat_t SuperLUStat;
100
 
    flops_t  *ops = SuperLUStat.ops;
101
 
    
102
 
    xsup    = Glu->xsup;
103
 
    supno   = Glu->supno;
104
 
    lsub    = Glu->lsub;
105
 
    xlsub   = Glu->xlsub;
106
 
    lusup   = Glu->lusup;
107
 
    xlusup  = Glu->xlusup;
108
 
    
109
 
    if ( first ) {
110
 
        maxsuper = sp_ienv(3);
111
 
        rowblk   = sp_ienv(4);
112
 
        colblk   = sp_ienv(5);
113
 
        first = 0;
114
 
    }
115
 
    ldaTmp = maxsuper + rowblk;
116
 
 
117
 
    /* 
118
 
     * For each nonz supernode segment of U[*,j] in topological order 
119
 
     */
120
 
    k = nseg - 1;
121
 
    for (ksub = 0; ksub < nseg; ksub++) { /* for each updating supernode */
122
 
 
123
 
        /* krep = representative of current k-th supernode
124
 
         * fsupc = first supernodal column
125
 
         * nsupc = no of columns in a supernode
126
 
         * nsupr = no of rows in a supernode
127
 
         */
128
 
        krep = segrep[k--];
129
 
        fsupc = xsup[supno[krep]];
130
 
        nsupc = krep - fsupc + 1;
131
 
        nsupr = xlsub[fsupc+1] - xlsub[fsupc];
132
 
        nrow = nsupr - nsupc;
133
 
        lptr = xlsub[fsupc];
134
 
        krep_ind = lptr + nsupc - 1;
135
 
 
136
 
        repfnz_col = repfnz;
137
 
        dense_col = dense;
138
 
        
139
 
        if ( nsupc >= colblk && nrow > rowblk ) { /* 2-D block update */
140
 
 
141
 
            TriTmp = tempv;
142
 
        
143
 
            /* Sequence through each column in panel -- triangular solves */
144
 
            for (jj = jcol; jj < jcol + w; jj++,
145
 
                 repfnz_col += m, dense_col += m, TriTmp += ldaTmp ) {
146
 
 
147
 
                kfnz = repfnz_col[krep];
148
 
                if ( kfnz == EMPTY ) continue;  /* Skip any zero segment */
149
 
            
150
 
                segsze = krep - kfnz + 1;
151
 
                luptr = xlusup[fsupc];
152
 
 
153
 
                ops[TRSV] += segsze * (segsze - 1);
154
 
                ops[GEMV] += 2 * nrow * segsze;
155
 
        
156
 
                /* Case 1: Update U-segment of size 1 -- col-col update */
157
 
                if ( segsze == 1 ) {
158
 
                    ukj = dense_col[lsub[krep_ind]];
159
 
                    luptr += nsupr*(nsupc-1) + nsupc;
160
 
 
161
 
                    for (i = lptr + nsupc; i < xlsub[fsupc+1]; i++) {
162
 
                        irow = lsub[i];
163
 
                        dense_col[irow] -= ukj * lusup[luptr];
164
 
                        ++luptr;
165
 
                    }
166
 
 
167
 
                } else if ( segsze <= 3 ) {
168
 
                    ukj = dense_col[lsub[krep_ind]];
169
 
                    ukj1 = dense_col[lsub[krep_ind - 1]];
170
 
                    luptr += nsupr*(nsupc-1) + nsupc-1;
171
 
                    luptr1 = luptr - nsupr;
172
 
 
173
 
                    if ( segsze == 2 ) {
174
 
                        ukj -= ukj1 * lusup[luptr1];
175
 
                        dense_col[lsub[krep_ind]] = ukj;
176
 
                        for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
177
 
                            irow = lsub[i];
178
 
                            luptr++; luptr1++;
179
 
                            dense_col[irow] -= (ukj*lusup[luptr]
180
 
                                                + ukj1*lusup[luptr1]);
181
 
                        }
182
 
                    } else {
183
 
                        ukj2 = dense_col[lsub[krep_ind - 2]];
184
 
                        luptr2 = luptr1 - nsupr;
185
 
                        ukj1 -= ukj2 * lusup[luptr2-1];
186
 
                        ukj = ukj - ukj1*lusup[luptr1] - ukj2*lusup[luptr2];
187
 
                        dense_col[lsub[krep_ind]] = ukj;
188
 
                        dense_col[lsub[krep_ind-1]] = ukj1;
189
 
                        for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
190
 
                            irow = lsub[i];
191
 
                            luptr++; luptr1++; luptr2++;
192
 
                            dense_col[irow] -= ( ukj*lusup[luptr]
193
 
                             + ukj1*lusup[luptr1] + ukj2*lusup[luptr2] );
194
 
                        }
195
 
                    }
196
 
 
197
 
                } else  {       /* segsze >= 4 */
198
 
                    
199
 
                    /* Copy U[*,j] segment from dense[*] to TriTmp[*], which
200
 
                       holds the result of triangular solves.    */
201
 
                    no_zeros = kfnz - fsupc;
202
 
                    isub = lptr + no_zeros;
203
 
                    for (i = 0; i < segsze; ++i) {
204
 
                        irow = lsub[isub];
205
 
                        TriTmp[i] = dense_col[irow]; /* Gather */
206
 
                        ++isub;
207
 
                    }
208
 
                    
209
 
                    /* 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, TriTmp, &incx );
216
 
#else
217
 
                    dtrsv_( "L", "N", "U", &segsze, &lusup[luptr], 
218
 
                           &nsupr, TriTmp, &incx );
219
 
#endif
220
 
#else           
221
 
                    dlsolve ( nsupr, segsze, &lusup[luptr], TriTmp );
222
 
#endif
223
 
                    
224
 
 
225
 
                } /* else ... */
226
 
            
227
 
            }  /* for jj ... end tri-solves */
228
 
 
229
 
            /* Block row updates; push all the way into dense[*] block */
230
 
            for ( r_ind = 0; r_ind < nrow; r_ind += rowblk ) {
231
 
                
232
 
                r_hi = SUPERLU_MIN(nrow, r_ind + rowblk);
233
 
                block_nrow = SUPERLU_MIN(rowblk, r_hi - r_ind);
234
 
                luptr = xlusup[fsupc] + nsupc + r_ind;
235
 
                isub1 = lptr + nsupc + r_ind;
236
 
                
237
 
                repfnz_col = repfnz;
238
 
                TriTmp = tempv;
239
 
                dense_col = dense;
240
 
                
241
 
                /* Sequence through each column in panel -- matrix-vector */
242
 
                for (jj = jcol; jj < jcol + w; jj++,
243
 
                     repfnz_col += m, dense_col += m, TriTmp += ldaTmp) {
244
 
                    
245
 
                    kfnz = repfnz_col[krep];
246
 
                    if ( kfnz == EMPTY ) continue; /* Skip any zero segment */
247
 
                    
248
 
                    segsze = krep - kfnz + 1;
249
 
                    if ( segsze <= 3 ) continue;   /* skip unrolled cases */
250
 
                    
251
 
                    /* Perform a block update, and scatter the result of
252
 
                       matrix-vector to dense[].                 */
253
 
                    no_zeros = kfnz - fsupc;
254
 
                    luptr1 = luptr + nsupr * no_zeros;
255
 
                    MatvecTmp = &TriTmp[maxsuper];
256
 
                    
257
 
#ifdef USE_VENDOR_BLAS
258
 
                    alpha = one; 
259
 
                    beta = zero;
260
 
#ifdef _CRAY
261
 
                    SGEMV(ftcs2, &block_nrow, &segsze, &alpha, &lusup[luptr1], 
262
 
                           &nsupr, TriTmp, &incx, &beta, MatvecTmp, &incy);
263
 
#else
264
 
                    dgemv_("N", &block_nrow, &segsze, &alpha, &lusup[luptr1], 
265
 
                           &nsupr, TriTmp, &incx, &beta, MatvecTmp, &incy);
266
 
#endif
267
 
#else
268
 
                    dmatvec(nsupr, block_nrow, segsze, &lusup[luptr1],
269
 
                           TriTmp, MatvecTmp);
270
 
#endif
271
 
                    
272
 
                    /* Scatter MatvecTmp[*] into SPA dense[*] temporarily
273
 
                     * such that MatvecTmp[*] can be re-used for the
274
 
                     * the next blok row update. dense[] will be copied into 
275
 
                     * global store after the whole panel has been finished.
276
 
                     */
277
 
                    isub = isub1;
278
 
                    for (i = 0; i < block_nrow; i++) {
279
 
                        irow = lsub[isub];
280
 
                        dense_col[irow] -= MatvecTmp[i];
281
 
                        MatvecTmp[i] = zero;
282
 
                        ++isub;
283
 
                    }
284
 
                    
285
 
                } /* for jj ... */
286
 
                
287
 
            } /* for each block row ... */
288
 
            
289
 
            /* Scatter the triangular solves into SPA dense[*] */
290
 
            repfnz_col = repfnz;
291
 
            TriTmp = tempv;
292
 
            dense_col = dense;
293
 
            
294
 
            for (jj = jcol; jj < jcol + w; jj++,
295
 
                 repfnz_col += m, dense_col += m, TriTmp += ldaTmp) {
296
 
                kfnz = repfnz_col[krep];
297
 
                if ( kfnz == EMPTY ) continue; /* Skip any zero segment */
298
 
                
299
 
                segsze = krep - kfnz + 1;
300
 
                if ( segsze <= 3 ) continue; /* skip unrolled cases */
301
 
                
302
 
                no_zeros = kfnz - fsupc;                
303
 
                isub = lptr + no_zeros;
304
 
                for (i = 0; i < segsze; i++) {
305
 
                    irow = lsub[isub];
306
 
                    dense_col[irow] = TriTmp[i];
307
 
                    TriTmp[i] = zero;
308
 
                    ++isub;
309
 
                }
310
 
                
311
 
            } /* for jj ... */
312
 
            
313
 
        } else { /* 1-D block modification */
314
 
            
315
 
            
316
 
            /* Sequence through each column in the panel */
317
 
            for (jj = jcol; jj < jcol + w; jj++,
318
 
                 repfnz_col += m, dense_col += m) {
319
 
                
320
 
                kfnz = repfnz_col[krep];
321
 
                if ( kfnz == EMPTY ) continue;  /* Skip any zero segment */
322
 
                
323
 
                segsze = krep - kfnz + 1;
324
 
                luptr = xlusup[fsupc];
325
 
 
326
 
                ops[TRSV] += segsze * (segsze - 1);
327
 
                ops[GEMV] += 2 * nrow * segsze;
328
 
                
329
 
                /* Case 1: Update U-segment of size 1 -- col-col update */
330
 
                if ( segsze == 1 ) {
331
 
                    ukj = dense_col[lsub[krep_ind]];
332
 
                    luptr += nsupr*(nsupc-1) + nsupc;
333
 
 
334
 
                    for (i = lptr + nsupc; i < xlsub[fsupc+1]; i++) {
335
 
                        irow = lsub[i];
336
 
                        dense_col[irow] -= ukj * lusup[luptr];
337
 
                        ++luptr;
338
 
                    }
339
 
 
340
 
                } else if ( segsze <= 3 ) {
341
 
                    ukj = dense_col[lsub[krep_ind]];
342
 
                    luptr += nsupr*(nsupc-1) + nsupc-1;
343
 
                    ukj1 = dense_col[lsub[krep_ind - 1]];
344
 
                    luptr1 = luptr - nsupr;
345
 
 
346
 
                    if ( segsze == 2 ) {
347
 
                        ukj -= ukj1 * lusup[luptr1];
348
 
                        dense_col[lsub[krep_ind]] = ukj;
349
 
                        for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
350
 
                            irow = lsub[i];
351
 
                            ++luptr;  ++luptr1;
352
 
                            dense_col[irow] -= (ukj*lusup[luptr]
353
 
                                                + ukj1*lusup[luptr1]);
354
 
                        }
355
 
                    } else {
356
 
                        ukj2 = dense_col[lsub[krep_ind - 2]];
357
 
                        luptr2 = luptr1 - nsupr;
358
 
                        ukj1 -= ukj2 * lusup[luptr2-1];
359
 
                        ukj = ukj - ukj1*lusup[luptr1] - ukj2*lusup[luptr2];
360
 
                        dense_col[lsub[krep_ind]] = ukj;
361
 
                        dense_col[lsub[krep_ind-1]] = ukj1;
362
 
                        for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
363
 
                            irow = lsub[i];
364
 
                            ++luptr; ++luptr1; ++luptr2;
365
 
                            dense_col[irow] -= ( ukj*lusup[luptr]
366
 
                             + ukj1*lusup[luptr1] + ukj2*lusup[luptr2] );
367
 
                        }
368
 
                    }
369
 
 
370
 
                } else  { /* segsze >= 4 */
371
 
                    /* 
372
 
                     * Perform a triangular solve and block update,
373
 
                     * then scatter the result of sup-col update to dense[].
374
 
                     */
375
 
                    no_zeros = kfnz - fsupc;
376
 
                    
377
 
                    /* Copy U[*,j] segment from dense[*] to tempv[*]: 
378
 
                     *    The result of triangular solve is in tempv[*];
379
 
                     *    The result of matrix vector update is in dense_col[*]
380
 
                     */
381
 
                    isub = lptr + no_zeros;
382
 
                    for (i = 0; i < segsze; ++i) {
383
 
                        irow = lsub[isub];
384
 
                        tempv[i] = dense_col[irow]; /* Gather */
385
 
                        ++isub;
386
 
                    }
387
 
                    
388
 
                    /* start effective triangle */
389
 
                    luptr += nsupr * no_zeros + no_zeros;
390
 
                    
391
 
#ifdef USE_VENDOR_BLAS
392
 
#ifdef _CRAY
393
 
                    STRSV( ftcs1, ftcs2, ftcs3, &segsze, &lusup[luptr], 
394
 
                           &nsupr, tempv, &incx );
395
 
#else
396
 
                    dtrsv_( "L", "N", "U", &segsze, &lusup[luptr], 
397
 
                           &nsupr, tempv, &incx );
398
 
#endif
399
 
                    
400
 
                    luptr += segsze;    /* Dense matrix-vector */
401
 
                    tempv1 = &tempv[segsze];
402
 
                    alpha = one;
403
 
                    beta = zero;
404
 
#ifdef _CRAY
405
 
                    SGEMV( ftcs2, &nrow, &segsze, &alpha, &lusup[luptr], 
406
 
                           &nsupr, tempv, &incx, &beta, tempv1, &incy );
407
 
#else
408
 
                    dgemv_( "N", &nrow, &segsze, &alpha, &lusup[luptr], 
409
 
                           &nsupr, tempv, &incx, &beta, tempv1, &incy );
410
 
#endif
411
 
#else
412
 
                    dlsolve ( nsupr, segsze, &lusup[luptr], tempv );
413
 
                    
414
 
                    luptr += segsze;        /* Dense matrix-vector */
415
 
                    tempv1 = &tempv[segsze];
416
 
                    dmatvec (nsupr, nrow, segsze, &lusup[luptr], tempv, tempv1);
417
 
#endif
418
 
                    
419
 
                    /* Scatter tempv[*] into SPA dense[*] temporarily, such
420
 
                     * that tempv[*] can be used for the triangular solve of
421
 
                     * the next column of the panel. They will be copied into 
422
 
                     * ucol[*] after the whole panel has been finished.
423
 
                     */
424
 
                    isub = lptr + no_zeros;
425
 
                    for (i = 0; i < segsze; i++) {
426
 
                        irow = lsub[isub];
427
 
                        dense_col[irow] = tempv[i];
428
 
                        tempv[i] = zero;
429
 
                        isub++;
430
 
                    }
431
 
                    
432
 
                    /* Scatter the update from tempv1[*] into SPA dense[*] */
433
 
                    /* Start dense rectangular L */
434
 
                    for (i = 0; i < nrow; i++) {
435
 
                        irow = lsub[isub];
436
 
                        dense_col[irow] -= tempv1[i];
437
 
                        tempv1[i] = zero;
438
 
                        ++isub; 
439
 
                    }
440
 
                    
441
 
                } /* else segsze>=4 ... */
442
 
                
443
 
            } /* for each column in the panel... */
444
 
            
445
 
        } /* else 1-D update ... */
446
 
 
447
 
    } /* for each updating supernode ... */
448
 
 
449
 
}
450
 
 
451
 
 
452