← 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/ssp_blas2.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/ssp_blas2.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
 * File name:           sp_blas2.c
 
12
 * Purpose:             Sparse BLAS 2, using some dense BLAS 2 operations.
 
13
 */
 
14
 
 
15
#include "ssp_defs.h"
 
16
#include "util.h"
 
17
 
 
18
 
 
19
/* 
 
20
 * Function prototypes 
 
21
 */
 
22
void susolve(int, int, float*, float*);
 
23
void slsolve(int, int, float*, float*);
 
24
void smatvec(int, int, int, float*, float*, float*);
 
25
 
 
26
 
 
27
int
 
28
sp_strsv(char *uplo, char *trans, char *diag, SuperMatrix *L, 
 
29
         SuperMatrix *U, float *x, int *info)
 
30
{
 
31
/*
 
32
 *   Purpose
 
33
 *   =======
 
34
 *
 
35
 *   sp_strsv() solves one of the systems of equations   
 
36
 *       A*x = b,   or   A'*x = b,
 
37
 *   where b and x are n element vectors and A is a sparse unit , or   
 
38
 *   non-unit, upper or lower triangular matrix.   
 
39
 *   No test for singularity or near-singularity is included in this   
 
40
 *   routine. Such tests must be performed before calling this routine.   
 
41
 *
 
42
 *   Parameters   
 
43
 *   ==========   
 
44
 *
 
45
 *   uplo   - (input) char*
 
46
 *            On entry, uplo specifies whether the matrix is an upper or   
 
47
 *             lower triangular matrix as follows:   
 
48
 *                uplo = 'U' or 'u'   A is an upper triangular matrix.   
 
49
 *                uplo = 'L' or 'l'   A is a lower triangular matrix.   
 
50
 *
 
51
 *   trans  - (input) char*
 
52
 *             On entry, trans specifies the equations to be solved as   
 
53
 *             follows:   
 
54
 *                trans = 'N' or 'n'   A*x = b.   
 
55
 *                trans = 'T' or 't'   A'*x = b.   
 
56
 *                trans = 'C' or 'c'   A'*x = b.   
 
57
 *
 
58
 *   diag   - (input) char*
 
59
 *             On entry, diag specifies whether or not A is unit   
 
60
 *             triangular as follows:   
 
61
 *                diag = 'U' or 'u'   A is assumed to be unit triangular.   
 
62
 *                diag = 'N' or 'n'   A is not assumed to be unit   
 
63
 *                                    triangular.   
 
64
 *           
 
65
 *   L       - (input) SuperMatrix*
 
66
 *             The factor L from the factorization Pr*A*Pc=L*U. Use
 
67
 *             compressed row subscripts storage for supernodes,
 
68
 *             i.e., L has types: Stype = SC, Dtype = S_, Mtype = TRLU.
 
69
 *
 
70
 *   U       - (input) SuperMatrix*
 
71
 *              The factor U from the factorization Pr*A*Pc=L*U.
 
72
 *              U has types: Stype = NC, Dtype = S_, Mtype = TRU.
 
73
 *    
 
74
 *   x       - (input/output) float*
 
75
 *             Before entry, the incremented array X must contain the n   
 
76
 *             element right-hand side vector b. On exit, X is overwritten 
 
77
 *             with the solution vector x.
 
78
 *
 
79
 *   info    - (output) int*
 
80
 *             If *info = -i, the i-th argument had an illegal value.
 
81
 *
 
82
 */
 
83
#ifdef _CRAY
 
84
    _fcd ftcs1 = _cptofcd("L", strlen("L")),
 
85
         ftcs2 = _cptofcd("N", strlen("N")),
 
86
         ftcs3 = _cptofcd("U", strlen("U"));
 
87
#endif
 
88
    SCformat *Lstore;
 
89
    NCformat *Ustore;
 
90
    float   *Lval, *Uval;
 
91
    int incx = 1, incy = 1;
 
92
    float alpha = 1.0, beta = 1.0;
 
93
    int nrow;
 
94
    int fsupc, nsupr, nsupc, luptr, istart, irow;
 
95
    int i, k, iptr, jcol;
 
96
    float *work;
 
97
    flops_t solve_ops;
 
98
    extern SuperLUStat_t SuperLUStat;
 
99
 
 
100
    /* Test the input parameters */
 
101
    *info = 0;
 
102
    if ( !lsame_(uplo,"L") && !lsame_(uplo, "U") ) *info = -1;
 
103
    else if ( !lsame_(trans, "N") && !lsame_(trans, "T") ) *info = -2;
 
104
    else if ( !lsame_(diag, "U") && !lsame_(diag, "N") ) *info = -3;
 
105
    else if ( L->nrow != L->ncol || L->nrow < 0 ) *info = -4;
 
106
    else if ( U->nrow != U->ncol || U->nrow < 0 ) *info = -5;
 
107
    if ( *info ) {
 
108
        i = -(*info);
 
109
        xerbla_("sp_strsv", &i);
 
110
        return 0;
 
111
    }
 
112
 
 
113
    Lstore = L->Store;
 
114
    Lval = Lstore->nzval;
 
115
    Ustore = U->Store;
 
116
    Uval = Ustore->nzval;
 
117
    solve_ops = 0;
 
118
 
 
119
    if ( !(work = floatCalloc(L->nrow)) )
 
120
        ABORT("Malloc fails for work in sp_strsv().");
 
121
    
 
122
    if ( lsame_(trans, "N") ) { /* Form x := inv(A)*x. */
 
123
        
 
124
        if ( lsame_(uplo, "L") ) {
 
125
            /* Form x := inv(L)*x */
 
126
            if ( L->nrow == 0 ) return 0; /* Quick return */
 
127
            
 
128
            for (k = 0; k <= Lstore->nsuper; k++) {
 
129
                fsupc = L_FST_SUPC(k);
 
130
                istart = L_SUB_START(fsupc);
 
131
                nsupr = L_SUB_START(fsupc+1) - istart;
 
132
                nsupc = L_FST_SUPC(k+1) - fsupc;
 
133
                luptr = L_NZ_START(fsupc);
 
134
                nrow = nsupr - nsupc;
 
135
 
 
136
                solve_ops += nsupc * (nsupc - 1);
 
137
                solve_ops += 2 * nrow * nsupc;
 
138
 
 
139
                if ( nsupc == 1 ) {
 
140
                    for (iptr=istart+1; iptr < L_SUB_START(fsupc+1); ++iptr) {
 
141
                        irow = L_SUB(iptr);
 
142
                        ++luptr;
 
143
                        x[irow] -= x[fsupc] * Lval[luptr];
 
144
                    }
 
145
                } else {
 
146
#ifdef USE_VENDOR_BLAS
 
147
#ifdef _CRAY
 
148
                    STRSV(ftcs1, ftcs2, ftcs3, &nsupc, &Lval[luptr], &nsupr,
 
149
                        &x[fsupc], &incx);
 
150
                
 
151
                    SGEMV(ftcs2, &nrow, &nsupc, &alpha, &Lval[luptr+nsupc], 
 
152
                        &nsupr, &x[fsupc], &incx, &beta, &work[0], &incy);
 
153
#else
 
154
                    strsv_("L", "N", "U", &nsupc, &Lval[luptr], &nsupr,
 
155
                        &x[fsupc], &incx);
 
156
                
 
157
                    sgemv_("N", &nrow, &nsupc, &alpha, &Lval[luptr+nsupc], 
 
158
                        &nsupr, &x[fsupc], &incx, &beta, &work[0], &incy);
 
159
#endif
 
160
#else
 
161
                    slsolve ( nsupr, nsupc, &Lval[luptr], &x[fsupc]);
 
162
                
 
163
                    smatvec ( nsupr, nsupr-nsupc, nsupc, &Lval[luptr+nsupc],
 
164
                        &x[fsupc], &work[0] );
 
165
#endif          
 
166
                
 
167
                    iptr = istart + nsupc;
 
168
                    for (i = 0; i < nrow; ++i, ++iptr) {
 
169
                        irow = L_SUB(iptr);
 
170
                        x[irow] -= work[i];     /* Scatter */
 
171
                        work[i] = 0.0;
 
172
 
 
173
                    }
 
174
                }
 
175
            } /* for k ... */
 
176
            
 
177
        } else {
 
178
            /* Form x := inv(U)*x */
 
179
            
 
180
            if ( U->nrow == 0 ) return 0; /* Quick return */
 
181
            
 
182
            for (k = Lstore->nsuper; k >= 0; k--) {
 
183
                fsupc = L_FST_SUPC(k);
 
184
                nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);
 
185
                nsupc = L_FST_SUPC(k+1) - fsupc;
 
186
                luptr = L_NZ_START(fsupc);
 
187
                
 
188
                solve_ops += nsupc * (nsupc + 1);
 
189
 
 
190
                if ( nsupc == 1 ) {
 
191
                    x[fsupc] /= Lval[luptr];
 
192
                    for (i = U_NZ_START(fsupc); i < U_NZ_START(fsupc+1); ++i) {
 
193
                        irow = U_SUB(i);
 
194
                        x[irow] -= x[fsupc] * Uval[i];
 
195
                    }
 
196
                } else {
 
197
#ifdef USE_VENDOR_BLAS
 
198
#ifdef _CRAY
 
199
                    STRSV(ftcs3, ftcs2, ftcs2, &nsupc, &Lval[luptr], &nsupr,
 
200
                       &x[fsupc], &incx);
 
201
#else
 
202
                    strsv_("U", "N", "N", &nsupc, &Lval[luptr], &nsupr,
 
203
                       &x[fsupc], &incx);
 
204
#endif
 
205
#else           
 
206
                    susolve ( nsupr, nsupc, &Lval[luptr], &x[fsupc] );
 
207
#endif          
 
208
 
 
209
                    for (jcol = fsupc; jcol < L_FST_SUPC(k+1); jcol++) {
 
210
                        solve_ops += 2*(U_NZ_START(jcol+1) - U_NZ_START(jcol));
 
211
                        for (i = U_NZ_START(jcol); i < U_NZ_START(jcol+1); 
 
212
                                i++) {
 
213
                            irow = U_SUB(i);
 
214
                            x[irow] -= x[jcol] * Uval[i];
 
215
                        }
 
216
                    }
 
217
                }
 
218
            } /* for k ... */
 
219
            
 
220
        }
 
221
    } else { /* Form x := inv(A')*x */
 
222
        
 
223
        if ( lsame_(uplo, "L") ) {
 
224
            /* Form x := inv(L')*x */
 
225
            if ( L->nrow == 0 ) return 0; /* Quick return */
 
226
            
 
227
            for (k = Lstore->nsuper; k >= 0; --k) {
 
228
                fsupc = L_FST_SUPC(k);
 
229
                istart = L_SUB_START(fsupc);
 
230
                nsupr = L_SUB_START(fsupc+1) - istart;
 
231
                nsupc = L_FST_SUPC(k+1) - fsupc;
 
232
                luptr = L_NZ_START(fsupc);
 
233
 
 
234
                solve_ops += 2 * (nsupr - nsupc) * nsupc;
 
235
 
 
236
                for (jcol = fsupc; jcol < L_FST_SUPC(k+1); jcol++) {
 
237
                    iptr = istart + nsupc;
 
238
                    for (i = L_NZ_START(jcol) + nsupc; 
 
239
                                i < L_NZ_START(jcol+1); i++) {
 
240
                        irow = L_SUB(iptr);
 
241
                        x[jcol] -= x[irow] * Lval[i];
 
242
                        iptr++;
 
243
                    }
 
244
                }
 
245
                
 
246
                if ( nsupc > 1 ) {
 
247
                    solve_ops += nsupc * (nsupc - 1);
 
248
#ifdef _CRAY
 
249
                    ftcs1 = _cptofcd("L", strlen("L"));
 
250
                    ftcs2 = _cptofcd("T", strlen("T"));
 
251
                    ftcs3 = _cptofcd("U", strlen("U"));
 
252
                    STRSV(ftcs1, ftcs2, ftcs3, &nsupc, &Lval[luptr], &nsupr,
 
253
                        &x[fsupc], &incx);
 
254
#else
 
255
                    strsv_("L", "T", "U", &nsupc, &Lval[luptr], &nsupr,
 
256
                        &x[fsupc], &incx);
 
257
#endif
 
258
                }
 
259
            }
 
260
        } else {
 
261
            /* Form x := inv(U')*x */
 
262
            if ( U->nrow == 0 ) return 0; /* Quick return */
 
263
            
 
264
            for (k = 0; k <= Lstore->nsuper; k++) {
 
265
                fsupc = L_FST_SUPC(k);
 
266
                nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);
 
267
                nsupc = L_FST_SUPC(k+1) - fsupc;
 
268
                luptr = L_NZ_START(fsupc);
 
269
 
 
270
                for (jcol = fsupc; jcol < L_FST_SUPC(k+1); jcol++) {
 
271
                    solve_ops += 2*(U_NZ_START(jcol+1) - U_NZ_START(jcol));
 
272
                    for (i = U_NZ_START(jcol); i < U_NZ_START(jcol+1); i++) {
 
273
                        irow = U_SUB(i);
 
274
                        x[jcol] -= x[irow] * Uval[i];
 
275
                    }
 
276
                }
 
277
 
 
278
                solve_ops += nsupc * (nsupc + 1);
 
279
 
 
280
                if ( nsupc == 1 ) {
 
281
                    x[fsupc] /= Lval[luptr];
 
282
                } else {
 
283
#ifdef _CRAY
 
284
                    ftcs1 = _cptofcd("U", strlen("U"));
 
285
                    ftcs2 = _cptofcd("T", strlen("T"));
 
286
                    ftcs3 = _cptofcd("N", strlen("N"));
 
287
                    STRSV( ftcs1, ftcs2, ftcs3, &nsupc, &Lval[luptr], &nsupr,
 
288
                            &x[fsupc], &incx);
 
289
#else
 
290
                    strsv_("U", "T", "N", &nsupc, &Lval[luptr], &nsupr,
 
291
                            &x[fsupc], &incx);
 
292
#endif
 
293
                }
 
294
            } /* for k ... */
 
295
        }
 
296
    }
 
297
 
 
298
    SuperLUStat.ops[SOLVE] += solve_ops;
 
299
    SUPERLU_FREE(work);
 
300
    return 0;
 
301
}
 
302
 
 
303
 
 
304
 
 
305
 
 
306
int
 
307
sp_sgemv(char *trans, float alpha, SuperMatrix *A, float *x, 
 
308
         int incx, float beta, float *y, int incy)
 
309
{
 
310
/*  Purpose   
 
311
    =======   
 
312
 
 
313
    sp_sgemv()  performs one of the matrix-vector operations   
 
314
       y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,   
 
315
    where alpha and beta are scalars, x and y are vectors and A is a
 
316
    sparse A->nrow by A->ncol matrix.   
 
317
 
 
318
    Parameters   
 
319
    ==========   
 
320
 
 
321
    TRANS  - (input) char*
 
322
             On entry, TRANS specifies the operation to be performed as   
 
323
             follows:   
 
324
                TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.   
 
325
                TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.   
 
326
                TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.   
 
327
 
 
328
    ALPHA  - (input) float
 
329
             On entry, ALPHA specifies the scalar alpha.   
 
330
 
 
331
    A      - (input) SuperMatrix*
 
332
             Matrix A with a sparse format, of dimension (A->nrow, A->ncol).
 
333
             Currently, the type of A can be:
 
334
                 Stype = NC or NCP; Dtype = S_; Mtype = GE. 
 
335
             In the future, more general A can be handled.
 
336
 
 
337
    X      - (input) float*, array of DIMENSION at least   
 
338
             ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'   
 
339
             and at least   
 
340
             ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.   
 
341
             Before entry, the incremented array X must contain the   
 
342
             vector x.   
 
343
 
 
344
    INCX   - (input) int
 
345
             On entry, INCX specifies the increment for the elements of   
 
346
             X. INCX must not be zero.   
 
347
 
 
348
    BETA   - (input) float
 
349
             On entry, BETA specifies the scalar beta. When BETA is   
 
350
             supplied as zero then Y need not be set on input.   
 
351
 
 
352
    Y      - (output) float*,  array of DIMENSION at least   
 
353
             ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'   
 
354
             and at least   
 
355
             ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.   
 
356
             Before entry with BETA non-zero, the incremented array Y   
 
357
             must contain the vector y. On exit, Y is overwritten by the 
 
358
             updated vector y.
 
359
             
 
360
    INCY   - (input) int
 
361
             On entry, INCY specifies the increment for the elements of   
 
362
             Y. INCY must not be zero.   
 
363
 
 
364
    ==== Sparse Level 2 Blas routine.   
 
365
*/
 
366
 
 
367
    /* Local variables */
 
368
    NCformat *Astore;
 
369
    float   *Aval;
 
370
    int info;
 
371
    float temp;
 
372
    int lenx, leny, i, j, irow;
 
373
    int iy, jx, jy, kx, ky;
 
374
    int notran;
 
375
 
 
376
    notran = lsame_(trans, "N");
 
377
    Astore = A->Store;
 
378
    Aval = Astore->nzval;
 
379
    
 
380
    /* Test the input parameters */
 
381
    info = 0;
 
382
    if ( !notran && !lsame_(trans, "T") && !lsame_(trans, "C")) info = 1;
 
383
    else if ( A->nrow < 0 || A->ncol < 0 ) info = 3;
 
384
    else if (incx == 0) info = 5;
 
385
    else if (incy == 0) info = 8;
 
386
    if (info != 0) {
 
387
        xerbla_("sp_sgemv ", &info);
 
388
        return 0;
 
389
    }
 
390
 
 
391
    /* Quick return if possible. */
 
392
    if (A->nrow == 0 || A->ncol == 0 || alpha == 0. && beta == 1.)
 
393
        return 0;
 
394
 
 
395
    /* Set  LENX  and  LENY, the lengths of the vectors x and y, and set 
 
396
       up the start points in  X  and  Y. */
 
397
    if (lsame_(trans, "N")) {
 
398
        lenx = A->ncol;
 
399
        leny = A->nrow;
 
400
    } else {
 
401
        lenx = A->nrow;
 
402
        leny = A->ncol;
 
403
    }
 
404
    if (incx > 0) kx = 0;
 
405
    else kx =  - (lenx - 1) * incx;
 
406
    if (incy > 0) ky = 0;
 
407
    else ky =  - (leny - 1) * incy;
 
408
 
 
409
    /* Start the operations. In this version the elements of A are   
 
410
       accessed sequentially with one pass through A. */
 
411
    /* First form  y := beta*y. */
 
412
    if (beta != 1.) {
 
413
        if (incy == 1) {
 
414
            if (beta == 0.)
 
415
                for (i = 0; i < leny; ++i) y[i] = 0.;
 
416
            else
 
417
                for (i = 0; i < leny; ++i) y[i] = beta * y[i];
 
418
        } else {
 
419
            iy = ky;
 
420
            if (beta == 0.)
 
421
                for (i = 0; i < leny; ++i) {
 
422
                    y[iy] = 0.;
 
423
                    iy += incy;
 
424
                }
 
425
            else
 
426
                for (i = 0; i < leny; ++i) {
 
427
                    y[iy] = beta * y[iy];
 
428
                    iy += incy;
 
429
                }
 
430
        }
 
431
    }
 
432
    
 
433
    if (alpha == 0.) return 0;
 
434
 
 
435
    if ( notran ) {
 
436
        /* Form  y := alpha*A*x + y. */
 
437
        jx = kx;
 
438
        if (incy == 1) {
 
439
            for (j = 0; j < A->ncol; ++j) {
 
440
                if (x[jx] != 0.) {
 
441
                    temp = alpha * x[jx];
 
442
                    for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) {
 
443
                        irow = Astore->rowind[i];
 
444
                        y[irow] += temp * Aval[i];
 
445
                    }
 
446
                }
 
447
                jx += incx;
 
448
            }
 
449
        } else {
 
450
            ABORT("Not implemented.");
 
451
        }
 
452
    } else {
 
453
        /* Form  y := alpha*A'*x + y. */
 
454
        jy = ky;
 
455
        if (incx == 1) {
 
456
            for (j = 0; j < A->ncol; ++j) {
 
457
                temp = 0.;
 
458
                for (i = Astore->colptr[j]; i < Astore->colptr[j+1]; ++i) {
 
459
                    irow = Astore->rowind[i];
 
460
                    temp += Aval[i] * x[irow];
 
461
                }
 
462
                y[jy] += alpha * temp;
 
463
                jy += incy;
 
464
            }
 
465
        } else {
 
466
            ABORT("Not implemented.");
 
467
        }
 
468
    }
 
469
    return 0;
 
470
} /* sp_sgemv */
 
471
 
 
472
 
 
473