← Back to branch summary

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

« back to all changes in this revision

Viewing changes to intern/opennl/superlu/smemory.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/smemory.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
#include "ssp_defs.h"
 
10
 
 
11
/* Constants */
 
12
#define NO_MEMTYPE  4      /* 0: lusup;
 
13
                              1: ucol;
 
14
                              2: lsub;
 
15
                              3: usub */
 
16
#define GluIntArray(n)   (5 * (n) + 5)
 
17
 
 
18
/* Internal prototypes */
 
19
void  *sexpand (int *, MemType,int, int, GlobalLU_t *);
 
20
int   sLUWorkInit (int, int, int, int **, float **, LU_space_t);
 
21
void  copy_mem_float (int, void *, void *);
 
22
void  sStackCompress (GlobalLU_t *);
 
23
void  sSetupSpace (void *, int, LU_space_t *);
 
24
void  *suser_malloc (int, int);
 
25
void  suser_free (int, int);
 
26
 
 
27
/* External prototypes (in memory.c - prec-indep) */
 
28
extern void    copy_mem_int    (int, void *, void *);
 
29
extern void    user_bcopy      (char *, char *, int);
 
30
 
 
31
/* Headers for 4 types of dynamatically managed memory */
 
32
typedef struct e_node {
 
33
    int size;      /* length of the memory that has been used */
 
34
    void *mem;     /* pointer to the new malloc'd store */
 
35
} ExpHeader;
 
36
 
 
37
typedef struct {
 
38
    int  size;
 
39
    int  used;
 
40
    int  top1;  /* grow upward, relative to &array[0] */
 
41
    int  top2;  /* grow downward */
 
42
    void *array;
 
43
} LU_stack_t;
 
44
 
 
45
/* Variables local to this file */
 
46
static ExpHeader *expanders = 0; /* Array of pointers to 4 types of memory */
 
47
static LU_stack_t stack;
 
48
static int no_expand;
 
49
 
 
50
/* Macros to manipulate stack */
 
51
#define StackFull(x)         ( x + stack.used >= stack.size )
 
52
#define NotDoubleAlign(addr) ( (long int)addr & 7 )
 
53
#define DoubleAlign(addr)    ( ((long int)addr + 7) & ~7L )
 
54
#define TempSpace(m, w)      ( (2*w + 4 + NO_MARKER) * m * sizeof(int) + \
 
55
                              (w + 1) * m * sizeof(float) )
 
56
#define Reduce(alpha)        ((alpha + 1) / 2)  /* i.e. (alpha-1)/2 + 1 */
 
57
 
 
58
 
 
59
 
 
60
 
 
61
/*
 
62
 * Setup the memory model to be used for factorization.
 
63
 *    lwork = 0: use system malloc;
 
64
 *    lwork > 0: use user-supplied work[] space.
 
65
 */
 
66
void sSetupSpace(void *work, int lwork, LU_space_t *MemModel)
 
67
{
 
68
    if ( lwork == 0 ) {
 
69
        *MemModel = SYSTEM; /* malloc/free */
 
70
    } else if ( lwork > 0 ) {
 
71
        *MemModel = USER;   /* user provided space */
 
72
        stack.used = 0;
 
73
        stack.top1 = 0;
 
74
        stack.top2 = (lwork/4)*4; /* must be word addressable */
 
75
        stack.size = stack.top2;
 
76
        stack.array = (void *) work;
 
77
    }
 
78
}
 
79
 
 
80
 
 
81
 
 
82
void *suser_malloc(int bytes, int which_end)
 
83
{
 
84
    void *buf;
 
85
    
 
86
    if ( StackFull(bytes) ) return (NULL);
 
87
 
 
88
    if ( which_end == HEAD ) {
 
89
        buf = (char*) stack.array + stack.top1;
 
90
        stack.top1 += bytes;
 
91
    } else {
 
92
        stack.top2 -= bytes;
 
93
        buf = (char*) stack.array + stack.top2;
 
94
    }
 
95
    
 
96
    stack.used += bytes;
 
97
    return buf;
 
98
}
 
99
 
 
100
 
 
101
void suser_free(int bytes, int which_end)
 
102
{
 
103
    if ( which_end == HEAD ) {
 
104
        stack.top1 -= bytes;
 
105
    } else {
 
106
        stack.top2 += bytes;
 
107
    }
 
108
    stack.used -= bytes;
 
109
}
 
110
 
 
111
 
 
112
 
 
113
/*
 
114
 * mem_usage consists of the following fields:
 
115
 *    - for_lu (float)
 
116
 *      The amount of space used in bytes for the L\U data structures.
 
117
 *    - total_needed (float)
 
118
 *      The amount of space needed in bytes to perform factorization.
 
119
 *    - expansions (int)
 
120
 *      Number of memory expansions during the LU factorization.
 
121
 */
 
122
int sQuerySpace(SuperMatrix *L, SuperMatrix *U, mem_usage_t *mem_usage)
 
123
{
 
124
    SCformat *Lstore;
 
125
    NCformat *Ustore;
 
126
    register int n, iword, dword, panel_size = sp_ienv(1);
 
127
 
 
128
    Lstore = L->Store;
 
129
    Ustore = U->Store;
 
130
    n = L->ncol;
 
131
    iword = sizeof(int);
 
132
    dword = sizeof(float);
 
133
 
 
134
    /* For LU factors */
 
135
    mem_usage->for_lu = (float)( (4*n + 3) * iword + Lstore->nzval_colptr[n] *
 
136
                                 dword + Lstore->rowind_colptr[n] * iword );
 
137
    mem_usage->for_lu += (float)( (n + 1) * iword +
 
138
                                 Ustore->colptr[n] * (dword + iword) );
 
139
 
 
140
    /* Working storage to support factorization */
 
141
    mem_usage->total_needed = mem_usage->for_lu +
 
142
        (float)( (2 * panel_size + 4 + NO_MARKER) * n * iword +
 
143
                (panel_size + 1) * n * dword );
 
144
 
 
145
    mem_usage->expansions = --no_expand;
 
146
 
 
147
    return 0;
 
148
} /* sQuerySpace */
 
149
 
 
150
/*
 
151
 * Allocate storage for the data structures common to all factor routines.
 
152
 * For those unpredictable size, make a guess as FILL * nnz(A).
 
153
 * Return value:
 
154
 *     If lwork = -1, return the estimated amount of space required, plus n;
 
155
 *     otherwise, return the amount of space actually allocated when
 
156
 *     memory allocation failure occurred.
 
157
 */
 
158
int
 
159
sLUMemInit(fact_t fact, void *work, int lwork, int m, int n, int annz,
 
160
          int panel_size, SuperMatrix *L, SuperMatrix *U, GlobalLU_t *Glu,
 
161
          int **iwork, float **dwork)
 
162
{
 
163
    int      info, iword, dword;
 
164
    SCformat *Lstore;
 
165
    NCformat *Ustore;
 
166
    int      *xsup, *supno;
 
167
    int      *lsub, *xlsub;
 
168
    float   *lusup;
 
169
    int      *xlusup;
 
170
    float   *ucol;
 
171
    int      *usub, *xusub;
 
172
    int      nzlmax, nzumax, nzlumax;
 
173
    int      FILL = sp_ienv(6);
 
174
    
 
175
    Glu->n    = n;
 
176
    no_expand = 0;
 
177
    iword     = sizeof(int);
 
178
    dword     = sizeof(float);
 
179
 
 
180
    if ( !expanders )   
 
181
        expanders = (ExpHeader*)SUPERLU_MALLOC(NO_MEMTYPE * sizeof(ExpHeader));
 
182
    if ( !expanders ) ABORT("SUPERLU_MALLOC fails for expanders");
 
183
    
 
184
    if ( fact != SamePattern_SameRowPerm ) {
 
185
        /* Guess for L\U factors */
 
186
        nzumax = nzlumax = FILL * annz;
 
187
        nzlmax = SUPERLU_MAX(1, FILL/4.) * annz;
 
188
 
 
189
        if ( lwork == -1 ) {
 
190
            return ( GluIntArray(n) * iword + TempSpace(m, panel_size)
 
191
                    + (nzlmax+nzumax)*iword + (nzlumax+nzumax)*dword + n );
 
192
        } else {
 
193
            sSetupSpace(work, lwork, &Glu->MemModel);
 
194
        }
 
195
        
 
196
#ifdef DEBUG               
 
197
        printf("sLUMemInit() called: annz %d, MemModel %d\n", 
 
198
                annz, Glu->MemModel);
 
199
#endif  
 
200
        
 
201
        /* Integer pointers for L\U factors */
 
202
        if ( Glu->MemModel == SYSTEM ) {
 
203
            xsup   = intMalloc(n+1);
 
204
            supno  = intMalloc(n+1);
 
205
            xlsub  = intMalloc(n+1);
 
206
            xlusup = intMalloc(n+1);
 
207
            xusub  = intMalloc(n+1);
 
208
        } else {
 
209
            xsup   = (int *)suser_malloc((n+1) * iword, HEAD);
 
210
            supno  = (int *)suser_malloc((n+1) * iword, HEAD);
 
211
            xlsub  = (int *)suser_malloc((n+1) * iword, HEAD);
 
212
            xlusup = (int *)suser_malloc((n+1) * iword, HEAD);
 
213
            xusub  = (int *)suser_malloc((n+1) * iword, HEAD);
 
214
        }
 
215
 
 
216
        lusup = (float *) sexpand( &nzlumax, LUSUP, 0, 0, Glu );
 
217
        ucol  = (float *) sexpand( &nzumax, UCOL, 0, 0, Glu );
 
218
        lsub  = (int *)    sexpand( &nzlmax, LSUB, 0, 0, Glu );
 
219
        usub  = (int *)    sexpand( &nzumax, USUB, 0, 1, Glu );
 
220
 
 
221
        while ( !lusup || !ucol || !lsub || !usub ) {
 
222
            if ( Glu->MemModel == SYSTEM ) {
 
223
                SUPERLU_FREE(lusup); 
 
224
                SUPERLU_FREE(ucol); 
 
225
                SUPERLU_FREE(lsub); 
 
226
                SUPERLU_FREE(usub);
 
227
            } else {
 
228
                suser_free((nzlumax+nzumax)*dword+(nzlmax+nzumax)*iword, HEAD);
 
229
            }
 
230
            nzlumax /= 2;
 
231
            nzumax /= 2;
 
232
            nzlmax /= 2;
 
233
            if ( nzlumax < annz ) {
 
234
                printf("Not enough memory to perform factorization.\n");
 
235
                return (smemory_usage(nzlmax, nzumax, nzlumax, n) + n);
 
236
            }
 
237
            lusup = (float *) sexpand( &nzlumax, LUSUP, 0, 0, Glu );
 
238
            ucol  = (float *) sexpand( &nzumax, UCOL, 0, 0, Glu );
 
239
            lsub  = (int *)    sexpand( &nzlmax, LSUB, 0, 0, Glu );
 
240
            usub  = (int *)    sexpand( &nzumax, USUB, 0, 1, Glu );
 
241
        }
 
242
        
 
243
    } else {
 
244
        /* fact == SamePattern_SameRowPerm */
 
245
        Lstore   = L->Store;
 
246
        Ustore   = U->Store;
 
247
        xsup     = Lstore->sup_to_col;
 
248
        supno    = Lstore->col_to_sup;
 
249
        xlsub    = Lstore->rowind_colptr;
 
250
        xlusup   = Lstore->nzval_colptr;
 
251
        xusub    = Ustore->colptr;
 
252
        nzlmax   = Glu->nzlmax;    /* max from previous factorization */
 
253
        nzumax   = Glu->nzumax;
 
254
        nzlumax  = Glu->nzlumax;
 
255
        
 
256
        if ( lwork == -1 ) {
 
257
            return ( GluIntArray(n) * iword + TempSpace(m, panel_size)
 
258
                    + (nzlmax+nzumax)*iword + (nzlumax+nzumax)*dword + n );
 
259
        } else if ( lwork == 0 ) {
 
260
            Glu->MemModel = SYSTEM;
 
261
        } else {
 
262
            Glu->MemModel = USER;
 
263
            stack.top2 = (lwork/4)*4; /* must be word-addressable */
 
264
            stack.size = stack.top2;
 
265
        }
 
266
        
 
267
        lsub  = expanders[LSUB].mem  = Lstore->rowind;
 
268
        lusup = expanders[LUSUP].mem = Lstore->nzval;
 
269
        usub  = expanders[USUB].mem  = Ustore->rowind;
 
270
        ucol  = expanders[UCOL].mem  = Ustore->nzval;;
 
271
        expanders[LSUB].size         = nzlmax;
 
272
        expanders[LUSUP].size        = nzlumax;
 
273
        expanders[USUB].size         = nzumax;
 
274
        expanders[UCOL].size         = nzumax;  
 
275
    }
 
276
 
 
277
    Glu->xsup    = xsup;
 
278
    Glu->supno   = supno;
 
279
    Glu->lsub    = lsub;
 
280
    Glu->xlsub   = xlsub;
 
281
    Glu->lusup   = lusup;
 
282
    Glu->xlusup  = xlusup;
 
283
    Glu->ucol    = ucol;
 
284
    Glu->usub    = usub;
 
285
    Glu->xusub   = xusub;
 
286
    Glu->nzlmax  = nzlmax;
 
287
    Glu->nzumax  = nzumax;
 
288
    Glu->nzlumax = nzlumax;
 
289
    
 
290
    info = sLUWorkInit(m, n, panel_size, iwork, dwork, Glu->MemModel);
 
291
    if ( info )
 
292
        return ( info + smemory_usage(nzlmax, nzumax, nzlumax, n) + n);
 
293
    
 
294
    ++no_expand;
 
295
    return 0;
 
296
    
 
297
} /* sLUMemInit */
 
298
 
 
299
/* Allocate known working storage. Returns 0 if success, otherwise
 
300
   returns the number of bytes allocated so far when failure occurred. */
 
301
int
 
302
sLUWorkInit(int m, int n, int panel_size, int **iworkptr, 
 
303
            float **dworkptr, LU_space_t MemModel)
 
304
{
 
305
    int    isize, dsize, extra;
 
306
    float *old_ptr;
 
307
    int    maxsuper = sp_ienv(3),
 
308
           rowblk   = sp_ienv(4);
 
309
 
 
310
    isize = ( (2 * panel_size + 3 + NO_MARKER ) * m + n ) * sizeof(int);
 
311
    dsize = (m * panel_size +
 
312
             NUM_TEMPV(m,panel_size,maxsuper,rowblk)) * sizeof(float);
 
313
    
 
314
    if ( MemModel == SYSTEM ) 
 
315
        *iworkptr = (int *) intCalloc(isize/sizeof(int));
 
316
    else
 
317
        *iworkptr = (int *) suser_malloc(isize, TAIL);
 
318
    if ( ! *iworkptr ) {
 
319
        fprintf(stderr, "sLUWorkInit: malloc fails for local iworkptr[]\n");
 
320
        return (isize + n);
 
321
    }
 
322
 
 
323
    if ( MemModel == SYSTEM )
 
324
        *dworkptr = (float *) SUPERLU_MALLOC(dsize);
 
325
    else {
 
326
        *dworkptr = (float *) suser_malloc(dsize, TAIL);
 
327
        if ( NotDoubleAlign(*dworkptr) ) {
 
328
            old_ptr = *dworkptr;
 
329
            *dworkptr = (float*) DoubleAlign(*dworkptr);
 
330
            *dworkptr = (float*) ((double*)*dworkptr - 1);
 
331
            extra = (char*)old_ptr - (char*)*dworkptr;
 
332
#ifdef DEBUG        
 
333
            printf("sLUWorkInit: not aligned, extra %d\n", extra);
 
334
#endif      
 
335
            stack.top2 -= extra;
 
336
            stack.used += extra;
 
337
        }
 
338
    }
 
339
    if ( ! *dworkptr ) {
 
340
        fprintf(stderr, "malloc fails for local dworkptr[].");
 
341
        return (isize + dsize + n);
 
342
    }
 
343
        
 
344
    return 0;
 
345
}
 
346
 
 
347
 
 
348
/*
 
349
 * Set up pointers for real working arrays.
 
350
 */
 
351
void
 
352
sSetRWork(int m, int panel_size, float *dworkptr,
 
353
         float **dense, float **tempv)
 
354
{
 
355
    float zero = 0.0;
 
356
 
 
357
    int maxsuper = sp_ienv(3),
 
358
        rowblk   = sp_ienv(4);
 
359
    *dense = dworkptr;
 
360
    *tempv = *dense + panel_size*m;
 
361
    sfill (*dense, m * panel_size, zero);
 
362
    sfill (*tempv, NUM_TEMPV(m,panel_size,maxsuper,rowblk), zero);     
 
363
}
 
364
        
 
365
/*
 
366
 * Free the working storage used by factor routines.
 
367
 */
 
368
void sLUWorkFree(int *iwork, float *dwork, GlobalLU_t *Glu)
 
369
{
 
370
    if ( Glu->MemModel == SYSTEM ) {
 
371
        SUPERLU_FREE (iwork);
 
372
        SUPERLU_FREE (dwork);
 
373
    } else {
 
374
        stack.used -= (stack.size - stack.top2);
 
375
        stack.top2 = stack.size;
 
376
/*      sStackCompress(Glu);  */
 
377
    }
 
378
    
 
379
    SUPERLU_FREE (expanders);   
 
380
    expanders = 0;
 
381
}
 
382
 
 
383
/* Expand the data structures for L and U during the factorization.
 
384
 * Return value:   0 - successful return
 
385
 *               > 0 - number of bytes allocated when run out of space
 
386
 */
 
387
int
 
388
sLUMemXpand(int jcol,
 
389
           int next,          /* number of elements currently in the factors */
 
390
           MemType mem_type,  /* which type of memory to expand  */
 
391
           int *maxlen,       /* modified - maximum length of a data structure */
 
392
           GlobalLU_t *Glu    /* modified - global LU data structures */
 
393
           )
 
394
{
 
395
    void   *new_mem;
 
396
    
 
397
#ifdef DEBUG    
 
398
    printf("sLUMemXpand(): jcol %d, next %d, maxlen %d, MemType %d\n",
 
399
           jcol, next, *maxlen, mem_type);
 
400
#endif    
 
401
 
 
402
    if (mem_type == USUB) 
 
403
        new_mem = sexpand(maxlen, mem_type, next, 1, Glu);
 
404
    else
 
405
        new_mem = sexpand(maxlen, mem_type, next, 0, Glu);
 
406
    
 
407
    if ( !new_mem ) {
 
408
        int    nzlmax  = Glu->nzlmax;
 
409
        int    nzumax  = Glu->nzumax;
 
410
        int    nzlumax = Glu->nzlumax;
 
411
        fprintf(stderr, "Can't expand MemType %d: jcol %d\n", mem_type, jcol);
 
412
        return (smemory_usage(nzlmax, nzumax, nzlumax, Glu->n) + Glu->n);
 
413
    }
 
414
 
 
415
    switch ( mem_type ) {
 
416
      case LUSUP:
 
417
        Glu->lusup   = (float *) new_mem;
 
418
        Glu->nzlumax = *maxlen;
 
419
        break;
 
420
      case UCOL:
 
421
        Glu->ucol   = (float *) new_mem;
 
422
        Glu->nzumax = *maxlen;
 
423
        break;
 
424
      case LSUB:
 
425
        Glu->lsub   = (int *) new_mem;
 
426
        Glu->nzlmax = *maxlen;
 
427
        break;
 
428
      case USUB:
 
429
        Glu->usub   = (int *) new_mem;
 
430
        Glu->nzumax = *maxlen;
 
431
        break;
 
432
    }
 
433
    
 
434
    return 0;
 
435
    
 
436
}
 
437
 
 
438
 
 
439
 
 
440
void
 
441
copy_mem_float(int howmany, void *old, void *new)
 
442
{
 
443
    register int i;
 
444
    float *dold = old;
 
445
    float *dnew = new;
 
446
    for (i = 0; i < howmany; i++) dnew[i] = dold[i];
 
447
}
 
448
 
 
449
/*
 
450
 * Expand the existing storage to accommodate more fill-ins.
 
451
 */
 
452
void
 
453
*sexpand (
 
454
         int *prev_len,   /* length used from previous call */
 
455
         MemType type,    /* which part of the memory to expand */
 
456
         int len_to_copy, /* size of the memory to be copied to new store */
 
457
         int keep_prev,   /* = 1: use prev_len;
 
458
                             = 0: compute new_len to expand */
 
459
         GlobalLU_t *Glu  /* modified - global LU data structures */
 
460
        )
 
461
{
 
462
    float    EXPAND = 1.5;
 
463
    float    alpha;
 
464
    void     *new_mem, *old_mem;
 
465
    int      new_len, tries, lword, extra, bytes_to_copy;
 
466
 
 
467
    alpha = EXPAND;
 
468
 
 
469
    if ( no_expand == 0 || keep_prev ) /* First time allocate requested */
 
470
        new_len = *prev_len;
 
471
    else {
 
472
        new_len = alpha * *prev_len;
 
473
    }
 
474
    
 
475
    if ( type == LSUB || type == USUB ) lword = sizeof(int);
 
476
    else lword = sizeof(float);
 
477
 
 
478
    if ( Glu->MemModel == SYSTEM ) {
 
479
        new_mem = (void *) SUPERLU_MALLOC(new_len * lword);
 
480
/*      new_mem = (void *) calloc(new_len, lword); */
 
481
        if ( no_expand != 0 ) {
 
482
            tries = 0;
 
483
            if ( keep_prev ) {
 
484
                if ( !new_mem ) return (NULL);
 
485
            } else {
 
486
                while ( !new_mem ) {
 
487
                    if ( ++tries > 10 ) return (NULL);
 
488
                    alpha = Reduce(alpha);
 
489
                    new_len = alpha * *prev_len;
 
490
                    new_mem = (void *) SUPERLU_MALLOC(new_len * lword); 
 
491
/*                  new_mem = (void *) calloc(new_len, lword); */
 
492
                }
 
493
            }
 
494
            if ( type == LSUB || type == USUB ) {
 
495
                copy_mem_int(len_to_copy, expanders[type].mem, new_mem);
 
496
            } else {
 
497
                copy_mem_float(len_to_copy, expanders[type].mem, new_mem);
 
498
            }
 
499
            SUPERLU_FREE (expanders[type].mem);
 
500
        }
 
501
        expanders[type].mem = (void *) new_mem;
 
502
        
 
503
    } else { /* MemModel == USER */
 
504
        if ( no_expand == 0 ) {
 
505
            new_mem = suser_malloc(new_len * lword, HEAD);
 
506
            if ( NotDoubleAlign(new_mem) &&
 
507
                (type == LUSUP || type == UCOL) ) {
 
508
                old_mem = new_mem;
 
509
                new_mem = (void *)DoubleAlign(new_mem);
 
510
                extra = (char*)new_mem - (char*)old_mem;
 
511
#ifdef DEBUG            
 
512
                printf("expand(): not aligned, extra %d\n", extra);
 
513
#endif          
 
514
                stack.top1 += extra;
 
515
                stack.used += extra;
 
516
            }
 
517
            expanders[type].mem = (void *) new_mem;
 
518
        }
 
519
        else {
 
520
            tries = 0;
 
521
            extra = (new_len - *prev_len) * lword;
 
522
            if ( keep_prev ) {
 
523
                if ( StackFull(extra) ) return (NULL);
 
524
            } else {
 
525
                while ( StackFull(extra) ) {
 
526
                    if ( ++tries > 10 ) return (NULL);
 
527
                    alpha = Reduce(alpha);
 
528
                    new_len = alpha * *prev_len;
 
529
                    extra = (new_len - *prev_len) * lword;          
 
530
                }
 
531
            }
 
532
 
 
533
            if ( type != USUB ) {
 
534
                new_mem = (void*)((char*)expanders[type + 1].mem + extra);
 
535
                bytes_to_copy = (char*)stack.array + stack.top1
 
536
                    - (char*)expanders[type + 1].mem;
 
537
                user_bcopy(expanders[type+1].mem, new_mem, bytes_to_copy);
 
538
 
 
539
                if ( type < USUB ) {
 
540
                    Glu->usub = expanders[USUB].mem =
 
541
                        (void*)((char*)expanders[USUB].mem + extra);
 
542
                }
 
543
                if ( type < LSUB ) {
 
544
                    Glu->lsub = expanders[LSUB].mem =
 
545
                        (void*)((char*)expanders[LSUB].mem + extra);
 
546
                }
 
547
                if ( type < UCOL ) {
 
548
                    Glu->ucol = expanders[UCOL].mem =
 
549
                        (void*)((char*)expanders[UCOL].mem + extra);
 
550
                }
 
551
                stack.top1 += extra;
 
552
                stack.used += extra;
 
553
                if ( type == UCOL ) {
 
554
                    stack.top1 += extra;   /* Add same amount for USUB */
 
555
                    stack.used += extra;
 
556
                }
 
557
                
 
558
            } /* if ... */
 
559
 
 
560
        } /* else ... */
 
561
    }
 
562
 
 
563
    expanders[type].size = new_len;
 
564
    *prev_len = new_len;
 
565
    if ( no_expand ) ++no_expand;
 
566
    
 
567
    return (void *) expanders[type].mem;
 
568
    
 
569
} /* sexpand */
 
570
 
 
571
 
 
572
/*
 
573
 * Compress the work[] array to remove fragmentation.
 
574
 */
 
575
void
 
576
sStackCompress(GlobalLU_t *Glu)
 
577
{
 
578
    register int iword, dword, ndim;
 
579
    char    *last, *fragment;
 
580
    int      *ifrom, *ito;
 
581
    float   *dfrom, *dto;
 
582
    int      *xlsub, *lsub, *xusub, *usub, *xlusup;
 
583
    float   *ucol, *lusup;
 
584
    
 
585
    iword = sizeof(int);
 
586
    dword = sizeof(float);
 
587
    ndim = Glu->n;
 
588
 
 
589
    xlsub  = Glu->xlsub;
 
590
    lsub   = Glu->lsub;
 
591
    xusub  = Glu->xusub;
 
592
    usub   = Glu->usub;
 
593
    xlusup = Glu->xlusup;
 
594
    ucol   = Glu->ucol;
 
595
    lusup  = Glu->lusup;
 
596
    
 
597
    dfrom = ucol;
 
598
    dto = (float *)((char*)lusup + xlusup[ndim] * dword);
 
599
    copy_mem_float(xusub[ndim], dfrom, dto);
 
600
    ucol = dto;
 
601
 
 
602
    ifrom = lsub;
 
603
    ito = (int *) ((char*)ucol + xusub[ndim] * iword);
 
604
    copy_mem_int(xlsub[ndim], ifrom, ito);
 
605
    lsub = ito;
 
606
    
 
607
    ifrom = usub;
 
608
    ito = (int *) ((char*)lsub + xlsub[ndim] * iword);
 
609
    copy_mem_int(xusub[ndim], ifrom, ito);
 
610
    usub = ito;
 
611
    
 
612
    last = (char*)usub + xusub[ndim] * iword;
 
613
    fragment = (char*) (((char*)stack.array + stack.top1) - last);
 
614
    stack.used -= (long int) fragment;
 
615
    stack.top1 -= (long int) fragment;
 
616
 
 
617
    Glu->ucol = ucol;
 
618
    Glu->lsub = lsub;
 
619
    Glu->usub = usub;
 
620
    
 
621
#ifdef DEBUG
 
622
    printf("sStackCompress: fragment %d\n", fragment);
 
623
    /* for (last = 0; last < ndim; ++last)
 
624
        print_lu_col("After compress:", last, 0);*/
 
625
#endif    
 
626
    
 
627
}
 
628
 
 
629
/*
 
630
 * Allocate storage for original matrix A
 
631
 */
 
632
void
 
633
sallocateA(int n, int nnz, float **a, int **asub, int **xa)
 
634
{
 
635
    *a    = (float *) floatMalloc(nnz);
 
636
    *asub = (int *) intMalloc(nnz);
 
637
    *xa   = (int *) intMalloc(n+1);
 
638
}
 
639
 
 
640
 
 
641
float *floatMalloc(int n)
 
642
{
 
643
    float *buf;
 
644
    buf = (float *) SUPERLU_MALLOC(n * sizeof(float)); 
 
645
    if ( !buf ) {
 
646
        ABORT("SUPERLU_MALLOC failed for buf in floatMalloc()\n");
 
647
    }
 
648
    return (buf);
 
649
}
 
650
 
 
651
float *floatCalloc(int n)
 
652
{
 
653
    float *buf;
 
654
    register int i;
 
655
    float zero = 0.0;
 
656
    buf = (float *) SUPERLU_MALLOC(n * sizeof(float));
 
657
    if ( !buf ) {
 
658
        ABORT("SUPERLU_MALLOC failed for buf in floatCalloc()\n");
 
659
    }
 
660
    for (i = 0; i < n; ++i) buf[i] = zero;
 
661
    return (buf);
 
662
}
 
663
 
 
664
 
 
665
int smemory_usage(const int nzlmax, const int nzumax, 
 
666
                  const int nzlumax, const int n)
 
667
{
 
668
    register int iword, dword;
 
669
 
 
670
    iword   = sizeof(int);
 
671
    dword   = sizeof(float);
 
672
    
 
673
    return (10 * n * iword +
 
674
            nzlmax * iword + nzumax * (iword + dword) + nzlumax * dword);
 
675
 
 
676
}