?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixing pointers up to Dklu_solve.
[JKLU.git] / src / main / java / edu / ufl / cise / klu / tdouble / Dklu_diagnostics.java
blobda6d522e9cb1d5c0509bd0d8ab53c72c9a3bd5d9
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_tsolve.klu_tsolve;
32 import static edu.ufl.cise.klu.tdouble.Dklu_solve.klu_solve;
33
34 /**
35 * Linear algebraic diagnostics.
36 */
37 public class Dklu_diagnostics extends Dklu_internal
38 {
39
40 /**
41 * Compute the reciprocal pivot growth factor. In MATLAB notation:
42 *
43 * rgrowth = min (max (abs ((R \ A (p,q)) - F))) ./ max (abs (U)))
44 *
45 * Takes O(|A|+|U|) time.
46 *
47 * @param Ap
48 * @param Ai
49 * @param Ax
50 * @param Symbolic
51 * @param Numeric
52 * @param Common
53 * @return TRUE if successful, FALSE otherwise
54 */
55 public static int klu_rgrowth(int[] Ap, int[] Ai, double[] Ax,
56 KLU_symbolic Symbolic, KLU_numeric Numeric, KLU_common Common)
57 {
58 double temp, max_ai, max_ui, min_block_rgrowth ;
59 double aik ;
60 int[] Q, Pinv ;
61 int[] Ulen, Uip ;
62 double[] LU ;
63 double[] Aentry, Ux, Ukk ;
64 double[] Rs ;
65 int i, newrow, oldrow, k1, k2, nk, j, oldcol, k, pend ;
66 int Uip_offset = 0, Ulen_offset = 0, Ukk_offset = 0 ;
67 int[] len = new int[1] ;
68 int[] Ui_offset = new int[1] ;
69 int[] Ux_offset = new int[1] ;
70
71 /* ---------------------------------------------------------------------- */
72 /* check inputs */
73 /* ---------------------------------------------------------------------- */
74
75 if (Common == null)
76 {
77 return (FALSE) ;
78 }
79
80 if (Symbolic == null || Ap == null || Ai == null || Ax == null)
81 {
82 Common.status = KLU_INVALID ;
83 return (FALSE) ;
84 }
85
86 if (Numeric == null)
87 {
88 /* treat this as a singular matrix */
89 Common.rgrowth = 0 ;
90 Common.status = KLU_SINGULAR ;
91 return (TRUE) ;
92 }
93 Common.status = KLU_OK ;
94
95 /* ---------------------------------------------------------------------- */
96 /* compute the reciprocal pivot growth */
97 /* ---------------------------------------------------------------------- */
98
99 Aentry = (double[]) Ax ;
100 Pinv = Numeric.Pinv ;
101 Rs = Numeric.Rs ;
102 Q = Symbolic.Q ;
103 Common.rgrowth = 1 ;
104
105 for (i = 0 ; i < Symbolic.nblocks ; i++)
106 {
107 k1 = Symbolic.R[i] ;
108 k2 = Symbolic.R[i+1] ;
109 nk = k2 - k1 ;
110 if (nk == 1)
111 {
112 continue ; /* skip singleton blocks */
113 }
114 LU = Numeric.LUbx[i] ;
115 Uip = Numeric.Uip ;
116 Uip_offset += k1 ;
117 Ulen = Numeric.Ulen ;
118 Ulen_offset += k1 ;
119 Ukk = Numeric.Udiag ;
120 Ukk_offset += k1 ;
121 min_block_rgrowth = 1 ;
122 for (j = 0 ; j < nk ; j++)
123 {
124 max_ai = 0 ;
125 max_ui = 0 ;
126 oldcol = Q[j + k1] ;
127 pend = Ap [oldcol + 1] ;
128 for (k = Ap [oldcol] ; k < pend ; k++)
129 {
130 oldrow = Ai [k] ;
131 newrow = Pinv [oldrow] ;
132 if (newrow < k1)
133 {
134 continue ; /* skip entry outside the block */
135 }
136 ASSERT (newrow < k2) ;
137 if (Rs != null)
138 {
139 //SCALE_DIV_ASSIGN (aik, Aentry [k], Rs [newrow]) ;
140 aik = Aentry [k] / Rs [newrow] ;
141 }
142 else
143 {
144 aik = Aentry [k] ;
145 }
146 temp = ABS ( aik ) ;
147 //ABS (temp, aik) ;
148 if (temp > max_ai)
149 {
150 max_ai = temp ;
151 }
152 }
153
154 Ux = GET_POINTER (LU, Uip, Uip_offset, Ulen, Ulen_offset,
155 Ui_offset, Ux_offset, j, len) ;
156 for (k = 0 ; k < len[0] ; k++)
157 {
158 temp = ABS (Ux [Ux_offset[0] + k]) ;
159 //ABS (temp, Ux [k]) ;
160 if (temp > max_ui)
161 {
162 max_ui = temp ;
163 }
164 }
165 /* consider the diagonal element */
166 temp = ABS (Ukk [Ukk_offset + j]) ;
167 //ABS (temp, Ukk [j]) ;
168 if (temp > max_ui)
169 {
170 max_ui = temp ;
171 }
172
173 /* if max_ui is 0, skip the column */
174 if (SCALAR_IS_ZERO (max_ui))
175 {
176 continue ;
177 }
178 temp = max_ai / max_ui ;
179 if (temp < min_block_rgrowth)
180 {
181 min_block_rgrowth = temp ;
182 }
183 }
184
185 if (min_block_rgrowth < Common.rgrowth)
186 {
187 Common.rgrowth = min_block_rgrowth ;
188 }
189 }
190 return (TRUE) ;
191 }
192
193 /**
194 * Estimate the condition number. Uses Higham and Tisseur's algorithm
195 * (A block algorithm for matrix 1-norm estimation, with applications to
196 * 1-norm pseudospectra, SIAM J. Matrix Anal. Appl., 21(4):1185-1201, 2000.
197 *
198 * Takes about O(|A|+5*(|L|+|U|)) time
199 *
200 * @param Ap
201 * @param Ax
202 * @param Symbolic
203 * @param Numeric
204 * @param Common
205 * @return TRUE if successful, FALSE otherwise
206 */
207 public static int klu_condest(int[] Ap, double[] Ax, KLU_symbolic Symbolic,
208 KLU_numeric Numeric, KLU_common Common)
209 {
210 double xj, Xmax, csum, anorm, ainv_norm, est_old, est_new, abs_value ;
211 double[] Udiag, Aentry, X, S ;
212 int i, j, jmax, jnew, pend, n ;
213 int unchanged ;
214
215 /* ---------------------------------------------------------------------- */
216 /* check inputs */
217 /* ---------------------------------------------------------------------- */
218
219 if (Common == null)
220 {
221 return (FALSE) ;
222 }
223 if (Symbolic == null || Ap == null || Ax == null)
224 {
225 Common.status = KLU_INVALID ;
226 return (FALSE) ;
227 }
228 abs_value = 0 ;
229 if (Numeric == null)
230 {
231 /* treat this as a singular matrix */
232 Common.condest = 1 / abs_value ;
233 Common.status = KLU_SINGULAR ;
234 return (TRUE) ;
235 }
236 Common.status = KLU_OK ;
237
238 /* ---------------------------------------------------------------------- */
239 /* get inputs */
240 /* ---------------------------------------------------------------------- */
241
242 n = Symbolic.n ;
243 Udiag = Numeric.Udiag ;
244
245 /* ---------------------------------------------------------------------- */
246 /* check if diagonal of U has a zero on it */
247 /* ---------------------------------------------------------------------- */
248
249 for (i = 0 ; i < n ; i++)
250 {
251 //ABS (abs_value, Udiag [i]) ;
252 abs_value = ABS (Udiag [i]) ;
253 if (SCALAR_IS_ZERO (abs_value))
254 {
255 Common.condest = 1 / abs_value ;
256 Common.status = KLU_SINGULAR ;
257 return (TRUE) ;
258 }
259 }
260
261 /* ---------------------------------------------------------------------- */
262 /* compute 1-norm (maximum column sum) of the matrix */
263 /* ---------------------------------------------------------------------- */
264
265 anorm = 0.0 ;
266 Aentry = (double[]) Ax ;
267 for (i = 0 ; i < n ; i++)
268 {
269 pend = Ap [i + 1] ;
270 csum = 0.0 ;
271 for (j = Ap [i] ; j < pend ; j++)
272 {
273 //ABS (abs_value, Aentry [j]) ;
274 abs_value = ABS (Aentry [j]) ;
275 csum += abs_value ;
276 }
277 if (csum > anorm)
278 {
279 anorm = csum ;
280 }
281 }
282
283 /* ---------------------------------------------------------------------- */
284 /* compute estimate of 1-norm of inv (A) */
285 /* ---------------------------------------------------------------------- */
286
287 /* get workspace (size 2*n double's) */
288 X = Numeric.Xwork ; /* size n space used in KLU_solve, tsolve */
289 //X += n ; /* X is size n */
290 int X_offset = n ;
291 //S = X + n ; /* S is size n */
292 S = X ;
293 int S_offset = 2*n ;
294
295 for (i = 0 ; i < n ; i++)
296 {
297 CLEAR (S, S_offset + i) ;
298 CLEAR (X, X_offset + i) ;
299 //REAL (X [i]) = 1.0 / ((double) n) ;
300 X [X_offset + i] = 1.0 / ((double) n);
301 }
302 jmax = 0 ;
303
304 ainv_norm = 0.0 ;
305 for (i = 0 ; i < 5 ; i++)
306 {
307 if (i > 0)
308 {
309 /* X [jmax] is the largest entry in X */
310 for (j = 0 ; j < n ; j++)
311 {
312 CLEAR (X, X_offset + j) ;
313 }
314 //REAL (X [jmax]) = 1 ;
315 X [X_offset + jmax] = 1 ;
316 }
317
318 klu_solve (Symbolic, Numeric, n, 1, (double[]) X, X_offset, Common) ;
319 est_old = ainv_norm ;
320 ainv_norm = 0.0 ;
321
322 for (j = 0 ; j < n ; j++)
323 {
324 /* ainv_norm += ABS (X [j]) ;*/
325 //ABS (abs_value, X [j]) ;
326 abs_value = ABS (X [X_offset + j]) ;
327 ainv_norm += abs_value ;
328 }
329
330 unchanged = TRUE ;
331
332 for (j = 0 ; j < n ; j++)
333 {
334 double s = (X [X_offset + j] >= 0) ? 1 : -1 ;
335 if (s != S [S_offset + j]) // s != REAL (S [j])
336 {
337 S [S_offset + j] = s ;
338 unchanged = FALSE ;
339 }
340 }
341
342 if (i > 0 && (ainv_norm <= est_old || unchanged == 1))
343 {
344 break ;
345 }
346
347 for (j = 0 ; j < n ; j++)
348 {
349 X [j] = S [S_offset + j] ;
350 }
351
352 /* do a transpose solve */
353 klu_tsolve (Symbolic, Numeric, n, 1, X, X_offset, Common) ;
354
355 /* jnew = the position of the largest entry in X */
356 jnew = 0 ;
357 Xmax = 0 ;
358 for (j = 0 ; j < n ; j++)
359 {
360 //ABS (xj, X [j]) ;
361 xj = ABS (X [X_offset + j]) ;
362 if (xj > Xmax)
363 {
364 Xmax = xj ;
365 jnew = j ;
366 }
367 }
368 if (i > 0 && jnew == jmax)
369 {
370 /* the position of the largest entry did not change
371 * from the previous iteration */
372 break ;
373 }
374 jmax = jnew ;
375 }
376
377 /* ---------------------------------------------------------------------- */
378 /* compute another estimate of norm(inv(A),1), and take the largest one */
379 /* ---------------------------------------------------------------------- */
380
381 for (j = 0 ; j < n ; j++)
382 {
383 CLEAR (X, X_offset + j) ;
384 if (j % 2 != 0)
385 {
386 //REAL (X [j]) = 1 + ((double) j) / ((double) (n-1)) ;
387 X [X_offset + j] = 1 + ((double) j) / ((double) (n-1)) ;
388 }
389 else
390 {
391 //REAL (X [j]) = -1 - ((double) j) / ((double) (n-1)) ;
392 X [X_offset + j] = -1 - ((double) j) / ((double) (n-1)) ;
393 }
394 }
395
396 klu_solve (Symbolic, Numeric, n, 1, (double[]) X, X_offset, Common) ;
397
398 est_new = 0.0 ;
399 for (j = 0 ; j < n ; j++)
400 {
401 /* est_new += ABS (X [j]) ;*/
402 //ABS (abs_value, X [j]) ;
403 abs_value = ABS (X [X_offset + j]) ;
404 est_new += abs_value ;
405 }
406 est_new = 2 * est_new / (3 * n) ;
407 ainv_norm = MAX (est_new, ainv_norm) ;
408
409 /* ---------------------------------------------------------------------- */
410 /* compute estimate of condition number */
411 /* ---------------------------------------------------------------------- */
412
413 Common.condest = ainv_norm * anorm ;
414 return (TRUE) ;
415 }
416
417 /**
418 * Compute the flop count for the LU factorization (in Common.flops)
419 *
420 * @param Symbolic
421 * @param Numeric
422 * @param Common
423 * @return TRUE if successful, FALSE otherwise
424 */
425 public static int klu_flops(KLU_symbolic Symbolic, KLU_numeric Numeric,
426 KLU_common Common)
427 {
428 double flops = 0 ;
429 int[] R, Uip, Llen, Ulen ;
430 /*int[]*/double[] Ui ;
431 double[][] LUbx ;
432 double[] LU ;
433 int k, ulen, p, nk, block, nblocks, k1 ;
434 int Llen_offset = 0, Uip_offset = 0, Ulen_offset = 0 ;
435
436 /* ---------------------------------------------------------------------- */
437 /* check inputs */
438 /* ---------------------------------------------------------------------- */
439
440 if (Common == null)
441 {
442 return (FALSE) ;
443 }
444 Common.flops = EMPTY ;
445 if (Numeric == null || Symbolic == null)
446 {
447 Common.status = KLU_INVALID ;
448 return (FALSE) ;
449 }
450 Common.status = KLU_OK ;
451
452 /* ---------------------------------------------------------------------- */
453 /* get the contents of the Symbolic object */
454 /* ---------------------------------------------------------------------- */
455
456 R = Symbolic.R ;
457 nblocks = Symbolic.nblocks ;
458
459 /* ---------------------------------------------------------------------- */
460 /* get the contents of the Numeric object */
461 /* ---------------------------------------------------------------------- */
462
463 LUbx = (double[][]) Numeric.LUbx ;
464
465 /* ---------------------------------------------------------------------- */
466 /* compute the flop count */
467 /* ---------------------------------------------------------------------- */
468
469 for (block = 0 ; block < nblocks ; block++)
470 {
471 k1 = R [block] ;
472 nk = R [block+1] - k1 ;
473 if (nk > 1)
474 {
475 Llen = Numeric.Llen ;
476 Llen_offset += k1 ;
477 Uip = Numeric.Uip ;
478 Uip_offset += k1 ;
479 Ulen = Numeric.Ulen ;
480 Ulen_offset += k1 ;
481 LU = LUbx [block] ;
482 int[] Ui_offset = new int[1] ;
483 for (k = 0 ; k < nk ; k++)
484 {
485 /* compute kth column of U, and update kth column of A */
486 Ui = GET_I_POINTER (LU, Uip, Uip_offset, Ui_offset, k) ;
487 ulen = Ulen [Ulen_offset + k] ;
488 for (p = 0 ; p < ulen ; p++)
489 {
490 flops += 2 * Llen [Llen_offset + (int) Ui [Ui_offset[0] + p]] ;
491 }
492 /* gather and divide by pivot to get kth column of L */
493 flops += Llen [Llen_offset + k] ;
494 }
495 }
496 }
497 Common.flops = flops ;
498 return (TRUE) ;
499 }
500
501 /**
502 * Compute a really cheap estimate of the reciprocal of the condition number,
503 * condition number, min(abs(diag(U))) / max(abs(diag(U))). If U has a zero
504 * pivot, or a NaN pivot, rcond will be zero. Takes O(n) time.
505 *
506 * @param Symbolic
507 * @param Numeric
508 * @param Common result in Common.rcond
509 * @return TRUE if successful, FALSE otherwise
510 */
511 public static int klu_rcond(KLU_symbolic Symbolic, KLU_numeric Numeric,
512 KLU_common Common)
513 {
514 double ukk, umin = 0, umax = 0 ;
515 double[] Udiag ;
516 int j, n ;
517
518 /* ---------------------------------------------------------------------- */
519 /* check inputs */
520 /* ---------------------------------------------------------------------- */
521
522 if (Common == null)
523 {
524 return (FALSE) ;
525 }
526 if (Symbolic == null)
527 {
528 Common.status = KLU_INVALID ;
529 return (FALSE) ;
530 }
531 if (Numeric == null)
532 {
533 Common.rcond = 0 ;
534 Common.status = KLU_SINGULAR ;
535 return (TRUE) ;
536 }
537 Common.status = KLU_OK ;
538
539 /* ---------------------------------------------------------------------- */
540 /* compute rcond */
541 /* ---------------------------------------------------------------------- */
542
543 n = Symbolic.n ;
544 Udiag = Numeric.Udiag ;
545 for (j = 0 ; j < n ; j++)
546 {
547 /* get the magnitude of the pivot */
548 //ABS (ukk, Udiag [j]) ;
549 ukk = ABS (Udiag [j]) ;
550 if (SCALAR_IS_NAN (ukk) || SCALAR_IS_ZERO (ukk))
551 {
552 /* if NaN, or zero, the rcond is zero */
553 Common.rcond = 0 ;
554 Common.status = KLU_SINGULAR ;
555 return (TRUE) ;
556 }
557 if (j == 0)
558 {
559 /* first pivot entry */
560 umin = ukk ;
561 umax = ukk ;
562 }
563 else
564 {
565 /* subsequent pivots */
566 umin = MIN (umin, ukk) ;
567 umax = MAX (umax, ukk) ;
568 }
569 }
570
571 Common.rcond = umin / umax ;
572 if (SCALAR_IS_NAN (Common.rcond) || SCALAR_IS_ZERO (Common.rcond))
573 {
574 /* this can occur if umin or umax are Inf or NaN */
575 Common.rcond = 0 ;
576 Common.status = KLU_SINGULAR ;
577 }
578 return (TRUE) ;
579 }
580
581 }
A sparse LU factorization algorithm in Java