← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

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