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