?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
[maven-release-plugin] prepare for next development iteration
[JKLU.git] / src / main / java / edu / ufl / cise / klu / tdouble / Dklu_factor.java
blob26b631271f58cf3edf200a152979bd9bf18fbc70
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 Iwork = Numeric.Iwork ; /* 5*maxblock for KLU_factor */
100 //Pblock = Iwork + 5*((int) Symbolic.maxblock) ; /* 1*maxblock for Pblock */
101 Pblock = new int [Symbolic.maxblock] ;
102 Common.nrealloc = 0 ;
103 scale = Common.scale ;
104 max_lnz_block = 1 ;
105 max_unz_block = 1 ;
106
107 /* compute the inverse of P from symbolic analysis. Will be updated to
108 * become the inverse of the numerical factorization when the factorization
109 * is done, for use in KLU_refactor */
110 if (!NDEBUG)
111 {
112 for (k = 0 ; k < n ; k++)
113 {
114 Pinv [k] = EMPTY ;
115 }
116 }
117 for (k = 0 ; k < n ; k++)
118 {
119 ASSERT (P [k] >= 0 && P [k] < n) ;
120 Pinv [P [k]] = k ;
121 }
122 if (!NDEBUG)
123 {
124 for (k = 0 ; k < n ; k++) ASSERT (Pinv [k] != EMPTY) ;
125 }
126
127 lnz = 0 ;
128 unz = 0 ;
129 Common.noffdiag = 0 ;
130 Offp [0] = 0 ;
131
132 /* ---------------------------------------------------------------------- */
133 /* optionally check input matrix and compute scale factors */
134 /* ---------------------------------------------------------------------- */
135
136 if (scale >= 0)
137 {
138 /* use Pnum as workspace. NOTE: scale factors are not yet permuted
139 * according to the final pivot row ordering, so Rs [oldrow] is the
140 * scale factor for A (oldrow,:), for the user's matrix A. Pnum is
141 * used as workspace in KLU_scale. When the factorization is done,
142 * the scale factors are permuted according to the final pivot row
143 * permutation, so that Rs [k] is the scale factor for the kth row of
144 * A(p,q) where p and q are the final row and column permutations. */
145 klu_scale (scale, n, Ap, Ai, Ax, Rs, Pnum, Common) ;
146 if (Common.status < KLU_OK)
147 {
148 /* matrix is invalid */
149 return ;
150 }
151 }
152
153 if (!NDEBUG)
154 {
155 if (scale > 0)
156 {
157 for (k = 0 ; k < n ; k++) PRINTF ("Rs [%d] %g\n", k, Rs [k]) ;
158 }
159 }
160
161 /* ---------------------------------------------------------------------- */
162 /* factor each block using klu */
163 /* ---------------------------------------------------------------------- */
164
165 for (block = 0 ; block < nblocks ; block++)
166 {
167
168 /* ------------------------------------------------------------------ */
169 /* the block is from rows/columns k1 to k2-1 */
170 /* ------------------------------------------------------------------ */
171
172 k1 = R [block] ;
173 k2 = R [block+1] ;
174 nk = k2 - k1 ;
175 PRINTF ("FACTOR BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1,k2-1,nk) ;
176
177 if (nk == 1)
178 {
179
180 /* -------------------------------------------------------------- */
181 /* singleton case */
182 /* -------------------------------------------------------------- */
183
184 poff = Offp [k1] ;
185 oldcol = Q [k1] ;
186 pend = Ap [oldcol+1] ;
187 //CLEAR (s) ;
188 s = 0.0;
189
190 if (scale <= 0)
191 {
192 /* no scaling */
193 for (p = Ap [oldcol] ; p < pend ; p++)
194 {
195 oldrow = Ai [p] ;
196 newrow = Pinv [oldrow] ;
197 if (newrow < k1)
198 {
199 Offi [poff] = oldrow ;
200 Offx [poff] = Ax [p] ;
201 poff++ ;
202 }
203 else
204 {
205 ASSERT (newrow == k1) ;
206 PRINTF ("singleton block %d", block) ;
207 PRINT_ENTRY (Ax [p]) ;
208 s = Ax [p] ;
209 }
210 }
211 }
212 else
213 {
214 /* row scaling. NOTE: scale factors are not yet permuted
215 * according to the pivot row permutation, so Rs [oldrow] is
216 * used below. When the factorization is done, the scale
217 * factors are permuted, so that Rs [newrow] will be used in
218 * klu_solve, klu_tsolve, and klu_rgrowth */
219 for (p = Ap [oldcol] ; p < pend ; p++)
220 {
221 oldrow = Ai [p] ;
222 newrow = Pinv [oldrow] ;
223 if (newrow < k1)
224 {
225 Offi [poff] = oldrow ;
226 //SCALE_DIV_ASSIGN (Offx [poff], Ax [p], Rs [oldrow]) ;
227 Offx [poff] = Ax [p] / Rs [oldrow] ;
228 poff++ ;
229 }
230 else
231 {
232 ASSERT (newrow == k1) ;
233 PRINTF ("singleton block %d ", block) ;
234 PRINT_ENTRY (Ax[p]) ;
235 s = Ax [p] / Rs [oldrow] ;
236 //SCALE_DIV_ASSIGN (s, Ax [p], Rs [oldrow]) ;
237 }
238 }
239 }
240
241 Udiag [k1] = s ;
242
243 if (IS_ZERO (s))
244 {
245 /* singular singleton */
246 Common.status = KLU_SINGULAR ;
247 Common.numerical_rank = k1 ;
248 Common.singular_col = oldcol ;
249 if (Common.halt_if_singular == 1)
250 {
251 return ;
252 }
253 }
254
255 Offp [k1+1] = poff ;
256 Pnum [k1] = P [k1] ;
257 lnz++ ;
258 unz++ ;
259
260 }
261 else
262 {
263
264 /* -------------------------------------------------------------- */
265 /* construct and factorize the kth block */
266 /* -------------------------------------------------------------- */
267
268 if (Lnz [block] < 0)
269 {
270 /* COLAMD was used - no estimate of fill-in */
271 /* use 10 times the nnz in A, plus n */
272 lsize = -(Common.initmem) ;
273 }
274 else
275 {
276 lsize = Common.initmem_amd * Lnz [block] + nk ;
277 }
278
279 /* allocates 1 arrays: LUbx [block] */
280 Numeric.LUsize [block] = klu_kernel_factor (
281 nk, Ap, Ai, Ax, Q,
282 lsize, LUbx, block, Udiag, k1, Llen, k1,
283 Ulen, k1, Lip, k1, Uip, k1, Pblock, lnz_block, unz_block,
284 X, Iwork, k1, Pinv, Rs, Offp, Offi, Offx, Common) ;
285
286 if (Common.status < KLU_OK ||
287 (Common.status == KLU_SINGULAR &&
288 Common.halt_if_singular == 1))
289 {
290 /* out of memory, invalid inputs, or singular */
291 return ;
292 }
293
294 PRINTF ("\n----------------------- L %d:\n", block) ;
295 if (!NDEBUG) ASSERT (klu_valid_LU (nk, TRUE, Lip, k1, Llen, k1, LUbx [block])) ;
296 PRINTF ("\n----------------------- U %d:\n", block) ;
297 if (!NDEBUG) ASSERT (klu_valid_LU (nk, FALSE, Uip, k1, Ulen, k1, LUbx [block])) ;
298
299 /* -------------------------------------------------------------- */
300 /* get statistics */
301 /* -------------------------------------------------------------- */
302
303 lnz += lnz_block[0] ;
304 unz += unz_block[0] ;
305 max_lnz_block = MAX (max_lnz_block, lnz_block[0]) ;
306 max_unz_block = MAX (max_unz_block, unz_block[0]) ;
307
308 if (Lnz [block] == EMPTY)
309 {
310 /* revise estimate for subsequent factorization */
311 Lnz [block] = MAX (lnz_block[0], unz_block[0]) ;
312 }
313
314 /* -------------------------------------------------------------- */
315 /* combine the klu row ordering with the symbolic pre-ordering */
316 /* -------------------------------------------------------------- */
317
318 PRINTF ("Pnum, 1-based:\n") ;
319 for (k = 0 ; k < nk ; k++)
320 {
321 ASSERT (k + k1 < n) ;
322 ASSERT (Pblock [k] + k1 < n) ;
323 Pnum [k + k1] = P [Pblock [k] + k1] ;
324 PRINTF ("Pnum (%d + %d + 1 = %d) = %d + 1 = %d\n",
325 k, k1, k+k1+1, Pnum [k+k1], Pnum [k+k1]+1) ;
326 }
327
328 /* the local pivot row permutation Pblock is no longer needed */
329 }
330 }
331 ASSERT (nzoff == Offp [n]) ;
332 PRINTF ("\n------------------- Off diagonal entries:\n") ;
333 if (!NDEBUG) ASSERT (klu_valid (n, Offp, Offi, Offx)) ;
334
335 Numeric.lnz = lnz ;
336 Numeric.unz = unz ;
337 Numeric.max_lnz_block = max_lnz_block ;
338 Numeric.max_unz_block = max_unz_block ;
339
340 /* compute the inverse of Pnum */
341 if (!NDEBUG)
342 {
343 for (k = 0 ; k < n ; k++)
344 {
345 Pinv [k] = EMPTY ;
346 }
347 }
348 for (k = 0 ; k < n ; k++)
349 {
350 ASSERT (Pnum [k] >= 0 && Pnum [k] < n) ;
351 Pinv [Pnum [k]] = k ;
352 }
353 if (!NDEBUG)
354 {
355 for (k = 0 ; k < n ; k++) ASSERT (Pinv [k] != EMPTY) ;
356 }
357
358 /* permute scale factors Rs according to pivotal row order */
359 if (scale > 0)
360 {
361 for (k = 0 ; k < n ; k++)
362 {
363 //REAL (X [k]) = Rs [Pnum [k]] ;
364 X [k] = Rs [Pnum [k]] ;
365 }
366 for (k = 0 ; k < n ; k++)
367 {
368 Rs [k] = X [k] ; //REAL (X [k]) ;
369 }
370 }
371
372 PRINTF ("\n------------------- Off diagonal entries, old:\n") ;
373 if (!NDEBUG) ASSERT (klu_valid (n, Offp, Offi, Offx)) ;
374
375 /* apply the pivot row permutations to the off-diagonal entries */
376 for (p = 0 ; p < nzoff ; p++)
377 {
378 ASSERT (Offi [p] >= 0 && Offi [p] < n) ;
379 Offi [p] = Pinv [Offi [p]] ;
380 }
381
382 PRINTF ("\n------------------- Off diagonal entries, new:\n") ;
383 if (!NDEBUG) ASSERT (klu_valid (n, Offp, Offi, Offx)) ;
384
385 if (!NDEBUG)
386 {
387 PRINTF ("\n ############# KLU_BTF_FACTOR done, nblocks %d\n",
388 nblocks);
389 double ss ;
390 Udiag = Numeric.Udiag ;
391 for (block = 0 ; block < nblocks && Common.status == KLU_OK ; block++)
392 {
393 k1 = R [block] ;
394 k2 = R [block+1] ;
395 nk = k2 - k1 ;
396 PRINTF ("\n======================KLU_factor output: k1 %d k2 %d nk %d\n",k1,k2,nk) ;
397 if (nk == 1)
398 {
399 PRINTF ("singleton ") ;
400 /* ENTRY_PRINT (singleton [block]) ; */
401 ss = Udiag [k1] ;
402 PRINT_ENTRY (ss) ;
403 }
404 else
405 {
406 double[] LU ;
407 Lip = Numeric.Lip ;
408 int Lip_offset = k1 ;
409 Llen = Numeric.Llen ;
410 int Llen_offset = k1 ;
411 LU = (double[]) Numeric.LUbx [block] ;
412 PRINTF ("\n---- L block %d\n", block);
413 if (!NDEBUG) ASSERT (klu_valid_LU (nk, TRUE, Lip, Lip_offset, Llen, Llen_offset, LU)) ;
414 Uip = Numeric.Uip ;
415 int Uip_offset = k1 ;
416 Ulen = Numeric.Ulen ;
417 int Ulen_offset = k1 ;
418 PRINTF ("\n---- U block %d\n", block) ;
419 if (!NDEBUG) ASSERT (klu_valid_LU (nk, FALSE, Uip, Uip_offset, Ulen, Ulen_offset, LU)) ;
420 }
421 }
422 }
423 }
424
425 /**
426 *
427 * @param Ap size n+1, column pointers
428 * @param Ai size nz, row indices
429 * @param Ax
430 * @param Symbolic
431 * @param Common
432 * @return null if error, or a valid KLU_numeric object if successful
433 */
434 public static KLU_numeric klu_factor(int[] Ap, int[] Ai, double[] Ax,
435 KLU_symbolic Symbolic, KLU_common Common)
436 {
437 int n, nzoff, nblocks, maxblock, k ;
438 int[] ok = new int [] {TRUE} ;
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 PRINTF ("KLU_factor: n %d nzoff %d nblocks %d maxblock %d\n",
466 n, nzoff, nblocks, maxblock) ;
467
468 /* ---------------------------------------------------------------------- */
469 /* get control parameters and make sure they are in the proper range */
470 /* ---------------------------------------------------------------------- */
471
472 Common.initmem_amd = MAX (1.0, Common.initmem_amd) ;
473 Common.initmem = MAX (1.0, Common.initmem) ;
474 Common.tol = MIN (Common.tol, 1.0) ;
475 Common.tol = MAX (0.0, Common.tol) ;
476 Common.memgrow = MAX (1.0, Common.memgrow) ;
477
478 /* ---------------------------------------------------------------------- */
479 /* allocate the Numeric object */
480 /* ---------------------------------------------------------------------- */
481
482 /* this will not cause int overflow (already checked by KLU_symbolic) */
483 n1 = n + 1 ;
484 nzoff1 = nzoff + 1 ;
485
486 //Numeric = klu_malloc (sizeof (KLU_numeric), 1, Common) ;
487 try
488 {
489 Numeric = new KLU_numeric();
490 }
491 catch (OutOfMemoryError e)
492 {
493 Common.status = KLU_OUT_OF_MEMORY ;
494 return (null) ;
495 }
496 Numeric.n = n ;
497 Numeric.nblocks = nblocks ;
498 Numeric.nzoff = nzoff ;
499 Numeric.Pnum = klu_malloc_int (n, Common) ;
500 Numeric.Offp = klu_malloc_int (n1, Common) ;
501 Numeric.Offi = klu_malloc_int (nzoff1, Common) ;
502 Numeric.Offx = klu_malloc_dbl (nzoff1, Common) ;
503
504 Numeric.Lip = klu_malloc_int (n, Common) ;
505 Numeric.Uip = klu_malloc_int (n, Common) ;
506 Numeric.Llen = klu_malloc_int (n, Common) ;
507 Numeric.Ulen = klu_malloc_int (n, Common) ;
508
509 Numeric.LUsize = klu_malloc_int (nblocks, Common) ;
510
511 //Numeric.LUbx = klu_malloc (nblocks, sizeof (double[]), Common) ;
512 Numeric.LUbx = new double [nblocks][] ;
513 if (Numeric.LUbx != null)
514 {
515 for (k = 0 ; k < nblocks ; k++)
516 {
517 Numeric.LUbx [k] = null ;
518 }
519 }
520
521 Numeric.Udiag = klu_malloc_dbl (n, Common) ;
522
523 if (Common.scale > 0)
524 {
525 Numeric.Rs = klu_malloc_dbl (n, Common) ;
526 }
527 else
528 {
529 /* no scaling */
530 Numeric.Rs = null ;
531 }
532
533 Numeric.Pinv = klu_malloc_int (n, Common) ;
534
535 /* allocate permanent workspace for factorization and solve. Note that the
536 * solver will use an Xwork of size 4n, whereas the factorization codes use
537 * an Xwork of size n and integer space (Iwork) of size 6n. KLU_condest
538 * uses an Xwork of size 2n. Total size is:
539 *
540 * n*sizeof(double) + max (6*maxblock*sizeof(Int), 3*n*sizeof(double))
541 */
542 //s = klu_mult_size_t (n, sizeof (double), ok) ;
543 s = n ;
544 //n3 = klu_mult_size_t (n, 3 * sizeof (double), ok) ;
545 n3 = 3 * n ;
546 //b6 = klu_mult_size_t (maxblock, 6 * sizeof (Int), ok) ;
547 b6 = 6 * maxblock ;
548 Numeric.worksize = klu_add_size_t (s, MAX (n3, b6), ok) ;
549 try
550 {
551 if (ok[0] == 0) throw new OutOfMemoryError() ;
552
553 //Numeric.Work = klu_malloc (Numeric.worksize, 1, Common) ;
554 Numeric.Work = new double [Numeric.worksize] ;
555 Numeric.Xwork = Numeric.Work ;
556 //Numeric.Iwork = (Int[]) ((double[]) Numeric.Xwork + n) ;
557 Numeric.Iwork = new int [b6] ;
558 }
559 catch (OutOfMemoryError e)
560 {
561 /* out of memory or problem too large */
562 Common.status = ok[0] == 1 ? KLU_OUT_OF_MEMORY : KLU_TOO_LARGE ;
563 //klu_free_numeric (Numeric, Common) ;
564 Numeric = null;
565 return (null) ;
566 }
567
568 /* ---------------------------------------------------------------------- */
569 /* factorize the blocks */
570 /* ---------------------------------------------------------------------- */
571
572 factor2 (Ap, Ai, Ax, Symbolic, Numeric, Common) ;
573
574 /* ---------------------------------------------------------------------- */
575 /* return or free the Numeric object */
576 /* ---------------------------------------------------------------------- */
577
578 if (Common.status < KLU_OK)
579 {
580 /* out of memory or inputs invalid */
581 //klu_free_numeric (Numeric, Common) ;
582 Numeric = null ;
583 }
584 else if (Common.status == KLU_SINGULAR)
585 {
586 if (Common.halt_if_singular == 1)
587 {
588 /* Matrix is singular, and the Numeric object is only partially
589 * defined because we halted early. This is the default case for
590 * a singular matrix. */
591 //klu_free_numeric (Numeric, Common) ;
592 Numeric = null ;
593
594 }
595 }
596 else if (Common.status == KLU_OK)
597 {
598 /* successful non-singular factorization */
599 Common.numerical_rank = n ;
600 Common.singular_col = n ;
601 }
602 return (Numeric) ;
603 }
604
605 }
A sparse LU factorization algorithm in Java