← Back to branch summary

~ubuntu-branches/ubuntu/natty/suitesparse/natty

« back to all changes in this revision

Viewing changes to KLU/Source/klu_factor.c

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2006-12-22 10:16:15 UTC
  • Revision ID: james.westby@ubuntu.com-20061222101615-2ohaj8902oix2rnk
Tags: upstream-2.3.1
Import upstream version 2.3.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
KLU/Source/klu_factor.c
 
1
/* ========================================================================== */
 
2
/* === klu_factor =========================================================== */
 
3
/* ========================================================================== */
 
4
 
 
5
/* Factor the matrix, after ordering and analyzing it with klu_analyze
 
6
 * or klu_analyze_given.
 
7
 */
 
8
 
 
9
#include "klu_internal.h"
 
10
 
 
11
/* ========================================================================== */
 
12
/* === klu_factor2 ========================================================== */
 
13
/* ========================================================================== */
 
14
 
 
15
static void klu_factor2
 
16
(
 
17
    /* inputs, not modified */
 
18
    int Ap [ ],         /* size n+1, column pointers */
 
19
    int Ai [ ],         /* size nz, row indices */
 
20
    Entry Ax [ ],
 
21
    klu_symbolic *Symbolic,
 
22
 
 
23
    /* inputs, modified on output: */
 
24
    klu_numeric *Numeric,
 
25
    klu_common *Common
 
26
)
 
27
{
 
28
    double lsize ;
 
29
    double *Lnz, *Rs ;
 
30
    int *P, *Q, *R, *Pnum, *Offp, *Offi, *Pblock, *Pinv, *Iwork ;
 
31
    int **Lbip, **Ubip, **Lblen, **Ublen ;
 
32
    Entry *Offx, *Singleton, *X, s ;
 
33
    Unit **LUbx, **Udiag ;
 
34
    int k1, k2, nk, k, block, oldcol, pend, oldrow, n, lnz, unz, p, newrow,
 
35
        nblocks, poff, nzoff, lnz_block, unz_block, scale, max_lnz_block,
 
36
        max_unz_block ;
 
37
 
 
38
    /* ---------------------------------------------------------------------- */
 
39
    /* initializations */
 
40
    /* ---------------------------------------------------------------------- */
 
41
 
 
42
    /* get the contents of the Symbolic object */
 
43
    n = Symbolic->n ;
 
44
    P = Symbolic->P ;
 
45
    Q = Symbolic->Q ;
 
46
    R = Symbolic->R ;
 
47
    Lnz = Symbolic->Lnz ;
 
48
    nblocks = Symbolic->nblocks ;
 
49
    nzoff = Symbolic->nzoff ;
 
50
 
 
51
    Pnum = Numeric->Pnum ;
 
52
    Offp = Numeric->Offp ;
 
53
    Offi = Numeric->Offi ;
 
54
    Offx = (Entry *) Numeric->Offx ;
 
55
    Singleton = (Entry *) Numeric->Singleton ;
 
56
 
 
57
    Lbip = Numeric->Lbip ;
 
58
    Ubip = Numeric->Ubip ;
 
59
    Lblen = Numeric->Lblen ;
 
60
    Ublen = Numeric->Ublen ;
 
61
    LUbx = (Unit **) Numeric->LUbx ;
 
62
    Udiag = (Unit **) Numeric->Udiag ;
 
63
 
 
64
    Rs = Numeric->Rs ;
 
65
    Pinv = Numeric->Pinv ;
 
66
    X = (Entry *) Numeric->Xwork ;              /* X is of size n */
 
67
    Iwork = Numeric->Iwork ;                    /* 5*maxblock for klu_factor */
 
68
                                                /* 1*maxblock for Pblock */
 
69
    Pblock = Iwork + 5*((size_t) Symbolic->maxblock) ;
 
70
    Common->nrealloc = 0 ;
 
71
    scale = Common->scale ;
 
72
    max_lnz_block = 1 ;
 
73
    max_unz_block = 1 ;
 
74
 
 
75
    /* compute the inverse of P from symbolic analysis.  Will be updated to
 
76
     * become the inverse of the numerical factorization when the factorization
 
77
     * is done, for use in klu_refactor */
 
78
#ifndef NDEBUG
 
79
    for (k = 0 ; k < n ; k++)
 
80
    {
 
81
        Pinv [k] = EMPTY ;
 
82
    }
 
83
#endif
 
84
    for (k = 0 ; k < n ; k++)
 
85
    {
 
86
        ASSERT (P [k] >= 0 && P [k] < n) ;
 
87
        Pinv [P [k]] = k ;
 
88
    }
 
89
#ifndef NDEBUG
 
90
    for (k = 0 ; k < n ; k++) ASSERT (Pinv [k] != EMPTY) ;
 
91
#endif
 
92
 
 
93
    for (block = 0 ; block < nblocks ; block++)
 
94
    {
 
95
        /* Singleton [block] = 0 ; */
 
96
        CLEAR (Singleton [block]) ;
 
97
    }
 
98
 
 
99
    lnz = 0 ;
 
100
    unz = 0 ;
 
101
    Common->noffdiag = 0 ;
 
102
    Offp [0] = 0 ;
 
103
 
 
104
    /* ---------------------------------------------------------------------- */
 
105
    /* optionally check input matrix and compute scale factors */
 
106
    /* ---------------------------------------------------------------------- */
 
107
 
 
108
    if (scale >= 0)
 
109
    {
 
110
        /* use Pnum as workspace */
 
111
        KLU_scale (scale, n, Ap, Ai, (double *) Ax, Rs, Pnum, Common) ;
 
112
        if (Common->status < KLU_OK)
 
113
        {
 
114
            /* matrix is invalid */
 
115
            return ;
 
116
        }
 
117
    }
 
118
 
 
119
#ifndef NDEBUG
 
120
    if (scale > 0)
 
121
    {
 
122
        for (k = 0 ; k < n ; k++) PRINTF (("Rs [%d] %g\n", k, Rs [k])) ;
 
123
    }
 
124
#endif
 
125
 
 
126
    /* ---------------------------------------------------------------------- */
 
127
    /* factor each block using klu */
 
128
    /* ---------------------------------------------------------------------- */
 
129
 
 
130
    for (block = 0 ; block < nblocks ; block++)
 
131
    {
 
132
 
 
133
        /* ------------------------------------------------------------------ */
 
134
        /* the block is from rows/columns k1 to k2-1 */
 
135
        /* ------------------------------------------------------------------ */
 
136
 
 
137
        k1 = R [block] ;
 
138
        k2 = R [block+1] ;
 
139
        nk = k2 - k1 ;
 
140
        PRINTF (("FACTOR BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1,k2-1,nk)) ;
 
141
 
 
142
        if (nk == 1)
 
143
        {
 
144
 
 
145
            /* -------------------------------------------------------------- */
 
146
            /* singleton case */
 
147
            /* -------------------------------------------------------------- */
 
148
 
 
149
            poff = Offp [k1] ;
 
150
            oldcol = Q [k1] ;
 
151
            pend = Ap [oldcol+1] ;
 
152
            CLEAR (s) ;
 
153
 
 
154
            if (scale <= 0)
 
155
            {
 
156
                /* no scaling */
 
157
                for (p = Ap [oldcol] ; p < pend ; p++)
 
158
                {
 
159
                    oldrow = Ai [p] ;
 
160
                    newrow = Pinv [oldrow] ;
 
161
                    if (newrow < k1)
 
162
                    {
 
163
                        Offi [poff] = oldrow ;
 
164
                        Offx [poff] = Ax [p] ;
 
165
                        poff++ ;
 
166
                    }
 
167
                    else
 
168
                    {
 
169
                        ASSERT (newrow == k1) ;
 
170
                        PRINTF (("Singleton block %d", block)) ;
 
171
                        PRINT_ENTRY (Ax [p]) ;
 
172
                        s = Ax [p] ;
 
173
                    }
 
174
                }
 
175
            }
 
176
            else
 
177
            {
 
178
                /* row scaling */
 
179
                for (p = Ap [oldcol] ; p < pend ; p++)
 
180
                {
 
181
                    oldrow = Ai [p] ;
 
182
                    newrow = Pinv [oldrow] ;
 
183
                    if (newrow < k1)
 
184
                    {
 
185
                        Offi [poff] = oldrow ;
 
186
                        /*Offx [poff] = Ax [p] / Rs [oldrow] ; */
 
187
                        SCALE_DIV_ASSIGN (Offx [poff], Ax [p], Rs [oldrow]) ;
 
188
                        poff++ ;
 
189
                    }
 
190
                    else
 
191
                    {
 
192
                        ASSERT (newrow == k1) ;
 
193
                        PRINTF (("Singleton block %d ", block)) ;
 
194
                        PRINT_ENTRY (Ax[p]) ;
 
195
                        SCALE_DIV_ASSIGN (s, Ax [p], Rs [oldrow]) ;
 
196
                    }
 
197
                }
 
198
            }
 
199
 
 
200
            Singleton [block] = s ;
 
201
 
 
202
            if (IS_ZERO (s))
 
203
            {
 
204
                /* singular singleton */
 
205
                Common->status = KLU_SINGULAR ;
 
206
                Common->numerical_rank = k1 ;
 
207
                Common->singular_col = oldcol ;
 
208
                if (Common->halt_if_singular)
 
209
                {
 
210
                    return ;
 
211
                }
 
212
            }
 
213
 
 
214
            Offp [k1+1] = poff ;
 
215
            Pnum [k1] = P [k1] ;
 
216
            lnz++ ;
 
217
            unz++ ;
 
218
 
 
219
        }
 
220
        else
 
221
        {
 
222
 
 
223
            /* -------------------------------------------------------------- */
 
224
            /* construct and factorize the kth block */
 
225
            /* -------------------------------------------------------------- */
 
226
 
 
227
            if (Lnz [block] < 0)
 
228
            {
 
229
                /* COLAMD was used - no estimate of fill-in */
 
230
                /* use 10 times the nnz in A, plus n */
 
231
                lsize = -(Common->initmem) ;
 
232
            }
 
233
            else
 
234
            {
 
235
                lsize = Common->initmem_amd * Lnz [block] + nk ;
 
236
            }
 
237
 
 
238
            /* allocates 5 arrays:
 
239
             * Lbip [block], Ubip [block], Lblen [block], Ublen [block],
 
240
             * LUbx [block] */
 
241
            KLU_kernel_factor (nk, Ap, Ai, Ax, Q, lsize,
 
242
                    &LUbx [block], Udiag [block], Lblen [block], Ublen [block],
 
243
                    Lbip [block], Ubip [block], Pblock, &lnz_block, &unz_block,
 
244
                    X, Iwork,
 
245
                    /* BTF and scale-related arguments: */
 
246
                    k1, Pinv, Rs, Offp, Offi, Offx, Common) ;
 
247
 
 
248
            if (Common->status < KLU_OK ||
 
249
               (Common->status == KLU_SINGULAR && Common->halt_if_singular))
 
250
            {
 
251
                /* out of memory, invalid inputs, or singular */
 
252
                return ;
 
253
            }
 
254
 
 
255
            PRINTF (("\n----------------------- L %d:\n", block)) ;
 
256
            ASSERT (KLU_valid_LU (nk, TRUE, Lbip [block],
 
257
                    Lblen [block], LUbx [block])) ;
 
258
            PRINTF (("\n----------------------- U %d:\n", block)) ;
 
259
            ASSERT (KLU_valid_LU (nk, FALSE, Ubip [block],
 
260
                    Ublen [block], LUbx [block])) ;
 
261
 
 
262
            /* -------------------------------------------------------------- */
 
263
            /* get statistics */
 
264
            /* -------------------------------------------------------------- */
 
265
 
 
266
            lnz += lnz_block ;
 
267
            unz += unz_block ;
 
268
            max_lnz_block = MAX (max_lnz_block, lnz_block) ;
 
269
            max_unz_block = MAX (max_unz_block, unz_block) ;
 
270
 
 
271
            if (Lnz [block] == EMPTY)
 
272
            {
 
273
                /* revise estimate for subsequent factorization */
 
274
                Lnz [block] = MAX (lnz_block, unz_block) ;
 
275
            }
 
276
 
 
277
            /* -------------------------------------------------------------- */
 
278
            /* combine the klu row ordering with the symbolic pre-ordering */
 
279
            /* -------------------------------------------------------------- */
 
280
 
 
281
            PRINTF (("Pnum, 1-based:\n")) ;
 
282
            for (k = 0 ; k < nk ; k++)
 
283
            {
 
284
                ASSERT (k + k1 < n) ;
 
285
                ASSERT (Pblock [k] + k1 < n) ;
 
286
                Pnum [k + k1] = P [Pblock [k] + k1] ;
 
287
                PRINTF (("Pnum (%d + %d + 1 = %d) = %d + 1 = %d\n",
 
288
                    k, k1, k+k1+1, Pnum [k+k1], Pnum [k+k1]+1)) ;
 
289
            }
 
290
 
 
291
            /* the local pivot row permutation Pblock is no longer needed */
 
292
        }
 
293
    }
 
294
    ASSERT (nzoff == Offp [n]) ;
 
295
    PRINTF (("\n------------------- Off diagonal entries:\n")) ;
 
296
    ASSERT (KLU_valid (n, Offp, Offi, Offx)) ;
 
297
 
 
298
    Numeric->lnz = lnz ;
 
299
    Numeric->unz = unz ;
 
300
    Numeric->max_lnz_block = max_lnz_block ;
 
301
    Numeric->max_unz_block = max_unz_block ;
 
302
 
 
303
    /* Numeric->flops = EMPTY ;         TODO not yet computed */
 
304
 
 
305
    /* compute the inverse of Pnum */
 
306
#ifndef NDEBUG
 
307
    for (k = 0 ; k < n ; k++)
 
308
    {
 
309
        Pinv [k] = EMPTY ;
 
310
    }
 
311
#endif
 
312
    for (k = 0 ; k < n ; k++)
 
313
    {
 
314
        ASSERT (Pnum [k] >= 0 && Pnum [k] < n) ;
 
315
        Pinv [Pnum [k]] = k ;
 
316
    }
 
317
#ifndef NDEBUG
 
318
    for (k = 0 ; k < n ; k++) ASSERT (Pinv [k] != EMPTY) ;
 
319
#endif
 
320
 
 
321
    /* permute scale factors Rs according to pivotal row order */
 
322
    if (scale > 0)
 
323
    {
 
324
        for (k = 0 ; k < n ; k++)
 
325
        {
 
326
            REAL (X [k]) = Rs [Pnum [k]] ;
 
327
        }
 
328
        for (k = 0 ; k < n ; k++)
 
329
        {
 
330
            Rs [k] = REAL (X [k]) ;
 
331
        }
 
332
    }
 
333
 
 
334
    PRINTF (("\n------------------- Off diagonal entries, old:\n")) ;
 
335
    ASSERT (KLU_valid (n, Offp, Offi, Offx)) ;
 
336
 
 
337
    /* apply the pivot row permutations to the off-diagonal entries */
 
338
    for (p = 0 ; p < nzoff ; p++)
 
339
    {
 
340
        ASSERT (Offi [p] >= 0 && Offi [p] < n) ;
 
341
        Offi [p] = Pinv [Offi [p]] ;
 
342
    }
 
343
 
 
344
    PRINTF (("\n------------------- Off diagonal entries, new:\n")) ;
 
345
    ASSERT (KLU_valid (n, Offp, Offi, Offx)) ;
 
346
 
 
347
#ifndef NDEBUG
 
348
    {
 
349
        PRINTF (("\n ############# KLU_BTF_FACTOR done, nblocks %d\n",nblocks));
 
350
        Entry *singleton = Numeric->Singleton ;
 
351
        for (block = 0 ; block < nblocks && Common->status == KLU_OK ; block++)
 
352
        {
 
353
            k1 = R [block] ;
 
354
            k2 = R [block+1] ;
 
355
            nk = k2 - k1 ;
 
356
            PRINTF (("\n======================klu_factor output: k1 %d k2 %d nk %d\n",k1,k2,nk)) ;
 
357
            if (nk == 1)
 
358
            {
 
359
                PRINTF (("singleton  ")) ;
 
360
                /* ENTRY_PRINT (singleton [block]) ; */
 
361
                PRINT_ENTRY (singleton [block]) ;
 
362
            }
 
363
            else
 
364
            {
 
365
                int *Lip, *Uip, *Llen, *Ulen ;
 
366
                Unit *LU ;
 
367
                Lip = Lbip [block] ;
 
368
                Llen = Lblen [block] ;
 
369
                LU = (Unit *) Numeric->LUbx [block] ;
 
370
                PRINTF (("\n---- L block %d\n", block));
 
371
                ASSERT (KLU_valid_LU (nk, TRUE, Lip, Llen, LU)) ;
 
372
                Uip = Ubip [block] ;
 
373
                Ulen = Ublen [block] ;
 
374
                PRINTF (("\n---- U block %d\n", block)) ;
 
375
                ASSERT (KLU_valid_LU (nk, FALSE, Uip, Ulen, LU)) ;
 
376
            }
 
377
        }
 
378
    }
 
379
#endif
 
380
}
 
381
 
 
382
 
 
383
/* ========================================================================== */
 
384
/* === CLEAR_PTR ============================================================ */
 
385
/* ========================================================================== */
 
386
 
 
387
/* Set an array of pointers to NULL */
 
388
 
 
389
#define CLEAR_PTR(Ptr,size) \
 
390
{ \
 
391
    int ii ; \
 
392
    if (Ptr != NULL) \
 
393
    { \
 
394
        for (ii = 0 ; ii < size ; ii++) \
 
395
        { \
 
396
            Ptr [ii] = NULL ; \
 
397
        } \
 
398
    } \
 
399
}
 
400
 
 
401
 
 
402
/* ========================================================================== */
 
403
/* === klu_factor =========================================================== */
 
404
/* ========================================================================== */
 
405
 
 
406
klu_numeric *KLU_factor         /* returns NULL if error, or a valid
 
407
                                   klu_numeric object if successful */
 
408
(
 
409
    /* --- inputs --- */
 
410
    int Ap [ ],         /* size n+1, column pointers */
 
411
    int Ai [ ],         /* size nz, row indices */
 
412
    double Ax [ ],
 
413
    klu_symbolic *Symbolic,
 
414
    /* -------------- */
 
415
    klu_common *Common
 
416
)
 
417
{
 
418
    int n, nzoff, nblocks, maxblock, block, k1, k2, nk, ok = TRUE ;
 
419
    int *R ;
 
420
    int **Lbip, **Ubip, **Lblen, **Ublen ;
 
421
    klu_numeric *Numeric ;
 
422
    Unit **Udiag ;
 
423
    size_t n1, nzoff1, nunits, s, b6, n3, nk1 ;
 
424
 
 
425
    if (Common == NULL)
 
426
    {
 
427
        return (NULL) ;
 
428
    }
 
429
    Common->status = KLU_OK ;
 
430
    Common->numerical_rank = EMPTY ;
 
431
    Common->singular_col = EMPTY ;
 
432
 
 
433
    /* ---------------------------------------------------------------------- */
 
434
    /* get the contents of the Symbolic object */
 
435
    /* ---------------------------------------------------------------------- */
 
436
 
 
437
    /* check for a valid Symbolic object */
 
438
    if (Symbolic == NULL)
 
439
    {
 
440
        Common->status = KLU_INVALID ;
 
441
        return (NULL) ;
 
442
    }
 
443
 
 
444
    n = Symbolic->n ;
 
445
    nzoff = Symbolic->nzoff ;
 
446
    nblocks = Symbolic->nblocks ;
 
447
    maxblock = Symbolic->maxblock ;
 
448
    R = Symbolic->R ;
 
449
    PRINTF (("klu_factor:  n %d nzoff %d nblocks %d maxblock %d\n",
 
450
        n, nzoff, nblocks, maxblock)) ;
 
451
 
 
452
    /* ---------------------------------------------------------------------- */
 
453
    /* get control parameters and make sure they are in the proper range */
 
454
    /* ---------------------------------------------------------------------- */
 
455
 
 
456
    Common->initmem_amd = MAX (1.0, Common->initmem_amd) ;
 
457
    Common->initmem = MAX (1.0, Common->initmem) ;
 
458
    Common->tol = MIN (Common->tol, 1.0) ;
 
459
    Common->tol = MAX (0.0, Common->tol) ;
 
460
    Common->growth = MAX (1.0, Common->growth) ;
 
461
 
 
462
    /* ---------------------------------------------------------------------- */
 
463
    /* allocate the Numeric object  */
 
464
    /* ---------------------------------------------------------------------- */
 
465
 
 
466
    /* this will not cause size_t overflow (already checked by klu_symbolic) */
 
467
    n1 = ((size_t) n) + 1 ;
 
468
    nzoff1 = ((size_t) nzoff) + 1 ;
 
469
 
 
470
    Numeric = klu_malloc (sizeof (klu_numeric), 1, Common) ;
 
471
    if (Common->status < KLU_OK)
 
472
    {
 
473
        /* out of memory */
 
474
        Common->status = KLU_OUT_OF_MEMORY ;
 
475
        return (NULL) ;
 
476
    }
 
477
    Numeric->nblocks = nblocks ;
 
478
    Numeric->Pnum = klu_malloc (n, sizeof (int), Common) ;
 
479
    Numeric->Offp = klu_malloc (n1, sizeof (int), Common) ;
 
480
    Numeric->Offi = klu_malloc (nzoff1, sizeof (int), Common) ;
 
481
    Numeric->Offx = klu_malloc (nzoff1, sizeof (Entry), Common) ;
 
482
    Numeric->Singleton = klu_malloc (nblocks, sizeof (Entry), Common) ;
 
483
 
 
484
    Numeric->Lbip  = klu_malloc (nblocks, sizeof (int *), Common) ;
 
485
    Numeric->Ubip  = klu_malloc (nblocks, sizeof (int *), Common) ;
 
486
    Numeric->Lblen = klu_malloc (nblocks, sizeof (int *), Common) ;
 
487
    Numeric->Ublen = klu_malloc (nblocks, sizeof (int *), Common) ;
 
488
    Numeric->LUbx  = klu_malloc (nblocks, sizeof (Unit *), Common) ;
 
489
    Numeric->Udiag = klu_malloc (nblocks, sizeof (Unit *), Common) ;
 
490
 
 
491
    if (Common->scale > 0)
 
492
    {
 
493
        Numeric->Rs = klu_malloc (n, sizeof (double), Common) ;
 
494
    }
 
495
    else
 
496
    {
 
497
        /* no scaling */
 
498
        Numeric->Rs = NULL ;
 
499
    }
 
500
 
 
501
    Numeric->Pinv = klu_malloc (n, sizeof (int), Common) ;
 
502
 
 
503
    /* allocate permanent workspace for factorization and solve.
 
504
     * Note that the solver will use an Xwork of size 4n, whereas
 
505
     * the factorization codes use an Xwork of size n and integer space
 
506
     * (Iwork) of size 6n. klu_condest uses an Xwork of size 2n. */
 
507
 
 
508
    s = klu_mult_size_t (n, sizeof (Entry), &ok) ;
 
509
    n3 = klu_mult_size_t (n, 3 * sizeof (Entry), &ok) ;
 
510
    b6 = klu_mult_size_t (maxblock, 6 * sizeof (int), &ok) ;
 
511
    Numeric->worksize = klu_add_size_t (s, MAX (n3, b6), &ok) ;
 
512
    if (!ok)
 
513
    {
 
514
        /* problem too large (almost impossible to happen here) */
 
515
        Common->status = KLU_TOO_LARGE ;
 
516
        KLU_free_numeric (&Numeric, Common) ;
 
517
        return (NULL) ;
 
518
    }
 
519
 
 
520
/*
 
521
    Numeric->worksize = n * sizeof (Entry) +
 
522
        MAX (n * 3 *sizeof (Entry), 6*maxblock * sizeof (int)) ;
 
523
*/
 
524
 
 
525
    Numeric->Work = klu_malloc (Numeric->worksize, 1, Common) ;
 
526
    Numeric->Xwork = Numeric->Work ;
 
527
    Numeric->Iwork = (int *) ((Entry *) Numeric->Xwork + n) ;
 
528
 
 
529
    if (Common->status < KLU_OK)
 
530
    {
 
531
        /* out of memory */
 
532
        KLU_free_numeric (&Numeric, Common) ;
 
533
        return (NULL) ;
 
534
    }
 
535
 
 
536
    /* clear the pointer arrays, so that klu_free_numeric works OK */
 
537
    CLEAR_PTR (Numeric->Lbip,  nblocks) ;
 
538
    CLEAR_PTR (Numeric->Ubip,  nblocks) ;
 
539
    CLEAR_PTR (Numeric->Lblen, nblocks) ;
 
540
    CLEAR_PTR (Numeric->Ublen, nblocks) ;
 
541
    CLEAR_PTR (Numeric->LUbx,  nblocks) ;
 
542
    CLEAR_PTR (Numeric->Udiag, nblocks) ;
 
543
 
 
544
    /* allocate the column pointer arrays for each block */
 
545
    Lbip = Numeric->Lbip ;
 
546
    Ubip = Numeric->Ubip ;
 
547
    Lblen = Numeric->Lblen ;
 
548
    Ublen = Numeric->Ublen ;
 
549
    Udiag = (Unit **) Numeric->Udiag ;
 
550
 
 
551
    for (block = 0 ; block < nblocks ; block++)
 
552
    {
 
553
        k1 = R [block] ;
 
554
        k2 = R [block+1] ;
 
555
        nk = k2 - k1 ;
 
556
        nk1 = ((size_t) nk) + 1 ;   /* cannot overflow */
 
557
        if (nk > 1)
 
558
        {
 
559
            nunits = UNITS (Entry, nk) ;
 
560
            Lbip [block]  = klu_malloc (nk1, sizeof (int), Common) ;
 
561
            Ubip [block]  = klu_malloc (nk1, sizeof (int), Common) ;
 
562
            Lblen [block] = klu_malloc (nk1, sizeof (int), Common) ;
 
563
            Ublen [block] = klu_malloc (nk1, sizeof (int), Common) ;
 
564
            Udiag [block] = klu_malloc (nunits, sizeof (Unit), Common) ;
 
565
            if (Common->status < KLU_OK)
 
566
            {
 
567
                /* out of memory */
 
568
                KLU_free_numeric (&Numeric, Common) ;
 
569
                return (NULL) ;
 
570
            }
 
571
        }
 
572
    }
 
573
 
 
574
    /* ---------------------------------------------------------------------- */
 
575
    /* factorize the blocks */
 
576
    /* ---------------------------------------------------------------------- */
 
577
 
 
578
    klu_factor2 (Ap, Ai, (Entry *) Ax, Symbolic, Numeric, Common) ;
 
579
 
 
580
    /* ---------------------------------------------------------------------- */
 
581
    /* return or free the Numeric object */
 
582
    /* ---------------------------------------------------------------------- */
 
583
 
 
584
    if (Common->status < KLU_OK)
 
585
    {
 
586
        /* out of memory or inputs invalid */
 
587
        KLU_free_numeric (&Numeric, Common) ;
 
588
    }
 
589
    else if (Common->status == KLU_SINGULAR)
 
590
    {
 
591
        if (Common->halt_if_singular)
 
592
        {
 
593
            /* Matrix is singular, and the Numeric object is only partially
 
594
             * defined because we halted early.  This is the default case for
 
595
             * a singular matrix. */
 
596
            KLU_free_numeric (&Numeric, Common) ;
 
597
        }
 
598
    }
 
599
    else if (Common->status == KLU_OK)
 
600
    {
 
601
        /* successful non-singular factorization */
 
602
        Common->numerical_rank = n ;
 
603
        Common->singular_col = n ;
 
604
    }
 
605
    return (Numeric) ;
 
606
}