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