← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/linsolve/SuperLU/SRC/cutil.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/linsolve/SuperLU/SRC/cutil.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 <math.h>
23
 
#include "csp_defs.h"
24
 
 
25
 
void
26
 
cCreate_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz, 
27
 
                       complex *nzval, int *rowind, int *colptr,
28
 
                       Stype_t stype, Dtype_t dtype, Mtype_t mtype)
29
 
{
30
 
    NCformat *Astore;
31
 
 
32
 
    A->Stype = stype;
33
 
    A->Dtype = dtype;
34
 
    A->Mtype = mtype;
35
 
    A->nrow = m;
36
 
    A->ncol = n;
37
 
    A->Store = (void *) SUPERLU_MALLOC( sizeof(NCformat) );
38
 
    if ( !(A->Store) ) ABORT("SUPERLU_MALLOC fails for A->Store");
39
 
    Astore = A->Store;
40
 
    Astore->nnz = nnz;
41
 
    Astore->nzval = nzval;
42
 
    Astore->rowind = rowind;
43
 
    Astore->colptr = colptr;
44
 
}
45
 
 
46
 
void
47
 
cCreate_CompRow_Matrix(SuperMatrix *A, int m, int n, int nnz, 
48
 
                       complex *nzval, int *colind, int *rowptr,
49
 
                       Stype_t stype, Dtype_t dtype, Mtype_t mtype)
50
 
{
51
 
    NRformat *Astore;
52
 
 
53
 
    A->Stype = stype;
54
 
    A->Dtype = dtype;
55
 
    A->Mtype = mtype;
56
 
    A->nrow = m;
57
 
    A->ncol = n;
58
 
    A->Store = (void *) SUPERLU_MALLOC( sizeof(NRformat) );
59
 
    if ( !(A->Store) ) ABORT("SUPERLU_MALLOC fails for A->Store");
60
 
    Astore = A->Store;
61
 
    Astore->nnz = nnz;
62
 
    Astore->nzval = nzval;
63
 
    Astore->colind = colind;
64
 
    Astore->rowptr = rowptr;
65
 
}
66
 
 
67
 
/* Copy matrix A into matrix B. */
68
 
void
69
 
cCopy_CompCol_Matrix(SuperMatrix *A, SuperMatrix *B)
70
 
{
71
 
    NCformat *Astore, *Bstore;
72
 
    int      ncol, nnz, i;
73
 
 
74
 
    B->Stype = A->Stype;
75
 
    B->Dtype = A->Dtype;
76
 
    B->Mtype = A->Mtype;
77
 
    B->nrow  = A->nrow;;
78
 
    B->ncol  = ncol = A->ncol;
79
 
    Astore   = (NCformat *) A->Store;
80
 
    Bstore   = (NCformat *) B->Store;
81
 
    Bstore->nnz = nnz = Astore->nnz;
82
 
    for (i = 0; i < nnz; ++i)
83
 
        ((complex *)Bstore->nzval)[i] = ((complex *)Astore->nzval)[i];
84
 
    for (i = 0; i < nnz; ++i) Bstore->rowind[i] = Astore->rowind[i];
85
 
    for (i = 0; i <= ncol; ++i) Bstore->colptr[i] = Astore->colptr[i];
86
 
}
87
 
 
88
 
 
89
 
void
90
 
cCreate_Dense_Matrix(SuperMatrix *X, int m, int n, complex *x, int ldx,
91
 
                    Stype_t stype, Dtype_t dtype, Mtype_t mtype)
92
 
{
93
 
    DNformat    *Xstore;
94
 
    
95
 
    X->Stype = stype;
96
 
    X->Dtype = dtype;
97
 
    X->Mtype = mtype;
98
 
    X->nrow = m;
99
 
    X->ncol = n;
100
 
    X->Store = (void *) SUPERLU_MALLOC( sizeof(DNformat) );
101
 
    if ( !(X->Store) ) ABORT("SUPERLU_MALLOC fails for X->Store");
102
 
    Xstore = (DNformat *) X->Store;
103
 
    Xstore->lda = ldx;
104
 
    Xstore->nzval = (complex *) x;
105
 
}
106
 
 
107
 
void
108
 
cCopy_Dense_Matrix(int M, int N, complex *X, int ldx,
109
 
                        complex *Y, int ldy)
110
 
{
111
 
/*
112
 
 *
113
 
 *  Purpose
114
 
 *  =======
115
 
 *
116
 
 *  Copies a two-dimensional matrix X to another matrix Y.
117
 
 */
118
 
    int    i, j;
119
 
    
120
 
    for (j = 0; j < N; ++j)
121
 
        for (i = 0; i < M; ++i)
122
 
            Y[i + j*ldy] = X[i + j*ldx];
123
 
}
124
 
 
125
 
void
126
 
cCreate_SuperNode_Matrix(SuperMatrix *L, int m, int n, int nnz, 
127
 
                        complex *nzval, int *nzval_colptr, int *rowind,
128
 
                        int *rowind_colptr, int *col_to_sup, int *sup_to_col,
129
 
                        Stype_t stype, Dtype_t dtype, Mtype_t mtype)
130
 
{
131
 
    SCformat *Lstore;
132
 
 
133
 
    L->Stype = stype;
134
 
    L->Dtype = dtype;
135
 
    L->Mtype = mtype;
136
 
    L->nrow = m;
137
 
    L->ncol = n;
138
 
    L->Store = (void *) SUPERLU_MALLOC( sizeof(SCformat) );
139
 
    if ( !(L->Store) ) ABORT("SUPERLU_MALLOC fails for L->Store");
140
 
    Lstore = L->Store;
141
 
    Lstore->nnz = nnz;
142
 
    Lstore->nsuper = col_to_sup[n];
143
 
    Lstore->nzval = nzval;
144
 
    Lstore->nzval_colptr = nzval_colptr;
145
 
    Lstore->rowind = rowind;
146
 
    Lstore->rowind_colptr = rowind_colptr;
147
 
    Lstore->col_to_sup = col_to_sup;
148
 
    Lstore->sup_to_col = sup_to_col;
149
 
 
150
 
}
151
 
 
152
 
 
153
 
/*
154
 
 * Convert a row compressed storage into a column compressed storage.
155
 
 */
156
 
void
157
 
cCompRow_to_CompCol(int m, int n, int nnz, 
158
 
                    complex *a, int *colind, int *rowptr,
159
 
                    complex **at, int **rowind, int **colptr)
160
 
{
161
 
    register int i, j, col, relpos;
162
 
    int *marker;
163
 
 
164
 
    /* Allocate storage for another copy of the matrix. */
165
 
    *at = (complex *) complexMalloc(nnz);
166
 
    *rowind = (int *) intMalloc(nnz);
167
 
    *colptr = (int *) intMalloc(n+1);
168
 
    marker = (int *) intCalloc(n);
169
 
    
170
 
    /* Get counts of each column of A, and set up column pointers */
171
 
    for (i = 0; i < m; ++i)
172
 
        for (j = rowptr[i]; j < rowptr[i+1]; ++j) ++marker[colind[j]];
173
 
    (*colptr)[0] = 0;
174
 
    for (j = 0; j < n; ++j) {
175
 
        (*colptr)[j+1] = (*colptr)[j] + marker[j];
176
 
        marker[j] = (*colptr)[j];
177
 
    }
178
 
 
179
 
    /* Transfer the matrix into the compressed column storage. */
180
 
    for (i = 0; i < m; ++i) {
181
 
        for (j = rowptr[i]; j < rowptr[i+1]; ++j) {
182
 
            col = colind[j];
183
 
            relpos = marker[col];
184
 
            (*rowind)[relpos] = i;
185
 
            (*at)[relpos] = a[j];
186
 
            ++marker[col];
187
 
        }
188
 
    }
189
 
 
190
 
    SUPERLU_FREE(marker);
191
 
}
192
 
 
193
 
 
194
 
void
195
 
cPrint_CompCol_Matrix(char *what, SuperMatrix *A)
196
 
{
197
 
    NCformat     *Astore;
198
 
    register int i,n;
199
 
    float       *dp;
200
 
    
201
 
    printf("\nCompCol matrix %s:\n", what);
202
 
    printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
203
 
    n = A->ncol;
204
 
    Astore = (NCformat *) A->Store;
205
 
    dp = (float *) Astore->nzval;
206
 
    printf("nrow %d, ncol %d, nnz %d\n", A->nrow,A->ncol,Astore->nnz);
207
 
    printf("nzval: ");
208
 
    for (i = 0; i < 2*Astore->colptr[n]; ++i) printf("%f  ", dp[i]);
209
 
    printf("\nrowind: ");
210
 
    for (i = 0; i < Astore->colptr[n]; ++i) printf("%d  ", Astore->rowind[i]);
211
 
    printf("\ncolptr: ");
212
 
    for (i = 0; i <= n; ++i) printf("%d  ", Astore->colptr[i]);
213
 
    printf("\n");
214
 
    fflush(stdout);
215
 
}
216
 
 
217
 
void
218
 
cPrint_SuperNode_Matrix(char *what, SuperMatrix *A)
219
 
{
220
 
    SCformat     *Astore;
221
 
    register int i, j, k, c, d, n, nsup;
222
 
    float       *dp;
223
 
    int *col_to_sup, *sup_to_col, *rowind, *rowind_colptr;
224
 
    
225
 
    printf("\nSuperNode matrix %s:\n", what);
226
 
    printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
227
 
    n = A->ncol;
228
 
    Astore = (SCformat *) A->Store;
229
 
    dp = (float *) Astore->nzval;
230
 
    col_to_sup = Astore->col_to_sup;
231
 
    sup_to_col = Astore->sup_to_col;
232
 
    rowind_colptr = Astore->rowind_colptr;
233
 
    rowind = Astore->rowind;
234
 
    printf("nrow %d, ncol %d, nnz %d, nsuper %d\n", 
235
 
           A->nrow,A->ncol,Astore->nnz,Astore->nsuper);
236
 
    printf("nzval:\n");
237
 
    for (k = 0; k <= Astore->nsuper; ++k) {
238
 
      c = sup_to_col[k];
239
 
      nsup = sup_to_col[k+1] - c;
240
 
      for (j = c; j < c + nsup; ++j) {
241
 
        d = Astore->nzval_colptr[j];
242
 
        for (i = rowind_colptr[c]; i < rowind_colptr[c+1]; ++i) {
243
 
          printf("%d\t%d\t%e\t%e\n", rowind[i], j, dp[d++], dp[d++]);
244
 
        }
245
 
      }
246
 
    }
247
 
#if 0
248
 
    for (i = 0; i < 2*Astore->nzval_colptr[n]; ++i) printf("%f  ", dp[i]);
249
 
#endif
250
 
    printf("\nnzval_colptr: ");
251
 
    for (i = 0; i <= n; ++i) printf("%d  ", Astore->nzval_colptr[i]);
252
 
    printf("\nrowind: ");
253
 
    for (i = 0; i < Astore->rowind_colptr[n]; ++i) 
254
 
        printf("%d  ", Astore->rowind[i]);
255
 
    printf("\nrowind_colptr: ");
256
 
    for (i = 0; i <= n; ++i) printf("%d  ", Astore->rowind_colptr[i]);
257
 
    printf("\ncol_to_sup: ");
258
 
    for (i = 0; i < n; ++i) printf("%d  ", col_to_sup[i]);
259
 
    printf("\nsup_to_col: ");
260
 
    for (i = 0; i <= Astore->nsuper+1; ++i) 
261
 
        printf("%d  ", sup_to_col[i]);
262
 
    printf("\n");
263
 
    fflush(stdout);
264
 
}
265
 
 
266
 
void
267
 
cPrint_Dense_Matrix(char *what, SuperMatrix *A)
268
 
{
269
 
    DNformat     *Astore;
270
 
    register int i;
271
 
    float       *dp;
272
 
    
273
 
    printf("\nDense matrix %s:\n", what);
274
 
    printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
275
 
    Astore = (DNformat *) A->Store;
276
 
    dp = (float *) Astore->nzval;
277
 
    printf("nrow %d, ncol %d, lda %d\n", A->nrow,A->ncol,Astore->lda);
278
 
    printf("\nnzval: ");
279
 
    for (i = 0; i < 2*A->nrow; ++i) printf("%f  ", dp[i]);
280
 
    printf("\n");
281
 
    fflush(stdout);
282
 
}
283
 
 
284
 
/*
285
 
 * Diagnostic print of column "jcol" in the U/L factor.
286
 
 */
287
 
void
288
 
cprint_lu_col(char *msg, int jcol, int pivrow, int *xprune, GlobalLU_t *Glu)
289
 
{
290
 
    int     i, k, fsupc;
291
 
    int     *xsup, *supno;
292
 
    int     *xlsub, *lsub;
293
 
    complex  *lusup;
294
 
    int     *xlusup;
295
 
    complex  *ucol;
296
 
    int     *usub, *xusub;
297
 
 
298
 
    xsup    = Glu->xsup;
299
 
    supno   = Glu->supno;
300
 
    lsub    = Glu->lsub;
301
 
    xlsub   = Glu->xlsub;
302
 
    lusup   = Glu->lusup;
303
 
    xlusup  = Glu->xlusup;
304
 
    ucol    = Glu->ucol;
305
 
    usub    = Glu->usub;
306
 
    xusub   = Glu->xusub;
307
 
    
308
 
    printf("%s", msg);
309
 
    printf("col %d: pivrow %d, supno %d, xprune %d\n", 
310
 
           jcol, pivrow, supno[jcol], xprune[jcol]);
311
 
    
312
 
    printf("\tU-col:\n");
313
 
    for (i = xusub[jcol]; i < xusub[jcol+1]; i++)
314
 
        printf("\t%d%10.4f, %10.4f\n", usub[i], ucol[i].r, ucol[i].i);
315
 
    printf("\tL-col in rectangular snode:\n");
316
 
    fsupc = xsup[supno[jcol]];  /* first col of the snode */
317
 
    i = xlsub[fsupc];
318
 
    k = xlusup[jcol];
319
 
    while ( i < xlsub[fsupc+1] && k < xlusup[jcol+1] ) {
320
 
        printf("\t%d\t%10.4f, %10.4f\n", lsub[i], lusup[k].r, lusup[k].i);
321
 
        i++; k++;
322
 
    }
323
 
    fflush(stdout);
324
 
}
325
 
 
326
 
 
327
 
/*
328
 
 * Check whether tempv[] == 0. This should be true before and after 
329
 
 * calling any numeric routines, i.e., "panel_bmod" and "column_bmod". 
330
 
 */
331
 
void ccheck_tempv(int n, complex *tempv)
332
 
{
333
 
    int i;
334
 
        
335
 
    for (i = 0; i < n; i++) {
336
 
        if ((tempv[i].r != 0.0) || (tempv[i].i != 0.0))
337
 
        {
338
 
            fprintf(stderr,"tempv[%d] = {%f, %f}\n", i, tempv[i].r, tempv[i].i);
339
 
            ABORT("ccheck_tempv");
340
 
        }
341
 
    }
342
 
}
343
 
 
344
 
 
345
 
void
346
 
cGenXtrue(int n, int nrhs, complex *x, int ldx)
347
 
{
348
 
    int  i, j;
349
 
    for (j = 0; j < nrhs; ++j)
350
 
        for (i = 0; i < n; ++i) {
351
 
            x[i + j*ldx].r = 1.0;
352
 
            x[i + j*ldx].i = 0.0;
353
 
        }
354
 
}
355
 
 
356
 
/*
357
 
 * Let rhs[i] = sum of i-th row of A, so the solution vector is all 1's
358
 
 */
359
 
void
360
 
cFillRHS(trans_t trans, int nrhs, complex *x, int ldx,
361
 
         SuperMatrix *A, SuperMatrix *B)
362
 
{
363
 
    NCformat *Astore;
364
 
    complex   *Aval;
365
 
    DNformat *Bstore;
366
 
    complex   *rhs;
367
 
    complex one = {1.0, 0.0};
368
 
    complex zero = {0.0, 0.0};
369
 
    int      ldc;
370
 
    char transc[1];
371
 
 
372
 
    Astore = A->Store;
373
 
    Aval   = (complex *) Astore->nzval;
374
 
    Bstore = B->Store;
375
 
    rhs    = Bstore->nzval;
376
 
    ldc    = Bstore->lda;
377
 
    
378
 
    if ( trans == NOTRANS ) *(unsigned char *)transc = 'N';
379
 
    else *(unsigned char *)transc = 'T';
380
 
 
381
 
    sp_cgemm(transc, "N", A->nrow, nrhs, A->ncol, one, A,
382
 
             x, ldx, zero, rhs, ldc);
383
 
 
384
 
}
385
 
 
386
 
/* 
387
 
 * Fills a complex precision array with a given value.
388
 
 */
389
 
void 
390
 
cfill(complex *a, int alen, complex dval)
391
 
{
392
 
    register int i;
393
 
    for (i = 0; i < alen; i++) a[i] = dval;
394
 
}
395
 
 
396
 
 
397
 
 
398
 
/* 
399
 
 * Check the inf-norm of the error vector 
400
 
 */
401
 
void cinf_norm_error(int nrhs, SuperMatrix *X, complex *xtrue)
402
 
{
403
 
    DNformat *Xstore;
404
 
    float err, xnorm;
405
 
    complex *Xmat, *soln_work;
406
 
    complex temp;
407
 
    int i, j;
408
 
 
409
 
    Xstore = X->Store;
410
 
    Xmat = Xstore->nzval;
411
 
 
412
 
    for (j = 0; j < nrhs; j++) {
413
 
      soln_work = &Xmat[j*Xstore->lda];
414
 
      err = xnorm = 0.0;
415
 
      for (i = 0; i < X->nrow; i++) {
416
 
        c_sub(&temp, &soln_work[i], &xtrue[i]);
417
 
        err = SUPERLU_MAX(err, c_abs(&temp));
418
 
        xnorm = SUPERLU_MAX(xnorm, c_abs(&soln_work[i]));
419
 
      }
420
 
      err = err / xnorm;
421
 
      printf("||X - Xtrue||/||X|| = %e\n", err);
422
 
    }
423
 
}
424
 
 
425
 
 
426
 
 
427
 
/* Print performance of the code. */
428
 
void
429
 
cPrintPerf(SuperMatrix *L, SuperMatrix *U, mem_usage_t *mem_usage,
430
 
           float rpg, float rcond, float *ferr,
431
 
           float *berr, char *equed, SuperLUStat_t *stat)
432
 
{
433
 
    SCformat *Lstore;
434
 
    NCformat *Ustore;
435
 
    double   *utime;
436
 
    flops_t  *ops;
437
 
    
438
 
    utime = stat->utime;
439
 
    ops   = stat->ops;
440
 
    
441
 
    if ( utime[FACT] != 0. )
442
 
        printf("Factor flops = %e\tMflops = %8.2f\n", ops[FACT],
443
 
               ops[FACT]*1e-6/utime[FACT]);
444
 
    printf("Identify relaxed snodes     = %8.2f\n", utime[RELAX]);
445
 
    if ( utime[SOLVE] != 0. )
446
 
        printf("Solve flops = %.0f, Mflops = %8.2f\n", ops[SOLVE],
447
 
               ops[SOLVE]*1e-6/utime[SOLVE]);
448
 
    
449
 
    Lstore = (SCformat *) L->Store;
450
 
    Ustore = (NCformat *) U->Store;
451
 
    printf("\tNo of nonzeros in factor L = %d\n", Lstore->nnz);
452
 
    printf("\tNo of nonzeros in factor U = %d\n", Ustore->nnz);
453
 
    printf("\tNo of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz);
454
 
        
455
 
    printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
456
 
           mem_usage->for_lu/1e6, mem_usage->total_needed/1e6,
457
 
           mem_usage->expansions);
458
 
        
459
 
    printf("\tFactor\tMflops\tSolve\tMflops\tEtree\tEquil\tRcond\tRefine\n");
460
 
    printf("PERF:%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f\n",
461
 
           utime[FACT], ops[FACT]*1e-6/utime[FACT],
462
 
           utime[SOLVE], ops[SOLVE]*1e-6/utime[SOLVE],
463
 
           utime[ETREE], utime[EQUIL], utime[RCOND], utime[REFINE]);
464
 
    
465
 
    printf("\tRpg\t\tRcond\t\tFerr\t\tBerr\t\tEquil?\n");
466
 
    printf("NUM:\t%e\t%e\t%e\t%e\t%s\n",
467
 
           rpg, rcond, ferr[0], berr[0], equed);
468
 
    
469
 
}
470
 
 
471
 
 
472
 
 
473
 
 
474
 
print_complex_vec(char *what, int n, complex *vec)
475
 
{
476
 
    int i;
477
 
    printf("%s: n %d\n", what, n);
478
 
    for (i = 0; i < n; ++i) printf("%d\t%f%f\n", i, vec[i].r, vec[i].i);
479
 
    return 0;
480
 
}
481