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

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/umfpack/umf_garbage_collection.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
* 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:
1
 
/* ========================================================================== */
2
 
/* === UMF_garbage_collection =============================================== */
3
 
/* ========================================================================== */
4
 
 
5
 
/* -------------------------------------------------------------------------- */
6
 
/* UMFPACK Version 4.1 (Apr. 30, 2003), Copyright (c) 2003 by Timothy A.      */
7
 
/* Davis.  All Rights Reserved.  See ../README for License.                   */
8
 
/* email: davis@cise.ufl.edu    CISE Department, Univ. of Florida.            */
9
 
/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
10
 
/* -------------------------------------------------------------------------- */
11
 
 
12
 
/*
13
 
    Compress the elements at the tail of Numeric->Memory, and delete the tuples.
14
 
    Elements are renumbered.  The new numbering space is compressed, and
15
 
    in the order of element creation (original elements of A first, followed
16
 
    by the new elements in the order that they were formed).
17
 
 
18
 
    Only called by UMF_get_memory.
19
 
 
20
 
    There are 5 ways in which garbage collection can be performed:
21
 
 
22
 
        Allocate a new working array for the current frontal matrix.  In this
23
 
        case, there are never any pivot rows/columns in the current frontal
24
 
        matrix (fnpiv = 0), and the old working array for the current frontal
25
 
        matrix can always be fully compacted, to fnrows-by-fncols.
26
 
 
27
 
            UMF_kernel : UMF_extend : UMF_grow_front : UMF_get_memory
28
 
            UMF_kernel : UMF_init_front : UMF_grow_front : UMF_get_memory
29
 
            UMF_kernel : UMF_start_front : UMF_grow_front : UMF_get_memory
30
 
 
31
 
        Allocate a new element.  In this case, UMF_grow_front may or may not
32
 
        be subsequently called, depending on Work->do_grow.  There are never
33
 
        any pivot rows/columns in the current frontal matrix (fnpiv=0), but one
34
 
        may be added if UMF_init_front is to be called just after
35
 
        UMF_create_element.  If do_grow is true, then the current front can be
36
 
        fully compacted, to fnrows-by-fncols.  Otherwise, it can only be
37
 
        partially compacted, to MAX (fnrows, fnrows_new + 1) -by-
38
 
        MAX (fncols, fncols_new + 1).
39
 
 
40
 
            UMF_kernel : UMF_create_element : UMF_get_memory
41
 
 
42
 
        Allocate rows of L and columns of U.  In this case, the current
43
 
        frontal matrix is only partially compacted, to (fnrows_new + 1)-by-
44
 
        (fncols_new + 1).  There are pivots in the frontal matrix (fnpiv > 0).
45
 
 
46
 
            UMF_kernel : UMF_store_lu : UMF_get_memory
47
 
*/
48
 
 
49
 
#include "umf_internal.h"
50
 
 
51
 
GLOBAL void UMF_garbage_collection
52
 
(
53
 
    NumericType *Numeric,
54
 
    WorkType *Work,
55
 
    Int drnew,      /* compact current front to drnew-by-dcnew */
56
 
    Int dcnew,
57
 
    Int do_Fcpos
58
 
)
59
 
{
60
 
    /* ---------------------------------------------------------------------- */
61
 
    /* local variables */
62
 
    /* ---------------------------------------------------------------------- */
63
 
 
64
 
    Int size, e, n_row, n_col, nrows, ncols, nrowsleft, ncolsleft, prevsize,
65
 
        csize, size2, i2, j2, i, j, cdeg, rdeg, *E, row, col,
66
 
        *Rows, *Cols, *Rows2, *Cols2, nel, e2, *Row_tuples, *Col_tuples,
67
 
        *Row_degree, *Col_degree ;
68
 
    Entry *C, *C1, *C3, *C2 ;
69
 
    Unit *psrc, *pdest, *p, *pnext ;
70
 
    Element *epsrc, *epdest ;
71
 
 
72
 
#ifndef NDEBUG
73
 
    Int nmark ;
74
 
#endif
75
 
 
76
 
    /* ---------------------------------------------------------------------- */
77
 
    /* get parameters */
78
 
    /* ---------------------------------------------------------------------- */
79
 
 
80
 
    Col_degree = Numeric->Cperm ;       /* for NON_PIVOTAL_COL macro */
81
 
    Row_degree = Numeric->Rperm ;       /* for NON_PIVOTAL_ROW macro */
82
 
    Row_tuples = Numeric->Uip ;
83
 
    Col_tuples = Numeric->Lip ;
84
 
    E = Work->E ;
85
 
    n_row = Work->n_row ;
86
 
    n_col = Work->n_col ;
87
 
 
88
 
    /* note that the tuple lengths (Col_tlen and Row_tlen) are updated, but */
89
 
    /* the tuple lists themselves are stale and are about to be destroyed */
90
 
    /* and recreated.  Do not attempt to scan them until they are recreated. */
91
 
 
92
 
#ifndef NDEBUG
93
 
    DEBUGm1 (("::::GARBAGE COLLECTION::::\n")) ;
94
 
    UMF_dump_memory (Numeric) ;
95
 
#endif
96
 
 
97
 
    Numeric->ngarbage++ ;
98
 
 
99
 
    /* ---------------------------------------------------------------------- */
100
 
    /* delete the tuple lists by marking the blocks as free */
101
 
    /* ---------------------------------------------------------------------- */
102
 
 
103
 
    /* do not modify Row_tlen and Col_tlen */
104
 
    /* those are needed for UMF_build_tuples */
105
 
 
106
 
    for (row = 0 ; row < n_row ; row++)
107
 
    {
108
 
        if (NON_PIVOTAL_ROW (row) && Row_tuples [row])
109
 
        {
110
 
            DEBUG2 (("row "ID" tuples "ID"\n", row, Row_tuples [row])) ;
111
 
            p = Numeric->Memory + Row_tuples [row] - 1 ;
112
 
            DEBUG2 (("Freeing tuple list row "ID", p-S "ID", size "ID"\n",
113
 
                row, (Int) (p-Numeric->Memory), p->header.size)) ;
114
 
            ASSERT (p->header.size > 0) ;
115
 
            ASSERT (p >= Numeric->Memory + Numeric->itail) ;
116
 
            ASSERT (p < Numeric->Memory + Numeric->size) ;
117
 
            p->header.size = -p->header.size ;
118
 
            Row_tuples [row] = 0 ;
119
 
        }
120
 
    }
121
 
 
122
 
    for (col = 0 ; col < n_col ; col++)
123
 
    {
124
 
        if (NON_PIVOTAL_COL (col) && Col_tuples [col])
125
 
        {
126
 
            DEBUG2 (("col "ID" tuples "ID"\n", col, Col_tuples [col])) ;
127
 
            p = Numeric->Memory + Col_tuples [col] - 1 ;
128
 
            DEBUG2 (("Freeing tuple list col "ID", p-S "ID", size "ID"\n",
129
 
                col, (Int) (p-Numeric->Memory), p->header.size)) ;
130
 
            ASSERT (p->header.size > 0) ;
131
 
            ASSERT (p >= Numeric->Memory + Numeric->itail) ;
132
 
            ASSERT (p < Numeric->Memory + Numeric->size) ;
133
 
            p->header.size = -p->header.size ;
134
 
            Col_tuples [col] = 0 ;
135
 
        }
136
 
    }
137
 
 
138
 
    /* ---------------------------------------------------------------------- */
139
 
    /* mark the elements, and compress the name space */
140
 
    /* ---------------------------------------------------------------------- */
141
 
 
142
 
    nel = Work->nel ;
143
 
    ASSERT (nel < Work->elen) ;
144
 
 
145
 
#ifndef NDEBUG
146
 
    nmark = 0 ;
147
 
    UMF_dump_current_front (Numeric, Work, FALSE) ;
148
 
    DEBUGm1 (("E [0] "ID"  \n", E [0])) ;
149
 
    ASSERT (IMPLIES (E [0],
150
 
                Work->Flublock == (Entry *) (Numeric->Memory + E [0]))) ;
151
 
    ASSERT (IMPLIES (Work->Flublock,
152
 
                Work->Flublock == (Entry *) (Numeric->Memory + E [0]))) ;
153
 
    ASSERT ((E [0] != 0) == (Work->Flublock != (Entry *) NULL)) ;
154
 
#endif
155
 
 
156
 
    e2 = 0 ;
157
 
 
158
 
    for (e = 0 ; e <= nel ; e++) /* for all elements in order of creation */
159
 
    {
160
 
        if (E [e])
161
 
        {
162
 
            psrc = Numeric->Memory + E [e] ;
163
 
            psrc-- ;            /* get the header of this block */
164
 
            if (e > 0)
165
 
            {
166
 
                e2++ ;  /* do not renumber element zero */
167
 
            }
168
 
            ASSERT (psrc->header.size > 0) ;
169
 
            psrc->header.size = e2  ;   /* store the new name in the header */
170
 
#ifndef NDEBUG
171
 
            nmark++ ;
172
 
#endif
173
 
            DEBUG7 ((ID":: Mark e "ID" at psrc-S "ID", new e "ID"\n",
174
 
                nmark, e, (Int) (psrc-Numeric->Memory), e2)) ;
175
 
            E [e] = 0 ;
176
 
            if (e == Work->prior_element)
177
 
            {
178
 
                Work->prior_element = e2 ;
179
 
            }
180
 
        }
181
 
    }
182
 
 
183
 
    /* all 1..e2 are now in use (element zero may or may not be in use) */
184
 
    Work->nel = e2 ;
185
 
    nel = Work->nel ;
186
 
 
187
 
#ifndef NDEBUG
188
 
    for (e = 0 ; e < Work->elen ; e++)
189
 
    {
190
 
        ASSERT (!E [e]) ;
191
 
    }
192
 
#endif
193
 
 
194
 
    /* ---------------------------------------------------------------------- */
195
 
    /* compress the elements */
196
 
    /* ---------------------------------------------------------------------- */
197
 
 
198
 
    /* point to tail marker block of size 1 + header */
199
 
    psrc = Numeric->Memory + Numeric->size - 2 ;
200
 
    pdest = psrc ;
201
 
    prevsize = psrc->header.prevsize ;
202
 
    DEBUG7 (("Starting the compression:\n")) ;
203
 
 
204
 
    while (prevsize > 0)
205
 
    {
206
 
 
207
 
        /* ------------------------------------------------------------------ */
208
 
        /* move up to the next element above the current header, and */
209
 
        /* get the element name and size */
210
 
        /* (if it is an element, the name will be positive) */
211
 
        /* ------------------------------------------------------------------ */
212
 
 
213
 
        size = prevsize ;
214
 
        psrc -= (size + 1) ;
215
 
        e = psrc->header.size ;
216
 
        prevsize = psrc->header.prevsize ;
217
 
        /* top block at tail has prevsize of 0 */
218
 
 
219
 
        /* a free block will have a negative size, so skip it */
220
 
        /* otherwise, if size >= 0, it holds the element name, not the size */
221
 
 
222
 
        DEBUG8 (("psrc-S: "ID" prevsize: "ID" size: "ID,
223
 
            (Int) (psrc-Numeric->Memory), prevsize, size)) ;
224
 
 
225
 
        if (e == 0)
226
 
        {
227
 
            /* -------------------------------------------------------------- */
228
 
            /* this is the current frontal matrix */
229
 
            /* -------------------------------------------------------------- */
230
 
 
231
 
            Entry *F1, *F2, *Fsrc, *Fdst ;
232
 
            Int c, r, k, dr, dc, gap, gap1, gap2, nb ;
233
 
 
234
 
            /* shift the frontal matrix down */
235
 
            F1 = (Entry *) (psrc + 1) ;
236
 
 
237
 
            /* get the size of the current front.  r and c could be zero */
238
 
            k = Work->fnpiv ;
239
 
            dr = Work->fnr_curr ;
240
 
            dc = Work->fnc_curr ;
241
 
            r = Work->fnrows ;
242
 
            c = Work->fncols ;
243
 
            nb = Work->nb ;
244
 
 
245
 
            ASSERT ((dr >= 0 && (dr % 2) == 1) || dr == 0) ;
246
 
            ASSERT (drnew >= 0) ;
247
 
            if (drnew % 2 == 0)
248
 
            {
249
 
                /* make sure leading frontal matrix dimension is always odd */
250
 
                drnew++ ;
251
 
            }
252
 
            drnew = MIN (dr, drnew) ;
253
 
            ASSERT ((drnew >= 0 && (drnew % 2) == 1) || drnew == 0) ;
254
 
 
255
 
            pnext = pdest ;
256
 
 
257
 
#ifndef NDEBUG
258
 
            DEBUGm2 (("move front: dr "ID" dc "ID" r "ID" drnew "ID" c "ID
259
 
                " dcnew " ID" k "ID"\n", dr, dc, r, drnew, c, dcnew, k)) ;
260
 
            DEBUG7 (("\n")) ;
261
 
            DEBUG7 ((ID":: Move current frontal matrix from: psrc-S: "ID" \n",
262
 
                nmark, (Int) (psrc-Numeric->Memory))) ;
263
 
            nmark-- ;
264
 
            ASSERT (E [e] == 0) ;
265
 
            ASSERT (Work->Flublock == F1) ;
266
 
            ASSERT (Work->Flblock  == Work->Flublock + nb*nb) ;
267
 
            ASSERT (Work->Fublock  == Work->Flblock  + dr*nb) ;
268
 
            ASSERT (Work->Fcblock  == Work->Fublock  + nb*dc) ;
269
 
            DEBUG7 (("C  block: ")) ;
270
 
            UMF_dump_dense (Work->Fcblock,  dr, r, c) ;
271
 
            DEBUG7 (("L  block: ")) ;
272
 
            UMF_dump_dense (Work->Flblock,  dr, r, k);
273
 
            DEBUG7 (("U' block: ")) ;
274
 
            UMF_dump_dense (Work->Fublock,  dc, c, k) ;
275
 
            DEBUG7 (("LU block: ")) ;
276
 
            UMF_dump_dense (Work->Flublock, nb, k, k) ;
277
 
            ASSERT (r <= drnew && c <= dcnew && drnew <= dr && dcnew <= dc) ;
278
 
#endif
279
 
 
280
 
            /* compact frontal matrix to drnew-by-dcnew before moving it */
281
 
 
282
 
            /* do not compact the LU block (nb-by-nb) */
283
 
 
284
 
            /* compact the columns of L (from dr-by-nb to drnew-by-nb) */
285
 
            Fsrc = Work->Flblock ;
286
 
            Fdst = Work->Flblock ;
287
 
            ASSERT (Fdst == F1 + nb*nb) ;
288
 
            gap1 = dr - r ;
289
 
            gap2 = drnew - r ;
290
 
            ASSERT (gap1 >= 0) ;
291
 
            for (j = 0 ; j < k ; j++)
292
 
            {
293
 
                for (i = 0 ; i < r ; i++)
294
 
                {
295
 
                    *Fdst++ = *Fsrc++ ;
296
 
                }
297
 
                Fsrc += gap1 ;
298
 
                Fdst += gap2 ;
299
 
            }
300
 
            ASSERT (Fdst == F1 + nb*nb + drnew*k) ;
301
 
            Fdst += drnew * (nb - k) ;
302
 
 
303
 
            /* compact the rows of U (U' from dc-by-nb to dcnew-by-nb) */
304
 
            Fsrc = Work->Fublock ;
305
 
            ASSERT (Fdst == F1 + nb*nb + drnew*nb) ;
306
 
            gap1 = dc - c ;
307
 
            gap2 = dcnew - c ;
308
 
            for (i = 0 ; i < k ; i++)
309
 
            {
310
 
                for (j = 0 ; j < c ; j++)
311
 
                {
312
 
                    *Fdst++ = *Fsrc++ ;
313
 
                }
314
 
                Fsrc += gap1 ;
315
 
                Fdst += gap2 ;
316
 
            }
317
 
            ASSERT (Fdst == F1 + nb*nb + drnew*nb + dcnew*k) ;
318
 
            Fdst += dcnew * (nb - k) ;
319
 
 
320
 
            /* compact the columns of C (from dr-by-dc to drnew-by-dcnew) */
321
 
            Fsrc = Work->Fcblock ;
322
 
            ASSERT (Fdst == F1 + nb*nb + drnew*nb + nb*dcnew) ;
323
 
            gap1 = dr - r ;
324
 
            gap2 = drnew - r ;
325
 
            for (j = 0 ; j < c ; j++)
326
 
            {
327
 
                for (i = 0 ; i < r ; i++)
328
 
                {
329
 
                    *Fdst++ = *Fsrc++ ;
330
 
                }
331
 
                Fsrc += gap1 ;
332
 
                Fdst += gap2 ;
333
 
            }
334
 
            ASSERT (Fdst == F1 + nb*nb + drnew*nb + nb*dcnew + drnew*c) ;
335
 
 
336
 
            /* recompute Fcpos, if necessary */
337
 
            if (do_Fcpos)
338
 
            {
339
 
                Int *Fcols, *Fcpos ;
340
 
                Fcols = Work->Fcols ;
341
 
                Fcpos = Work->Fcpos ;
342
 
                for (j = 0 ; j < c ; j++)
343
 
                {
344
 
                    col = Fcols [j] ;
345
 
                    ASSERT (col >= 0 && col < Work->n_col) ;
346
 
                    ASSERT (Fcpos [col] == j * dr) ;
347
 
                    Fcpos [col] = j * drnew ;
348
 
                }
349
 
#ifndef NDEBUG
350
 
                {
351
 
                    Int cnt = 0 ;
352
 
                    for (j = 0 ; j < Work->n_col ; j++)
353
 
                    {
354
 
                        if (Fcpos [j] != EMPTY) cnt++ ;
355
 
                    }
356
 
                    DEBUGm2 (("Recompute Fcpos cnt "ID" c "ID"\n", cnt, c)) ;
357
 
                    ASSERT (cnt == c) ;
358
 
                }
359
 
#endif
360
 
            }
361
 
 
362
 
#ifndef NDEBUG
363
 
            DEBUGm2 (("Compacted front, drnew "ID" dcnew "ID"\n", drnew, dcnew)) ;
364
 
            DEBUG7 (("C  block: ")) ;
365
 
            UMF_dump_dense (F1 + nb*nb + drnew*nb + nb*dcnew, drnew, r, c) ;
366
 
            DEBUG7 (("L  block: ")) ;
367
 
            UMF_dump_dense (F1 + nb*nb, drnew, r, k) ;
368
 
            DEBUG7 (("U  block: ")) ;
369
 
            UMF_dump_dense (F1 + nb*nb + drnew*nb, nb, k, c) ;
370
 
            DEBUG7 (("LU block: ")) ;
371
 
            UMF_dump_dense (F1, nb, k, k) ;
372
 
#endif
373
 
 
374
 
            /* Compacted dimensions of the new frontal matrix. */
375
 
            Work->fnr_curr = drnew ;
376
 
            Work->fnc_curr = dcnew ;
377
 
            Work->fcurr_size = (drnew + nb) * (dcnew + nb) ;
378
 
            size = UNITS (Entry, Work->fcurr_size) ;
379
 
 
380
 
            /* make sure the object doesn't evaporate.  The front can have
381
 
             * zero size (Work->fcurr_size = 0), but the size of the memory
382
 
             * block containing it cannot have zero size. */
383
 
            size = MAX (1, size) ;
384
 
 
385
 
            /* get the destination of frontal matrix */
386
 
            pnext->header.prevsize = size ;
387
 
            pdest -= (size + 1) ;
388
 
            F2 = (Entry *) (pdest + 1) ;
389
 
 
390
 
            ASSERT ((unsigned Int) psrc + 1 + size <= (unsigned Int) pnext) ;
391
 
            ASSERT (psrc <= pdest) ;
392
 
            ASSERT (F1 <= F2) ;
393
 
 
394
 
            /* move the C block first */
395
 
            Fsrc = F1 + nb*nb + drnew*nb + nb*dcnew + drnew*c ;
396
 
            Fdst = F2 + nb*nb + drnew*nb + nb*dcnew + drnew*c ;
397
 
            gap = drnew - r ;
398
 
            for (j = c-1 ; j >= 0 ; j--)
399
 
            {
400
 
                Fsrc -= gap ;
401
 
                Fdst -= gap ;
402
 
                /* move column j of C */
403
 
                for (i = r-1 ; i >= 0 ; i--)
404
 
                {
405
 
                    *--Fdst = *--Fsrc ;
406
 
                }
407
 
            }
408
 
            ASSERT (Fsrc == F1 + nb*nb + drnew*nb + nb*dcnew) ;
409
 
            ASSERT (Fdst == F2 + nb*nb + drnew*nb + nb*dcnew) ;
410
 
 
411
 
            /* move the U block */
412
 
            Fsrc -= dcnew * (nb - k) ;
413
 
            Fdst -= dcnew * (nb - k) ;
414
 
            ASSERT (Fsrc == F1 + nb*nb + drnew*nb + dcnew*k) ;
415
 
            ASSERT (Fdst == F2 + nb*nb + drnew*nb + dcnew*k) ;
416
 
            gap = dcnew - c ;
417
 
            for (i = k-1 ; i >= 0 ; i--)
418
 
            {
419
 
                Fsrc -= gap ;
420
 
                Fdst -= gap ;
421
 
                for (j = c-1 ; j >= 0 ; j--)
422
 
                {
423
 
                    *--Fdst = *--Fsrc ;
424
 
                }
425
 
            }
426
 
            ASSERT (Fsrc == F1 + nb*nb + drnew*nb) ;
427
 
            ASSERT (Fdst == F2 + nb*nb + drnew*nb) ;
428
 
 
429
 
            /* move the L block */
430
 
            Fsrc -= drnew * (nb - k) ;
431
 
            Fdst -= drnew * (nb - k) ;
432
 
            ASSERT (Fsrc == F1 + nb*nb + drnew*k) ;
433
 
            ASSERT (Fdst == F2 + nb*nb + drnew*k) ;
434
 
            gap = drnew - r ;
435
 
            for (j = k-1 ; j >= 0 ; j--)
436
 
            {
437
 
                Fsrc -= gap ;
438
 
                Fdst -= gap ;
439
 
                for (i = r-1 ; i >= 0 ; i--)
440
 
                {
441
 
                    *--Fdst = *--Fsrc ;
442
 
                }
443
 
            }
444
 
            ASSERT (Fsrc == F1 + nb*nb) ;
445
 
            ASSERT (Fdst == F2 + nb*nb) ;
446
 
 
447
 
            /* move the LU block */
448
 
            Fsrc -= nb * (nb - k) ;
449
 
            Fdst -= nb * (nb - k) ;
450
 
            ASSERT (Fsrc == F1 + nb*k) ;
451
 
            ASSERT (Fdst == F2 + nb*k) ;
452
 
            gap = nb - k ;
453
 
            for (j = k-1 ; j >= 0 ; j--)
454
 
            {
455
 
                Fsrc -= gap ;
456
 
                Fdst -= gap ;
457
 
                for (i = k-1 ; i >= 0 ; i--)
458
 
                {
459
 
                    *--Fdst = *--Fsrc ;
460
 
                }
461
 
            }
462
 
            ASSERT (Fsrc == F1) ;
463
 
            ASSERT (Fdst == F2) ;
464
 
 
465
 
            E [0] = (pdest + 1) - Numeric->Memory ;
466
 
 
467
 
            Work->Flublock = (Entry *) (Numeric->Memory + E [0]) ;
468
 
            ASSERT (Work->Flublock == F2) ;
469
 
            Work->Flblock  = Work->Flublock + nb * nb ;
470
 
            Work->Fublock  = Work->Flblock  + drnew * nb ;
471
 
            Work->Fcblock  = Work->Fublock  + nb * dcnew ;
472
 
 
473
 
            pdest->header.prevsize = 0 ;
474
 
            pdest->header.size = size ;
475
 
 
476
 
#ifndef NDEBUG
477
 
            DEBUG7 (("After moving compressed current frontal matrix:\n")) ;
478
 
            DEBUG7 (("C  block: ")) ;
479
 
            UMF_dump_dense (Work->Fcblock,  drnew, r, c) ;
480
 
            DEBUG7 (("L  block: ")) ;
481
 
            UMF_dump_dense (Work->Flblock,  drnew, r, k);
482
 
            DEBUG7 (("U' block: ")) ;
483
 
            UMF_dump_dense (Work->Fublock,  dcnew, c, k) ;
484
 
            DEBUG7 (("LU block: ")) ;
485
 
            UMF_dump_dense (Work->Flublock, nb, k, k) ;
486
 
#endif
487
 
 
488
 
        }
489
 
        else if (e > 0)
490
 
        {
491
 
 
492
 
            /* -------------------------------------------------------------- */
493
 
            /* this is an element, compress and move from psrc down to pdest */
494
 
            /* -------------------------------------------------------------- */
495
 
 
496
 
#ifndef NDEBUG
497
 
            DEBUG7 (("\n")) ;
498
 
            DEBUG7 ((ID":: Move element "ID": from: "ID" \n",
499
 
                nmark, e, (Int) (psrc-Numeric->Memory))) ;
500
 
            nmark-- ;
501
 
            ASSERT (e <= nel) ;
502
 
            ASSERT (E [e] == 0) ;
503
 
#endif
504
 
 
505
 
            /* -------------------------------------------------------------- */
506
 
            /* get the element scalars, and pointers to C, Rows, and Cols: */
507
 
            /* -------------------------------------------------------------- */
508
 
 
509
 
            p = psrc + 1 ;
510
 
            GET_ELEMENT (epsrc, p, Cols, Rows, ncols, nrows, C) ;
511
 
            nrowsleft = epsrc->nrowsleft ;
512
 
            ncolsleft = epsrc->ncolsleft ;
513
 
            cdeg = epsrc->cdeg ;
514
 
            rdeg = epsrc->rdeg ;
515
 
 
516
 
#ifndef NDEBUG
517
 
            DEBUG7 ((" nrows "ID" nrowsleft "ID"\n", nrows, nrowsleft)) ;
518
 
            DEBUG7 ((" ncols "ID" ncolsleft "ID"\n", ncols, ncolsleft)) ;
519
 
            DEBUG8 ((" Rows:")) ;
520
 
            for (i = 0 ; i < nrows ; i++) DEBUG8 ((" "ID, Rows [i])) ;
521
 
            DEBUG8 (("\n Cols:")) ;
522
 
            for (j = 0 ; j < ncols ; j++) DEBUG8 ((" "ID, Cols [j])) ;
523
 
            DEBUG8 (("\n")) ;
524
 
#endif
525
 
 
526
 
            /* -------------------------------------------------------------- */
527
 
            /* determine the layout of the new element */
528
 
            /* -------------------------------------------------------------- */
529
 
 
530
 
            csize = nrowsleft * ncolsleft ;
531
 
            size2 = UNITS (Element, 1)
532
 
                  + UNITS (Int, nrowsleft + ncolsleft)
533
 
                  + UNITS (Entry, csize) ;
534
 
 
535
 
            DEBUG7 (("Old size "ID" New size "ID"\n", size, size2)) ;
536
 
 
537
 
            pnext = pdest ;
538
 
            pnext->header.prevsize = size2 ;
539
 
            pdest -= (size2 + 1) ;
540
 
 
541
 
            ASSERT (size2 <= size) ;
542
 
            ASSERT ((unsigned Int) psrc + 1 + size <= (unsigned Int) pnext) ;
543
 
            ASSERT (psrc <= pdest) ;
544
 
 
545
 
            p = pdest + 1 ;
546
 
            epdest = (Element *) p ;
547
 
            p += UNITS (Element, 1) ;
548
 
            Cols2 = (Int *) p ;
549
 
            Rows2 = Cols2 + ncolsleft ;
550
 
            p += UNITS (Int, nrowsleft + ncolsleft) ;
551
 
            C2 = (Entry *) p ;
552
 
 
553
 
            ASSERT (epdest >= epsrc) ;
554
 
            ASSERT (Rows2 >= Rows) ;
555
 
            ASSERT (Cols2 >= Cols) ;
556
 
            ASSERT (C2 >= C) ;
557
 
            ASSERT (p + UNITS (Entry, csize) == pnext) ;
558
 
 
559
 
            /* -------------------------------------------------------------- */
560
 
            /* move the contribution block */
561
 
            /* -------------------------------------------------------------- */
562
 
 
563
 
            /* overlap = psrc + size + 1 > pdest ; */
564
 
 
565
 
            if (nrowsleft < nrows || ncolsleft < ncols)
566
 
            {
567
 
 
568
 
                /* ---------------------------------------------------------- */
569
 
                /* compress contribution block in place prior to moving it */
570
 
                /* ---------------------------------------------------------- */
571
 
 
572
 
                DEBUG7 (("Compress C in place prior to move:\n"));
573
 
#ifndef NDEBUG
574
 
                UMF_dump_dense (C, nrows, nrows, ncols) ;
575
 
#endif
576
 
                C1 = C ;
577
 
                C3 = C ;
578
 
                for (j = 0 ; j < ncols ; j++)
579
 
                {
580
 
                    if (Cols [j] >= 0)
581
 
                    {
582
 
                        for (i = 0 ; i < nrows ; i++)
583
 
                        {
584
 
                            if (Rows [i] >= 0)
585
 
                            {
586
 
                                *C3++ = C1 [i] ;
587
 
                            }
588
 
                        }
589
 
                    }
590
 
                    C1 += nrows ;
591
 
                }
592
 
                ASSERT (C3-C == csize) ;
593
 
                DEBUG8 (("Newly compressed contrib. block (all in use):\n")) ;
594
 
#ifndef NDEBUG
595
 
                UMF_dump_dense (C, nrowsleft, nrowsleft, ncolsleft) ;
596
 
#endif
597
 
            }
598
 
 
599
 
            /* shift the contribution block down */
600
 
            C += csize ;
601
 
            C2 += csize ;
602
 
            for (i = 0 ; i < csize ; i++)
603
 
            {
604
 
                *--C2 = *--C ;
605
 
            }
606
 
 
607
 
            /* -------------------------------------------------------------- */
608
 
            /* move the row indices */
609
 
            /* -------------------------------------------------------------- */
610
 
 
611
 
            i2 = nrowsleft ;
612
 
            for (i = nrows - 1 ; i >= 0 ; i--)
613
 
            {
614
 
                ASSERT (Rows2+i2 >= Rows+i) ;
615
 
                if (Rows [i] >= 0)
616
 
                {
617
 
                    Rows2 [--i2] = Rows [i] ;
618
 
                }
619
 
            }
620
 
            ASSERT (i2 == 0) ;
621
 
 
622
 
            j2 = ncolsleft ;
623
 
            for (j = ncols - 1 ; j >= 0 ; j--)
624
 
            {
625
 
                ASSERT (Cols2+j2 >= Cols+j) ;
626
 
                if (Cols [j] >= 0)
627
 
                {
628
 
                    Cols2 [--j2] = Cols [j] ;
629
 
                }
630
 
            }
631
 
            ASSERT (j2 == 0) ;
632
 
 
633
 
            /* -------------------------------------------------------------- */
634
 
            /* construct the new header */
635
 
            /* -------------------------------------------------------------- */
636
 
 
637
 
            /* E [0...e] is now valid */
638
 
            E [e] = (pdest + 1) - Numeric->Memory ;
639
 
            epdest = (Element *) (pdest + 1) ;
640
 
 
641
 
            epdest->next = EMPTY ;      /* destroys the son list */
642
 
            epdest->ncols = ncolsleft ;
643
 
            epdest->nrows = nrowsleft ;
644
 
            epdest->ncolsleft = ncolsleft ;
645
 
            epdest->nrowsleft = nrowsleft ;
646
 
            epdest->rdeg = rdeg ;
647
 
            epdest->cdeg = cdeg ;
648
 
 
649
 
            ASSERT (size2 <= size) ;
650
 
            pdest->header.prevsize = 0 ;
651
 
            pdest->header.size = size2 ;
652
 
 
653
 
            DEBUG7 (("After moving it:\n")) ;
654
 
#ifndef NDEBUG
655
 
            UMF_dump_element (Numeric, Work, e, FALSE) ;
656
 
#endif
657
 
        }
658
 
 
659
 
#ifndef NDEBUG
660
 
        else
661
 
        {
662
 
            DEBUG8 ((" free\n")) ;
663
 
        }
664
 
#endif
665
 
        DEBUG7 (("psrc "ID"  tail "ID"\n",
666
 
        (Int) (psrc-Numeric->Memory), Numeric->itail)) ;
667
 
    }
668
 
 
669
 
    ASSERT (psrc == Numeric->Memory + Numeric->itail) ;
670
 
    ASSERT (nmark == 0) ;
671
 
 
672
 
    /* ---------------------------------------------------------------------- */
673
 
    /* final tail pointer */
674
 
    /* ---------------------------------------------------------------------- */
675
 
 
676
 
    ASSERT (pdest >= Numeric->Memory + Numeric->itail) ;
677
 
    Numeric->itail = pdest - Numeric->Memory ;
678
 
    pdest->header.prevsize = 0 ;
679
 
    Numeric->ibig = EMPTY ;
680
 
    Numeric->tail_usage = Numeric->size - Numeric->itail ;
681
 
 
682
 
    /* ---------------------------------------------------------------------- */
683
 
    /* clear the unused E [nel+1 .. Work->elen - 1] */
684
 
    /* ---------------------------------------------------------------------- */
685
 
 
686
 
    for (e = nel+1 ; e < Work->elen ; e++)
687
 
    {
688
 
        E [e] = 0 ;
689
 
    }
690
 
 
691
 
#ifndef NDEBUG
692
 
    UMF_dump_packed_memory (Numeric, Work) ;
693
 
#endif
694
 
 
695
 
    DEBUG8 (("::::GARBAGE COLLECTION DONE::::\n")) ;
696
 
}