43 Int i, pos, jnew, head, l_length ;
56 ASSERT (jnew >= 0 && jnew < k) ;
62 PRINTF ((
"[ start dfs at %d : new %d\n", j, jnew)) ;
65 (Lpend [jnew] ==
EMPTY) ? Llen [jnew] : Lpend [jnew] ;
70 Li = (
Int *) (LU + Lip [jnew]) ;
71 for (pos = --Ap_pos [head] ; pos >= 0 ; --pos)
107 PRINTF ((
" end dfs at %d ] head : %d\n", j, head)) ;
111 *plength = l_length ;
157 Int i, p, pend, oldcol, kglobal, top, l_length ;
161 Lik = (
Int *) (LU + lup);
168 oldcol = Q [kglobal] ;
169 pend = Ap [oldcol+1] ;
170 for (p = Ap [oldcol] ; p < pend ; p++)
172 i = PSinv [Ai [p]] - k1 ;
173 if (i < 0) continue ;
176 PRINTF ((
"\n ===== DFS at node %d in b, inew: %d\n", i, Pinv [i])) ;
181 top =
dfs (i, k, Pinv, Llen, Lip, Stack, Flag,
182 Lpend, top, LU, Lik, &l_length, Ap_pos) ;
195 Llen [k] = l_length ;
235 Int i, p, pend, oldcol, kglobal, poff, oldrow ;
242 poff = Offp [kglobal] ;
243 oldcol = Q [kglobal] ;
244 pend = Ap [oldcol+1] ;
249 for (p = Ap [oldcol] ; p < pend ; p++)
252 i = PSinv [oldrow] - k1 ;
257 Offi [poff] = oldrow ;
271 for (p = Ap [oldcol] ; p < pend ; p++)
274 i = PSinv [oldrow] - k1 ;
280 Offi [poff] = oldrow ;
292 Offp [kglobal+1] = poff ;
327 Int p, s, j, jnew, len ;
330 for (s = top ; s <
n ; s++)
338 ASSERT (Lip [jnew] <= Lip [jnew+1]) ;
339 for (p = 0 ; p < len ; p++)
375 Entry x, pivot, *Lx ;
376 double abs_pivot, xabs ;
377 Int p, i, ppivrow, pdiag, pivrow, *Li, last_row_index, firstrow, len ;
383 if (
Common->halt_if_singular)
387 for (firstrow = *p_firstrow ; firstrow <
n ; firstrow++)
389 PRINTF ((
"check %d\n", firstrow)) ;
390 if (Pinv [firstrow] < 0)
394 PRINTF ((
"Got pivotal row: %d\n", pivrow)) ;
398 ASSERT (pivrow >= 0 && pivrow <
n) ;
403 *p_firstrow = firstrow ;
412 last_row_index = Li [i] ;
419 for (p = 0 ; p < len ; p++)
437 if (xabs > abs_pivot)
445 ABS (xabs, X [last_row_index]) ;
446 if (xabs > abs_pivot)
453 if (last_row_index == diagrow)
455 if (xabs >=
tol * abs_pivot)
461 else if (pdiag !=
EMPTY)
464 ABS (xabs, Lx [pdiag]) ;
465 if (xabs >=
tol * abs_pivot)
473 if (ppivrow !=
EMPTY)
475 pivrow = Li [ppivrow] ;
476 pivot = Lx [ppivrow] ;
478 Li [ppivrow] = last_row_index ;
479 Lx [ppivrow] = X [last_row_index] ;
483 pivrow = last_row_index ;
484 pivot = X [last_row_index] ;
486 CLEAR (X [last_row_index]) ;
490 *p_abs_pivot = abs_pivot ;
491 ASSERT (pivrow >= 0 && pivrow <
n) ;
500 for (p = 0 ; p < Llen [k] ; p++)
503 DIV (Lx [p], Lx [p], pivot) ;
539 Int p, i, j, p2, phead, ptail, llen, ulen ;
543 for (p = 0 ; p < ulen ; p++)
547 PRINTF ((
"%d is pruned: %d. Lpend[j] %d Lip[j+1] %d\n",
548 j, Lpend [j] !=
EMPTY, Lpend [j], Lip [j+1])) ;
549 if (Lpend [j] ==
EMPTY)
553 for (p2 = 0 ; p2 < llen ; p2++)
555 if (pivrow == Li [p2])
559 PRINTF ((
"==== PRUNE: col j %d of L\n", j)) ;
562 for (p3 = 0 ; p3 < Llen [j] ; p3++)
564 PRINTF ((
"before: %i pivotal: %d\n", Li [p3],
565 Pinv [Li [p3]] >= 0)) ;
574 while (phead < ptail)
586 Li [phead] = Li [ptail] ;
589 Lx [phead] = Lx [ptail] ;
605 for (p3 = 0 ; p3 < Llen [j] ; p3++)
607 if (p3 == Lpend [j])
PRINTF ((
"----\n")) ;
608 PRINTF ((
"after: %i pivotal: %d\n", Li [p3],
609 Pinv [Li [p3]] >= 0)) ;
674 double abs_pivot, xsize, nunits,
tol,
memgrow ;
678 Int k, p, i, j, pivrow, kbar, diagrow, firstrow, lup, top,
scale, len ;
696 PRINTF ((
"input: lusize %d \n", lusize)) ;
707 for (k = 0 ; k <
n ; k++)
720 for (k = 0 ; k <
n ; k++)
723 Pinv [k] =
FLIP (k) ;
734 for (k = 0 ; k <
n ; k++)
736 PRINTF ((
"Initial P [%d] = %d\n", k,
P [k])) ;
744 for (k = 0 ; k <
n ; k++)
747 PRINTF ((
"\n\n==================================== k: %d\n", k)) ;
758 PRINTF ((
"lup %d lusize %g lup+nunits: %g\n", lup, (
double) lusize,
760 xsize = ((double) lup) + nunits ;
761 if (xsize > (
double) lusize)
764 xsize = (
memgrow * ((double) lusize) + 4*
n + 1) ;
767 PRINTF ((
"Matrix is too large (Int overflow)\n")) ;
771 newlusize = (size_t) (
memgrow * lusize + 2*
n + 1 ) ;
778 PRINTF ((
"Matrix is too large (LU)\n")) ;
782 PRINTF ((
"inc LU to %d done\n", lusize)) ;
796 for (i = 0 ; i <
n ; i++)
805 Lpend, Ap_pos, LU, lup, Llen, Lip, k1, PSinv) ;
808 PRINTF ((
"--- in U:\n")) ;
809 for (p = top ; p <
n ; p++)
811 PRINTF ((
"pattern of X for U: %d : %d pivot row: %d\n",
812 p, Stack [p], Pinv [Stack [p]])) ;
813 ASSERT (Flag [Stack [p]] == k) ;
815 PRINTF ((
"--- in L:\n")) ;
816 Li = (
Int *) (LU + Lip [k]);
817 for (p = 0 ; p < Llen [k] ; p++)
819 PRINTF ((
"pattern of X in L: %d : %d pivot row: %d\n",
820 p, Li [p], Pinv [Li [p]])) ;
821 ASSERT (Flag [Li [p]] == k) ;
824 for (i = 0 ; i <
n ; i++)
827 if (Flag [i] == k) p++ ;
836 k1, PSinv, Rs,
scale, Offp, Offi, Offx) ;
845 for (p = top ; p <
n ; p++)
847 PRINTF ((
"X for U %d : ", Stack [p])) ;
850 Li = (
Int *) (LU + Lip [k]) ;
851 for (p = 0 ; p < Llen [k] ; p++)
853 PRINTF ((
"X for L %d : ", Li [p])) ;
864 PRINTF ((
"k %d, diagrow = %d, UNFLIP (diagrow) = %d\n",
865 k, diagrow,
UNFLIP (diagrow))) ;
868 if (!
lpivot (diagrow, &pivrow, &pivot, &abs_pivot,
tol, X, LU, Lip,
869 Llen, k,
n, Pinv, &firstrow,
Common))
875 Common->numerical_rank = k+k1 ;
876 Common->singular_col = Q [k+k1] ;
878 if (
Common->halt_if_singular)
887 PRINTF ((
"\nk %d : Pivot row %d : ", k, pivrow)) ;
889 ASSERT (pivrow >= 0 && pivrow <
n) ;
890 ASSERT (Pinv [pivrow] < 0) ;
903 for (p = top, i = 0 ; p <
n ; p++, i++)
924 if (pivrow != diagrow)
928 PRINTF ((
">>>>>>>>>>>>>>>>> pivrow %d k %d off-diagonal\n",
930 if (Pinv [diagrow] < 0)
936 kbar =
FLIP (Pinv [pivrow]) ;
938 Pinv [diagrow] =
FLIP (kbar) ;
947 for (p = 0 ; p < len ; p++)
949 PRINTF ((
"Column %d of U: %d : ", k, Ui [p])) ;
953 for (p = 0 ; p < len ; p++)
955 PRINTF ((
"Column %d of L: %d : ", k, Li [p])) ;
964 prune (Lpend, Pinv, k, pivrow, LU, Uip, Lip, Ulen, Llen) ;
966 *lnz += Llen [k] + 1 ;
967 *unz += Ulen [k] + 1 ;
974 for (p = 0 ; p <
n ; p++)
976 Li = (
Int *) (LU + Lip [p]) ;
977 for (i = 0 ; i < Llen [p] ; i++)
979 Li [i] = Pinv [Li [i]] ;
984 for (i = 0 ; i <
n ; i++)
986 PRINTF ((
"P [%d] = %d Pinv [%d] = %d\n", i,
P [i], i, Pinv [i])) ;
988 for (i = 0 ; i <
n ; i++)
990 ASSERT (Pinv [i] >= 0 && Pinv [i] <
n) ;
1002 ASSERT ((
size_t) newlusize <= lusize) ;
1007 return (newlusize) ;
#define GET_POINTER(LU, Xip, Xlen, Xi, Xx, k, xlen)
static Int lsolve_symbolic(Int n, Int k, Int Ap [], Int Ai [], Int Q [], Int Pinv [], Int Stack [], Int Flag [], Int Lpend [], Int Ap_pos [], Unit LU [], Int lup, Int Llen [], Int Lip [], Int k1, Int PSinv [])
static void construct_column(Int k, Int Ap [], Int Ai [], Entry Ax [], Int Q [], Entry X [], Int k1, Int PSinv [], double Rs [], Int scale, Int Offp [], Int Offi [], Entry Offx [])
static void prune(Int Lpend [], Int Pinv [], Int k, Int pivrow, Unit *LU, Int Uip [], Int Lip [], Int Ulen [], Int Llen [])
static Int dfs(Int j, Int k, Int Pinv [], Int Llen [], Int Lip [], Int Stack [], Int Flag [], Int Lpend [], Int top, Unit LU [], Int *Lik, Int *plength, Int Ap_pos [])
static void lsolve_numeric(Int Pinv [], Unit *LU, Int Stack [], Int Lip [], Int top, Int n, Int Llen [], Entry X [])
static Int lpivot(Int diagrow, Int *p_pivrow, Entry *p_pivot, double *p_abs_pivot, double tol, Entry X [], Unit *LU, Int Lip [], Int Llen [], Int k, Int n, Int Pinv [], Int *p_firstrow, KLU_common *Common)
size_t KLU_kernel(Int n, Int Ap [], Int Ai [], Entry Ax [], Int Q [], size_t lusize, Int Pinv [], Int P [], Unit **p_LU, Entry Udiag [], Int Llen [], Int Ulen [], Int Lip [], Int Uip [], Int *lnz, Int *unz, Entry X [], Int Stack [], Int Flag [], Int Ap_pos [], Int Lpend [], Int k1, Int PSinv [], double Rs [], Int Offp [], Int Offi [], Entry Offx [], KLU_common *Common)
#define KLU_OUT_OF_MEMORY
#define ASSERT(expression)
#define MULT_SUB(c, a, b)