← Back to branch summary

~ubuntu-branches/ubuntu/lucid/python-scipy/lucid

« back to all changes in this revision

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