?  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_refactor.java
blobcc34ae9e04824f59d9a70090761074de246512b1
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_memory.klu_malloc_dbl;
32 import static edu.ufl.cise.klu.tdouble.Dklu_scale.klu_scale;
33 import static edu.ufl.cise.klu.tdouble.Dklu_dump.klu_valid;
34 import static edu.ufl.cise.klu.tdouble.Dklu_dump.klu_valid_LU;
35
36 /**
37 * Factor the matrix, after ordering and analyzing it with KLU_analyze, and
38 * factoring it once with KLU_factor. This routine cannot do any numerical
39 * pivoting. The pattern of the input matrix (Ap, Ai) must be identical to
40 * the pattern given to KLU_factor.
41 */
42 public class Dklu_refactor extends Dklu_internal {
43
44 /**
45 *
46 * @param Ap size n+1, column pointers
47 * @param Ai size nz, row indices
48 * @param Ax
49 * @param Symbolic
50 * @param Numeric
51 * @param Common
52 * @return true if successful, false otherwise
53 */
54 public static int klu_refactor(int[] Ap, int[] Ai, double[] Ax,
55 KLU_symbolic Symbolic, KLU_numeric Numeric, KLU_common Common)
56 {
57 double ukk, ujk, s ;
58 double[] Offx, Lx, Ux, X, Az, Udiag ;
59 double[] Rs ;
60 int[] Q, R, Pnum, Offp, Offi, Pinv, Lip, Uip, Llen, Ulen ;
61 /*int[]*/double[] Ui, Li ;
62 double[][] LUbx ;
63 double[] LU ;
64 int k1, k2, nk, k, block, oldcol, pend, oldrow, n, p, newrow, scale,
65 nblocks, poff, i, j, up, maxblock, nzoff ;
66 int[] ulen = new int[1] ;
67 int[] Ui_offset = new int[1] ;
68 int[] Ux_offset = new int[1] ;
69 int[] llen = new int[1] ;
70 int[] Li_offset = new int[1] ;
71 int[] Lx_offset = new int[1] ;
72
73 /* ---------------------------------------------------------------------- */
74 /* check inputs */
75 /* ---------------------------------------------------------------------- */
76
77 if (Common == null)
78 {
79 return (FALSE) ;
80 }
81 Common.status = KLU_OK ;
82
83 if (Numeric == null)
84 {
85 /* invalid Numeric object */
86 Common.status = KLU_INVALID ;
87 return (FALSE) ;
88 }
89
90 Common.numerical_rank = EMPTY ;
91 Common.singular_col = EMPTY ;
92
93 Az = (double[]) Ax ;
94
95 /* ---------------------------------------------------------------------- */
96 /* get the contents of the Symbolic object */
97 /* ---------------------------------------------------------------------- */
98
99 n = Symbolic.n ;
100 Q = Symbolic.Q ;
101 R = Symbolic.R ;
102 nblocks = Symbolic.nblocks ;
103 maxblock = Symbolic.maxblock ;
104
105 /* ---------------------------------------------------------------------- */
106 /* get the contents of the Numeric object */
107 /* ---------------------------------------------------------------------- */
108
109 Pnum = Numeric.Pnum ;
110 Offp = Numeric.Offp ;
111 Offi = Numeric.Offi ;
112 Offx = Numeric.Offx ;
113
114 LUbx = Numeric.LUbx ;
115
116 scale = Common.scale ;
117 if (scale > 0)
118 {
119 /* factorization was not scaled, but refactorization is scaled */
120 if (Numeric.Rs == null)
121 {
122 Numeric.Rs = klu_malloc_dbl (n, Common) ;
123 if (Common.status < KLU_OK)
124 {
125 Common.status = KLU_OUT_OF_MEMORY ;
126 return (FALSE) ;
127 }
128 }
129 }
130 else
131 {
132 /* no scaling for refactorization; ensure Numeric.Rs is freed. This
133 * does nothing if Numeric.Rs is already null. */
134 //Numeric.Rs = KLU_free (Numeric.Rs, n, sizeof (double), Common) ;
135 Numeric.Rs = null ;
136 }
137 Rs = Numeric.Rs ;
138
139 Pinv = Numeric.Pinv ;
140 X = Numeric.Xwork ;
141 Common.nrealloc = 0 ;
142 Udiag = Numeric.Udiag ;
143 nzoff = Symbolic.nzoff ;
144
145 /* ---------------------------------------------------------------------- */
146 /* check the input matrix compute the row scale factors, Rs */
147 /* ---------------------------------------------------------------------- */
148
149 /* do no scale, or check the input matrix, if scale < 0 */
150 if (scale >= 0)
151 {
152 /* check for out-of-range indices, but do not check for duplicates */
153 if (klu_scale (scale, n, Ap, Ai, Ax, Rs, null, Common) == 0)
154 {
155 return (FALSE) ;
156 }
157 }
158
159 /* ---------------------------------------------------------------------- */
160 /* clear workspace X */
161 /* ---------------------------------------------------------------------- */
162
163 for (k = 0 ; k < maxblock ; k++)
164 {
165 /* X [k] = 0 ; */
166 CLEAR (X, k) ;
167 }
168
169 poff = 0 ;
170
171 /* ---------------------------------------------------------------------- */
172 /* factor each block */
173 /* ---------------------------------------------------------------------- */
174
175 if (scale <= 0)
176 {
177
178 /* ------------------------------------------------------------------ */
179 /* no scaling */
180 /* ------------------------------------------------------------------ */
181
182 for (block = 0 ; block < nblocks ; block++)
183 {
184
185 /* -------------------------------------------------------------- */
186 /* the block is from rows/columns k1 to k2-1 */
187 /* -------------------------------------------------------------- */
188
189 k1 = R [block] ;
190 k2 = R [block+1] ;
191 nk = k2 - k1 ;
192
193 if (nk == 1)
194 {
195
196 /* ---------------------------------------------------------- */
197 /* singleton case */
198 /* ---------------------------------------------------------- */
199
200 oldcol = Q [k1] ;
201 pend = Ap [oldcol+1] ;
202 s = 0 ; //CLEAR (s) ;
203 for (p = Ap [oldcol] ; p < pend ; p++)
204 {
205 newrow = Pinv [Ai [p]] - k1 ;
206 if (newrow < 0 && poff < nzoff)
207 {
208 /* entry in off-diagonal block */
209 Offx [poff] = Az [p] ;
210 poff++ ;
211 }
212 else
213 {
214 /* singleton */
215 s = Az [p] ;
216 }
217 }
218 Udiag [k1] = s ;
219
220 }
221 else
222 {
223
224 /* ---------------------------------------------------------- */
225 /* construct and factor the kth block */
226 /* ---------------------------------------------------------- */
227
228 Lip = Numeric.Lip ;
229 int Lip_offset = k1 ;
230 Llen = Numeric.Llen ;
231 int Llen_offset = k1 ;
232 Uip = Numeric.Uip ;
233 int Uip_offset = k1 ;
234 Ulen = Numeric.Ulen ;
235 int Ulen_offset = k1 ;
236 LU = LUbx [block] ;
237
238 for (k = 0 ; k < nk ; k++)
239 {
240
241 /* ------------------------------------------------------ */
242 /* scatter kth column of the block into workspace X */
243 /* ------------------------------------------------------ */
244
245 oldcol = Q [k+k1] ;
246 pend = Ap [oldcol+1] ;
247 for (p = Ap [oldcol] ; p < pend ; p++)
248 {
249 newrow = Pinv [Ai [p]] - k1 ;
250 if (newrow < 0 && poff < nzoff)
251 {
252 /* entry in off-diagonal block */
253 Offx [poff] = Az [p] ;
254 poff++ ;
255 }
256 else
257 {
258 /* (newrow,k) is an entry in the block */
259 X [newrow] = Az [p] ;
260 }
261 }
262
263 /* ------------------------------------------------------ */
264 /* compute kth column of U, and update kth column of A */
265 /* ------------------------------------------------------ */
266
267 Ui = Ux = GET_POINTER (LU, Uip, Uip_offset, Ulen, Ulen_offset,
268 Ui_offset, Ux_offset, k, ulen) ;
269 for (up = 0 ; up < ulen[0] ; up++)
270 {
271 j = (int) Ui [Ui_offset[0] + up] ;
272 ujk = X [j] ;
273 /* X [j] = 0 ; */
274 CLEAR (X, j) ;
275 Ux [Ux_offset[0] + up] = ujk ;
276 Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset,
277 Li_offset, Lx_offset, j, llen) ;
278 for (p = 0 ; p < llen[0] ; p++)
279 {
280 //MULT_SUB (X [Li [p]], Lx [p], ujk) ;
281 X [(int) Li [Li_offset[0] + p]] -= Lx [Lx_offset[0] + p] * ujk ;
282 }
283 }
284 /* get the diagonal entry of U */
285 ukk = X [k] ;
286 /* X [k] = 0 ; */
287 CLEAR (X, k) ;
288 /* singular case */
289 if (IS_ZERO (ukk))
290 {
291 /* matrix is numerically singular */
292 Common.status = KLU_SINGULAR ;
293 if (Common.numerical_rank == EMPTY)
294 {
295 Common.numerical_rank = k+k1 ;
296 Common.singular_col = Q [k+k1] ;
297 }
298 if (Common.halt_if_singular != 0)
299 {
300 /* do not continue the factorization */
301 return (FALSE) ;
302 }
303 }
304 Udiag [k+k1] = ukk ;
305 /* gather and divide by pivot to get kth column of L */
306 Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset,
307 Li_offset, Lx_offset, k, llen) ;
308 for (p = 0 ; p < llen[0] ; p++)
309 {
310 i = (int) Li [Li_offset[0] + p] ;
311 //DIV (Lx [p], X [i], ukk) ;
312 Lx [Lx_offset[0] + p] = X [i] / ukk ;
313 CLEAR (X, i) ;
314 }
315
316 }
317 }
318 }
319
320 }
321 else
322 {
323
324 /* ------------------------------------------------------------------ */
325 /* scaling */
326 /* ------------------------------------------------------------------ */
327
328 for (block = 0 ; block < nblocks ; block++)
329 {
330
331 /* -------------------------------------------------------------- */
332 /* the block is from rows/columns k1 to k2-1 */
333 /* -------------------------------------------------------------- */
334
335 k1 = R [block] ;
336 k2 = R [block+1] ;
337 nk = k2 - k1 ;
338
339 if (nk == 1)
340 {
341
342 /* ---------------------------------------------------------- */
343 /* singleton case */
344 /* ---------------------------------------------------------- */
345
346 oldcol = Q [k1] ;
347 pend = Ap [oldcol+1] ;
348 s = 0 ; //CLEAR (s) ;
349 for (p = Ap [oldcol] ; p < pend ; p++)
350 {
351 oldrow = Ai [p] ;
352 newrow = Pinv [oldrow] - k1 ;
353 if (newrow < 0 && poff < nzoff)
354 {
355 /* entry in off-diagonal block */
356 Offx [poff] = Az [p] / Rs [oldrow] ;
357 //SCALE_DIV_ASSIGN (Offx [poff], Az [p], Rs [oldrow]) ;
358 poff++ ;
359 }
360 else
361 {
362 /* singleton */
363 s = Az [p] / Rs [oldrow] ;
364 //SCALE_DIV_ASSIGN (s, Az [p], Rs [oldrow]) ;
365 }
366 }
367 Udiag [k1] = s ;
368
369 }
370 else
371 {
372
373 /* ---------------------------------------------------------- */
374 /* construct and factor the kth block */
375 /* ---------------------------------------------------------- */
376
377 Lip = Numeric.Lip ;
378 int Lip_offset = k1 ;
379 Llen = Numeric.Llen ;
380 int Llen_offset = k1 ;
381 Uip = Numeric.Uip ;
382 int Uip_offset = k1 ;
383 Ulen = Numeric.Ulen ;
384 int Ulen_offset = k1 ;
385 LU = LUbx [block] ;
386
387 for (k = 0 ; k < nk ; k++)
388 {
389
390 /* ------------------------------------------------------ */
391 /* scatter kth column of the block into workspace X */
392 /* ------------------------------------------------------ */
393
394 oldcol = Q [k+k1] ;
395 pend = Ap [oldcol+1] ;
396 for (p = Ap [oldcol] ; p < pend ; p++)
397 {
398 oldrow = Ai [p] ;
399 newrow = Pinv [oldrow] - k1 ;
400 if (newrow < 0 && poff < nzoff)
401 {
402 /* entry in off-diagonal part */
403 //SCALE_DIV_ASSIGN (Offx [poff], Az [p], Rs [oldrow]);
404 Offx [poff] = Az [p] / Rs [oldrow] ;
405 poff++ ;
406 }
407 else
408 {
409 /* (newrow,k) is an entry in the block */
410 //SCALE_DIV_ASSIGN (X [newrow], Az [p], Rs [oldrow]) ;
411 X [newrow] = Az [p] / Rs [oldrow] ;
412 }
413 }
414
415 /* ------------------------------------------------------ */
416 /* compute kth column of U, and update kth column of A */
417 /* ------------------------------------------------------ */
418
419 Ui = Ux = GET_POINTER (LU, Uip, Uip_offset, Ulen, Ulen_offset,
420 Ui_offset, Ux_offset, k, ulen) ;
421 for (up = 0 ; up < ulen[0] ; up++)
422 {
423 j = (int) Ui [Ui_offset[0] + up] ;
424 ujk = X [j] ;
425 /* X [j] = 0 ; */
426 CLEAR (X, j) ;
427 Ux [Ux_offset[0] + up] = ujk ;
428 Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset,
429 Li_offset, Lx_offset, j, llen) ;
430 for (p = 0 ; p < llen[0] ; p++)
431 {
432 //MULT_SUB (X [Li [p]], Lx [p], ujk) ;
433 X [(int) Li [Li_offset[0] + p]] -= Lx [Lx_offset[0] + p] * ujk ;
434 }
435 }
436 /* get the diagonal entry of U */
437 ukk = X [k] ;
438 /* X [k] = 0 ; */
439 CLEAR (X, k) ;
440 /* singular case */
441 if (IS_ZERO (ukk))
442 {
443 /* matrix is numerically singular */
444 Common.status = KLU_SINGULAR ;
445 if (Common.numerical_rank == EMPTY)
446 {
447 Common.numerical_rank = k+k1 ;
448 Common.singular_col = Q [k+k1] ;
449 }
450 if (Common.halt_if_singular != 0)
451 {
452 /* do not continue the factorization */
453 return (FALSE) ;
454 }
455 }
456 Udiag [k+k1] = ukk ;
457 /* gather and divide by pivot to get kth column of L */
458 Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset,
459 Li_offset, Lx_offset, k, llen) ;
460 for (p = 0 ; p < llen[0] ; p++)
461 {
462 i = (int) Li [Li_offset[0] + p] ;
463 //DIV (Lx [p], X [i], ukk) ;
464 Lx [Lx_offset[0] + p] = X [i] / ukk ;
465 CLEAR (X, i) ;
466 }
467 }
468 }
469 }
470 }
471
472 /* ---------------------------------------------------------------------- */
473 /* permute scale factors Rs according to pivotal row order */
474 /* ---------------------------------------------------------------------- */
475
476 if (scale > 0)
477 {
478 for (k = 0 ; k < n ; k++)
479 {
480 X [k] = Rs [Pnum [k]] ;
481 //REAL (X [k]) = Rs [Pnum [k]] ;
482 }
483 for (k = 0 ; k < n ; k++)
484 {
485 Rs [k] = X [k] ;
486 //Rs [k] = REAL (X [k]) ;
487 }
488 }
489
490 if (!NDEBUG)
491 {
492 ASSERT (Offp [n] == poff) ;
493 ASSERT (Symbolic.nzoff == poff) ;
494 PRINTF (("\n------------------- Off diagonal entries, new:\n")) ;
495 if (!NDEBUG) ASSERT (klu_valid (n, Offp, Offi, Offx)) ;
496 if (Common.status == KLU_OK)
497 {
498 PRINTF ("\n ########### KLU_BTF_REFACTOR done, nblocks %d\n",
499 nblocks);
500 for (block = 0 ; block < nblocks ; block++)
501 {
502 k1 = R [block] ;
503 k2 = R [block+1] ;
504 nk = k2 - k1 ;
505 PRINTF (
506 "\n================KLU_refactor output: k1 %d k2 %d nk %d\n",
507 k1, k2, nk) ;
508 if (nk == 1)
509 {
510 PRINTF ("singleton ") ;
511 PRINT_ENTRY (Udiag [k1]) ;
512 }
513 else
514 {
515 Lip = Numeric.Lip ;
516 int Lip_offset = k1 ;
517 Llen = Numeric.Llen ;
518 int Llen_offset = k1 ;
519 LU = (double[]) Numeric.LUbx [block] ;
520 PRINTF ("\n---- L block %d\n", block) ;
521 if (!NDEBUG) ASSERT (klu_valid_LU (nk, TRUE, Lip, Lip_offset, Llen, Llen_offset, LU)) ;
522 Uip = Numeric.Uip ;
523 int Uip_offset = k1 ;
524 Ulen = Numeric.Ulen ;
525 int Ulen_offset = k1 ;
526 PRINTF ("\n---- U block %d\n", block) ;
527 if (!NDEBUG) ASSERT (klu_valid_LU (nk, FALSE, Uip, Uip_offset, Ulen, Ulen_offset, LU)) ;
528 }
529 }
530 }
531 }
532
533 return (TRUE) ;
534 }
535
536 }
A sparse LU factorization algorithm in Java