← Back to branch summary

~ubuntu-branches/ubuntu/gutsy/blender/gutsy-security

« back to all changes in this revision

Viewing changes to intern/opennl/superlu/util.c

  • Committer: Bazaar Package Importer
  • Author(s): Florian Ernst
  • Date: 2005-11-06 12:40:03 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20051106124003-3pgs7tcg5rox96xg
Tags: 2.37a-1.1
http://bugs.debian.org/333958
* Non-maintainer upload.
* Split out parts of 01_SConstruct_debian.dpatch again: root_build_dir
  really needs to get adjusted before the clean target runs - closes: #333958,
  see #288882 for reference

Show diffs side-by-side

added added

removed removed

Lines of Context:
intern/opennl/superlu/util.c
 
1
/*
 
2
 * -- SuperLU routine (version 3.0) --
 
3
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 
4
 * and Lawrence Berkeley National Lab.
 
5
 * October 15, 2003
 
6
 *
 
7
 */
 
8
/*
 
9
  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
 
10
 
 
11
  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
 
12
  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
 
13
 
 
14
  Permission is hereby granted to use or copy this program for any
 
15
  purpose, provided the above notices are retained on all copies.
 
16
  Permission to modify the code and to distribute modified code is
 
17
  granted, provided the above notices are retained, and a notice that
 
18
  the code was modified is included with the above copyright notice.
 
19
*/
 
20
 
 
21
#include <math.h>
 
22
#include "ssp_defs.h"
 
23
#include "util.h"
 
24
 
 
25
/* prototypes */
 
26
flops_t LUFactFlops(SuperLUStat_t *stat);
 
27
flops_t LUSolveFlops(SuperLUStat_t *stat);
 
28
float SpaSize(int n, int np, float sum_npw);
 
29
float DenseSize(int n, float sum_nw);
 
30
 
 
31
/* 
 
32
 * Global statistics variale
 
33
 */
 
34
 
 
35
void superlu_abort_and_exit(char* msg)
 
36
{
 
37
    fprintf(stderr, msg);
 
38
    exit (-1);
 
39
}
 
40
 
 
41
/*
 
42
 * Set the default values for the options argument.
 
43
 */
 
44
void set_default_options(superlu_options_t *options)
 
45
{
 
46
    options->Fact = DOFACT;
 
47
    options->Equil = YES;
 
48
    options->ColPerm = COLAMD;
 
49
    options->DiagPivotThresh = 1.0;
 
50
    options->Trans = NOTRANS;
 
51
    options->IterRefine = NOREFINE;
 
52
    options->SymmetricMode = NO;
 
53
    options->PivotGrowth = NO;
 
54
    options->ConditionNumber = NO;
 
55
    options->PrintStat = YES;
 
56
}
 
57
 
 
58
/* Deallocate the structure pointing to the actual storage of the matrix. */
 
59
void
 
60
Destroy_SuperMatrix_Store(SuperMatrix *A)
 
61
{
 
62
    SUPERLU_FREE ( A->Store );
 
63
}
 
64
 
 
65
void
 
66
Destroy_CompCol_Matrix(SuperMatrix *A)
 
67
{
 
68
    SUPERLU_FREE( ((NCformat *)A->Store)->rowind );
 
69
    SUPERLU_FREE( ((NCformat *)A->Store)->colptr );
 
70
    SUPERLU_FREE( ((NCformat *)A->Store)->nzval );
 
71
    SUPERLU_FREE( A->Store );
 
72
}
 
73
 
 
74
void
 
75
Destroy_CompRow_Matrix(SuperMatrix *A)
 
76
{
 
77
    SUPERLU_FREE( ((NRformat *)A->Store)->colind );
 
78
    SUPERLU_FREE( ((NRformat *)A->Store)->rowptr );
 
79
    SUPERLU_FREE( ((NRformat *)A->Store)->nzval );
 
80
    SUPERLU_FREE( A->Store );
 
81
}
 
82
 
 
83
void
 
84
Destroy_SuperNode_Matrix(SuperMatrix *A)
 
85
{
 
86
    SUPERLU_FREE ( ((SCformat *)A->Store)->rowind );
 
87
    SUPERLU_FREE ( ((SCformat *)A->Store)->rowind_colptr );
 
88
    SUPERLU_FREE ( ((SCformat *)A->Store)->nzval );
 
89
    SUPERLU_FREE ( ((SCformat *)A->Store)->nzval_colptr );
 
90
    SUPERLU_FREE ( ((SCformat *)A->Store)->col_to_sup );
 
91
    SUPERLU_FREE ( ((SCformat *)A->Store)->sup_to_col );
 
92
    SUPERLU_FREE ( A->Store );
 
93
}
 
94
 
 
95
/* A is of type Stype==NCP */
 
96
void
 
97
Destroy_CompCol_Permuted(SuperMatrix *A)
 
98
{
 
99
    SUPERLU_FREE ( ((NCPformat *)A->Store)->colbeg );
 
100
    SUPERLU_FREE ( ((NCPformat *)A->Store)->colend );
 
101
    SUPERLU_FREE ( A->Store );
 
102
}
 
103
 
 
104
/* A is of type Stype==DN */
 
105
void
 
106
Destroy_Dense_Matrix(SuperMatrix *A)
 
107
{
 
108
    DNformat* Astore = A->Store;
 
109
    SUPERLU_FREE (Astore->nzval);
 
110
    SUPERLU_FREE ( A->Store );
 
111
}
 
112
 
 
113
/*
 
114
 * Reset repfnz[] for the current column 
 
115
 */
 
116
void
 
117
resetrep_col (const int nseg, const int *segrep, int *repfnz)
 
118
{
 
119
    int i, irep;
 
120
    
 
121
    for (i = 0; i < nseg; i++) {
 
122
        irep = segrep[i];
 
123
        repfnz[irep] = EMPTY;
 
124
    }
 
125
}
 
126
 
 
127
 
 
128
/*
 
129
 * Count the total number of nonzeros in factors L and U,  and in the 
 
130
 * symmetrically reduced L. 
 
131
 */
 
132
void
 
133
countnz(const int n, int *xprune, int *nnzL, int *nnzU, GlobalLU_t *Glu)
 
134
{
 
135
    int          nsuper, fsupc, i, j;
 
136
    int          nnzL0, jlen, irep;
 
137
    int          *xsup, *xlsub;
 
138
 
 
139
    xsup   = Glu->xsup;
 
140
    xlsub  = Glu->xlsub;
 
141
    *nnzL  = 0;
 
142
    *nnzU  = (Glu->xusub)[n];
 
143
    nnzL0  = 0;
 
144
    nsuper = (Glu->supno)[n];
 
145
 
 
146
    if ( n <= 0 ) return;
 
147
 
 
148
    /* 
 
149
     * For each supernode
 
150
     */
 
151
    for (i = 0; i <= nsuper; i++) {
 
152
        fsupc = xsup[i];
 
153
        jlen = xlsub[fsupc+1] - xlsub[fsupc];
 
154
 
 
155
        for (j = fsupc; j < xsup[i+1]; j++) {
 
156
            *nnzL += jlen;
 
157
            *nnzU += j - fsupc + 1;
 
158
            jlen--;
 
159
        }
 
160
        irep = xsup[i+1] - 1;
 
161
        nnzL0 += xprune[irep] - xlsub[irep];
 
162
    }
 
163
    
 
164
    /* printf("\tNo of nonzeros in symm-reduced L = %d\n", nnzL0);*/
 
165
}
 
166
 
 
167
 
 
168
 
 
169
/*
 
170
 * Fix up the data storage lsub for L-subscripts. It removes the subscript
 
171
 * sets for structural pruning, and applies permuation to the remaining
 
172
 * subscripts.
 
173
 */
 
174
void
 
175
fixupL(const int n, const int *perm_r, GlobalLU_t *Glu)
 
176
{
 
177
    register int nsuper, fsupc, nextl, i, j, k, jstrt;
 
178
    int          *xsup, *lsub, *xlsub;
 
179
 
 
180
    if ( n <= 1 ) return;
 
181
 
 
182
    xsup   = Glu->xsup;
 
183
    lsub   = Glu->lsub;
 
184
    xlsub  = Glu->xlsub;
 
185
    nextl  = 0;
 
186
    nsuper = (Glu->supno)[n];
 
187
    
 
188
    /* 
 
189
     * For each supernode ...
 
190
     */
 
191
    for (i = 0; i <= nsuper; i++) {
 
192
        fsupc = xsup[i];
 
193
        jstrt = xlsub[fsupc];
 
194
        xlsub[fsupc] = nextl;
 
195
        for (j = jstrt; j < xlsub[fsupc+1]; j++) {
 
196
            lsub[nextl] = perm_r[lsub[j]]; /* Now indexed into P*A */
 
197
            nextl++;
 
198
        }
 
199
        for (k = fsupc+1; k < xsup[i+1]; k++) 
 
200
                xlsub[k] = nextl;       /* Other columns in supernode i */
 
201
 
 
202
    }
 
203
 
 
204
    xlsub[n] = nextl;
 
205
}
 
206
 
 
207
 
 
208
/*
 
209
 * Diagnostic print of segment info after panel_dfs().
 
210
 */
 
211
void print_panel_seg(int n, int w, int jcol, int nseg, 
 
212
                     int *segrep, int *repfnz)
 
213
{
 
214
    int j, k;
 
215
    
 
216
    for (j = jcol; j < jcol+w; j++) {
 
217
        printf("\tcol %d:\n", j);
 
218
        for (k = 0; k < nseg; k++)
 
219
            printf("\t\tseg %d, segrep %d, repfnz %d\n", k, 
 
220
                        segrep[k], repfnz[(j-jcol)*n + segrep[k]]);
 
221
    }
 
222
 
 
223
}
 
224
 
 
225
 
 
226
void
 
227
StatInit(SuperLUStat_t *stat)
 
228
{
 
229
    register int i, w, panel_size, relax;
 
230
 
 
231
    panel_size = sp_ienv(1);
 
232
    relax = sp_ienv(2);
 
233
    w = SUPERLU_MAX(panel_size, relax);
 
234
    stat->panel_histo = intCalloc(w+1);
 
235
    stat->utime = (double *) SUPERLU_MALLOC(NPHASES * sizeof(double));
 
236
    if (!stat->utime) ABORT("SUPERLU_MALLOC fails for stat->utime");
 
237
    stat->ops = (flops_t *) SUPERLU_MALLOC(NPHASES * sizeof(flops_t));
 
238
    if (!stat->ops) ABORT("SUPERLU_MALLOC fails for stat->ops");
 
239
    for (i = 0; i < NPHASES; ++i) {
 
240
        stat->utime[i] = 0.;
 
241
        stat->ops[i] = 0.;
 
242
    }
 
243
}
 
244
 
 
245
 
 
246
void
 
247
StatPrint(SuperLUStat_t *stat)
 
248
{
 
249
    double         *utime;
 
250
    flops_t        *ops;
 
251
 
 
252
    utime = stat->utime;
 
253
    ops   = stat->ops;
 
254
    printf("Factor time  = %8.2f\n", utime[FACT]);
 
255
    if ( utime[FACT] != 0.0 )
 
256
      printf("Factor flops = %e\tMflops = %8.2f\n", ops[FACT],
 
257
             ops[FACT]*1e-6/utime[FACT]);
 
258
 
 
259
    printf("Solve time   = %8.2f\n", utime[SOLVE]);
 
260
    if ( utime[SOLVE] != 0.0 )
 
261
      printf("Solve flops = %e\tMflops = %8.2f\n", ops[SOLVE],
 
262
             ops[SOLVE]*1e-6/utime[SOLVE]);
 
263
 
 
264
}
 
265
 
 
266
 
 
267
void
 
268
StatFree(SuperLUStat_t *stat)
 
269
{
 
270
    SUPERLU_FREE(stat->panel_histo);
 
271
    SUPERLU_FREE(stat->utime);
 
272
    SUPERLU_FREE(stat->ops);
 
273
}
 
274
 
 
275
 
 
276
flops_t
 
277
LUFactFlops(SuperLUStat_t *stat)
 
278
{
 
279
    return (stat->ops[FACT]);
 
280
}
 
281
 
 
282
flops_t
 
283
LUSolveFlops(SuperLUStat_t *stat)
 
284
{
 
285
    return (stat->ops[SOLVE]);
 
286
}
 
287
 
 
288
 
 
289
 
 
290
 
 
291
 
 
292
/* 
 
293
 * Fills an integer array with a given value.
 
294
 */
 
295
void ifill(int *a, int alen, int ival)
 
296
{
 
297
    register int i;
 
298
    for (i = 0; i < alen; i++) a[i] = ival;
 
299
}
 
300
 
 
301
 
 
302
 
 
303
/* 
 
304
 * Get the statistics of the supernodes 
 
305
 */
 
306
#define NBUCKS 10
 
307
static  int     max_sup_size;
 
308
 
 
309
void super_stats(int nsuper, int *xsup)
 
310
{
 
311
    register int nsup1 = 0;
 
312
    int          i, isize, whichb, bl, bh;
 
313
    int          bucket[NBUCKS];
 
314
 
 
315
    max_sup_size = 0;
 
316
 
 
317
    for (i = 0; i <= nsuper; i++) {
 
318
        isize = xsup[i+1] - xsup[i];
 
319
        if ( isize == 1 ) nsup1++;
 
320
        if ( max_sup_size < isize ) max_sup_size = isize;       
 
321
    }
 
322
 
 
323
    printf("    Supernode statistics:\n\tno of super = %d\n", nsuper+1);
 
324
    printf("\tmax supernode size = %d\n", max_sup_size);
 
325
    printf("\tno of size 1 supernodes = %d\n", nsup1);
 
326
 
 
327
    /* Histogram of the supernode sizes */
 
328
    ifill (bucket, NBUCKS, 0);
 
329
 
 
330
    for (i = 0; i <= nsuper; i++) {
 
331
        isize = xsup[i+1] - xsup[i];
 
332
        whichb = (float) isize / max_sup_size * NBUCKS;
 
333
        if (whichb >= NBUCKS) whichb = NBUCKS - 1;
 
334
        bucket[whichb]++;
 
335
    }
 
336
    
 
337
    printf("\tHistogram of supernode sizes:\n");
 
338
    for (i = 0; i < NBUCKS; i++) {
 
339
        bl = (float) i * max_sup_size / NBUCKS;
 
340
        bh = (float) (i+1) * max_sup_size / NBUCKS;
 
341
        printf("\tsnode: %d-%d\t\t%d\n", bl+1, bh, bucket[i]);
 
342
    }
 
343
 
 
344
}
 
345
 
 
346
 
 
347
float SpaSize(int n, int np, float sum_npw)
 
348
{
 
349
    return (sum_npw*8 + np*8 + n*4)/1024.;
 
350
}
 
351
 
 
352
float DenseSize(int n, float sum_nw)
 
353
{
 
354
    return (sum_nw*8 + n*8)/1024.;;
 
355
}
 
356
 
 
357
 
 
358
 
 
359
/*
 
360
 * Check whether repfnz[] == EMPTY after reset.
 
361
 */
 
362
void check_repfnz(int n, int w, int jcol, int *repfnz)
 
363
{
 
364
    int jj, k;
 
365
 
 
366
    for (jj = jcol; jj < jcol+w; jj++) 
 
367
        for (k = 0; k < n; k++)
 
368
            if ( repfnz[(jj-jcol)*n + k] != EMPTY ) {
 
369
                fprintf(stderr, "col %d, repfnz_col[%d] = %d\n", jj,
 
370
                        k, repfnz[(jj-jcol)*n + k]);
 
371
                ABORT("check_repfnz");
 
372
            }
 
373
}
 
374
 
 
375
 
 
376
/* Print a summary of the testing results. */
 
377
void
 
378
PrintSumm(char *type, int nfail, int nrun, int nerrs)
 
379
{
 
380
    if ( nfail > 0 )
 
381
        printf("%3s driver: %d out of %d tests failed to pass the threshold\n",
 
382
               type, nfail, nrun);
 
383
    else
 
384
        printf("All tests for %3s driver passed the threshold (%6d tests run)\n", type, nrun);
 
385
 
 
386
    if ( nerrs > 0 )
 
387
        printf("%6d error messages recorded\n", nerrs);
 
388
}
 
389
 
 
390
 
 
391
int print_int_vec(char *what, int n, int *vec)
 
392
{
 
393
    int i;
 
394
    printf("%s\n", what);
 
395
    for (i = 0; i < n; ++i) printf("%d\t%d\n", i, vec[i]);
 
396
    return 0;
 
397
}