30 Int *
P, *Q, *R, *Pnum, *Offp, *Offi, *Pblock, *Pinv, *Iwork,
31 *Lip, *Uip, *Llen, *Ulen ;
32 Entry *Offx, *X, s, *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,
48 nblocks = Symbolic->nblocks ;
49 nzoff = Symbolic->nzoff ;
51 Pnum = Numeric->Pnum ;
52 Offp = Numeric->Offp ;
53 Offi = Numeric->Offi ;
54 Offx = (
Entry *) Numeric->Offx ;
58 Llen = Numeric->Llen ;
59 Ulen = Numeric->Ulen ;
60 LUbx = (
Unit **) Numeric->LUbx ;
61 Udiag = Numeric->Udiag ;
64 Pinv = Numeric->Pinv ;
65 X = (
Entry *) Numeric->Xwork ;
66 Iwork = Numeric->Iwork ;
68 Pblock = Iwork + 5*((
size_t) Symbolic->maxblock) ;
78 for (k = 0 ; k <
n ; k++)
83 for (k = 0 ; k <
n ; k++)
121 for (k = 0 ; k <
n ; k++)
PRINTF ((
"Rs [%d] %g\n", k, Rs [k])) ;
129 for (block = 0 ; block < nblocks ; block++)
139 PRINTF ((
"FACTOR BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1,k2-1,nk)) ;
150 pend = Ap [oldcol+1] ;
156 for (p = Ap [oldcol] ; p < pend ; p++)
159 newrow = Pinv [oldrow] ;
162 Offi [poff] = oldrow ;
163 Offx [poff] = Ax [p] ;
169 PRINTF ((
"singleton block %d", block)) ;
182 for (p = Ap [oldcol] ; p < pend ; p++)
185 newrow = Pinv [oldrow] ;
188 Offi [poff] = oldrow ;
196 PRINTF ((
"singleton block %d ", block)) ;
209 Common->numerical_rank = k1 ;
210 Common->singular_col = oldcol ;
211 if (
Common->halt_if_singular)
234 lsize = -(
Common->initmem) ;
238 lsize =
Common->initmem_amd * Lnz [block] + nk ;
243 lsize, &LUbx [block], Udiag + k1, Llen + k1, Ulen + k1,
244 Lip + k1, Uip + k1, Pblock, &lnz_block, &unz_block,
245 X, Iwork, k1, Pinv, Rs, Offp, Offi, Offx,
Common) ;
254 PRINTF ((
"\n----------------------- L %d:\n", block)) ;
256 PRINTF ((
"\n----------------------- U %d:\n", block)) ;
265 max_lnz_block =
MAX (max_lnz_block, lnz_block) ;
266 max_unz_block =
MAX (max_unz_block, unz_block) ;
268 if (Lnz [block] ==
EMPTY)
271 Lnz [block] =
MAX (lnz_block, unz_block) ;
278 PRINTF ((
"Pnum, 1-based:\n")) ;
279 for (k = 0 ; k < nk ; k++)
283 Pnum [k + k1] =
P [Pblock [k] + k1] ;
284 PRINTF ((
"Pnum (%d + %d + 1 = %d) = %d + 1 = %d\n",
285 k, k1, k+k1+1, Pnum [k+k1], Pnum [k+k1]+1)) ;
292 PRINTF ((
"\n------------------- Off diagonal entries:\n")) ;
297 Numeric->max_lnz_block = max_lnz_block ;
298 Numeric->max_unz_block = max_unz_block ;
302 for (k = 0 ; k <
n ; k++)
307 for (k = 0 ; k <
n ; k++)
309 ASSERT (Pnum [k] >= 0 && Pnum [k] <
n) ;
310 Pinv [Pnum [k]] = k ;
319 for (k = 0 ; k <
n ; k++)
321 REAL (X [k]) = Rs [Pnum [k]] ;
323 for (k = 0 ; k <
n ; k++)
325 Rs [k] =
REAL (X [k]) ;
329 PRINTF ((
"\n------------------- Off diagonal entries, old:\n")) ;
333 for (p = 0 ; p < nzoff ; p++)
335 ASSERT (Offi [p] >= 0 && Offi [p] <
n) ;
336 Offi [p] = Pinv [Offi [p]] ;
339 PRINTF ((
"\n------------------- Off diagonal entries, new:\n")) ;
344 PRINTF ((
"\n ############# KLU_BTF_FACTOR done, nblocks %d\n",nblocks));
345 Entry ss, *Udiag = Numeric->Udiag ;
346 for (block = 0 ; block < nblocks &&
Common->status ==
KLU_OK ; block++)
351 PRINTF ((
"\n======================KLU_factor output: k1 %d k2 %d nk %d\n",k1,k2,nk)) ;
361 Int *Lip, *Uip, *Llen, *Ulen ;
363 Lip = Numeric->Lip + k1 ;
364 Llen = Numeric->Llen + k1 ;
365 LU = (
Unit *) Numeric->LUbx [block] ;
366 PRINTF ((
"\n---- L block %d\n", block));
368 Uip = Numeric->Uip + k1 ;
369 Ulen = Numeric->Ulen + k1 ;
370 PRINTF ((
"\n---- U block %d\n", block)) ;
396 Int n, nzoff, nblocks, maxblock, k, ok =
TRUE ;
399 size_t n1, nzoff1, s, b6, n3 ;
400 #ifdef KLU_ENABLE_OPENMP 417 if (Symbolic ==
NULL)
424 nzoff = Symbolic->nzoff ;
425 nblocks = Symbolic->nblocks ;
426 maxblock = Symbolic->maxblock ;
428 PRINTF ((
"KLU_factor: n %d nzoff %d nblocks %d maxblock %d\n",
429 n, nzoff, nblocks, maxblock)) ;
446 n1 = ((size_t)
n) + 1 ;
447 nzoff1 = ((size_t) nzoff) + 1 ;
457 Numeric->nblocks = nblocks ;
458 Numeric->nzoff = nzoff ;
472 if (Numeric->LUbx !=
NULL)
474 for (k = 0 ; k < nblocks ; k++)
476 Numeric->LUbx [k] =
NULL ;
503 #ifdef KLU_ENABLE_OPENMP 505 num_threads = omp_get_num_threads();
515 Numeric->Xwork = Numeric->Work ;
516 Numeric->Iwork = (
Int *) ((
Entry *) Numeric->Xwork +
n) ;
542 if (
Common->halt_if_singular)
#define KLU_kernel_factor
KLU_numeric * KLU_factor(Int Ap [], Int Ai [], double Ax [], KLU_symbolic *Symbolic, KLU_common *Common)
#define KLU_OUT_OF_MEMORY
#define ASSERT(expression)
static void factor2(Int Ap [], Int Ai [], Entry Ax [], KLU_symbolic *Symbolic, KLU_numeric *Numeric, KLU_common *Common)
#define SCALE_DIV_ASSIGN(a, c, s)