1 /* ========================================================================== */
2 /* === UMF_2by2 ============================================================= */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Copyright (c) Timothy A. Davis, CISE, */
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */
9 /* -------------------------------------------------------------------------- */
10
11 /* Not user-callable. Computes a row permutation P so that A (P,:) has a
12 * mostly zero-free diagonal, with large entries on the diagonal. It does this
13 * by swapping pairs of rows. Once a row is swapped it is not swapped again.
14 * This is a "cheap" assignment, not a complete max. transversal or
15 * bi-partite matching. It is only a partial matching. For most matrices
16 * for which this algorithm is used, however, the matching is complete (in
17 * UMFPACK this algorithm is used for matrices with roughly symmetric pattern,
18 * and these matrices typically have a mostly-zero-free diagonal to begin with.
19 * This algorithm is not meant to be used on arbitrary unsymmetric matrices
20 * (for those matrices, UMFPACK uses its unsymmetric strategy and does not
21 * use this algorithm).
22 *
23 * Even if incomplete, the matching is usually good enough for UMFPACK's
24 * symmetric strategy, which can easily pivot off the diagonal during numerical
25 * factorization if it finds a weak diagonal entry.
26 *
27 * The algorithms works as follows. First, row scaling factors are computed,
28 * and weak diagonal entries are found. A weak entry is a value A(k,k) whose
29 * absolute value is < tol * max (abs (A (:,k))). For each weak diagonal k in
30 * increasing order of degree in A+A', the algorithm finds an index j such
31 * that A (k,j) and A (j,k) are "large" (greater than or equal to tol times
32 * the largest magnitude in their columns). Row j must also not have already
33 * been swapped. Rows j and k are then swapped. If we come to a diagonal k
34 * that has already been swapped, then it is not modified. This case occurs
35 * for "oxo" pivots:
36 *
37 * k j
38 * k o x
39 * j x o
40 *
41 * which are swapped once to obtain
42 *
43 * k j
44 * j x o
45 * k o x
46 *
47 * These two rows are then not modified any further (A (j,j) was weak, but
48 * after one swap the permuted the jth diagonal entry is strong.
49 *
50 * This algorithm only works on square matrices (real, complex, or pattern-
51 * only). The numerical values are optional. If not present, each entry is
52 * treated as numerically acceptable (tol is ignored), and the algorithm
53 * operates by just using the pattern, not the values. Each column of the
54 * input matrix A must be sorted, with no duplicate entries. The matrix A
55 * can be optionally scaled prior to the numerical test. The matrix A (:,P)
56 * has the same diagonal entries as A (:,P), except in different order. So
57 * the output permutation P can also be used to swap the columns of A.
58 */
59
60 #include "umf_internal.h"
61
62 #ifndef NDEBUG
63 #include "umf_is_permutation.h"
64 #endif
65
66 /* x is "weak" if it is less than ctol. If x or ctol are NaN, then define
67 * x as not "weak". This is a rather arbitrary choice, made to simplify the
68 * computation. On all but a PC with Microsoft C/C++, this test becomes
69 * ((x) - ctol < 0). */
70 #define WEAK(x,ctol) (SCALAR_IS_LTZERO ((x)-(ctol)))
71
72 /* For flag value in Next [col] */
73 #define IS_WEAK -2
74
75 /* ========================================================================== */
76 /* === two_by_two =========================================================== */
77 /* ========================================================================== */
78
two_by_two(Int n2,Int Cp[],Int Ci[],Int Degree[],Int Next[],Int Ri[],Int P[],Int Rp[],Int Head[])79 PRIVATE Int two_by_two /* returns # unmatched weak diagonals */
80 (
81 /* input, not modified */
82 Int n2, /* C is n2-by-n2 */
83 Int Cp [ ], /* size n2+1, column pointers for C */
84 Int Ci [ ], /* size snz = Cp [n2], row indices for C */
85 Int Degree [ ], /* Degree [i] = degree of row i of C+C' */
86
87 /* input, not defined on output */
88 Int Next [ ], /* Next [k] == IS_WEAK if k is a weak diagonal */
89 Int Ri [ ], /* Ri [i] is the length of row i in C */
90
91 /* output, not defined on input */
92 Int P [ ],
93
94 /* workspace, not defined on input or output */
95 Int Rp [ ],
96 Int Head [ ]
97 )
98 {
99 Int deg, newcol, row, col, p, p2, unmatched, k, j, j2, j_best, best, jdiff,
100 jdiff_best, jdeg, jdeg_best, cp, cp1, cp2, rp, rp1, rp2, maxdeg,
101 mindeg ;
102
103 /* ---------------------------------------------------------------------- */
104 /* place weak diagonals in the degree lists */
105 /* ---------------------------------------------------------------------- */
106
107 for (deg = 0 ; deg < n2 ; deg++)
108 {
109 Head [deg] = EMPTY ;
110 }
111
112 maxdeg = 0 ;
113 mindeg = Int_MAX ;
114 for (newcol = n2-1 ; newcol >= 0 ; newcol--)
115 {
116 if (Next [newcol] == IS_WEAK)
117 {
118 /* add this column to the list of weak nodes */
119 DEBUGm1 ((" newcol " ID " has a weak diagonal deg " ID "\n",
120 newcol, deg)) ;
121 deg = Degree [newcol] ;
122 ASSERT (deg >= 0 && deg < n2) ;
123 Next [newcol] = Head [deg] ;
124 Head [deg] = newcol ;
125 maxdeg = MAX (maxdeg, deg) ;
126 mindeg = MIN (mindeg, deg) ;
127 }
128 }
129
130 /* ---------------------------------------------------------------------- */
131 /* construct R = C' (C = strong entries in pruned submatrix) */
132 /* ---------------------------------------------------------------------- */
133
134 /* Ri [0..n2-1] is the length of each row of R */
135 /* use P as temporary pointer into the row form of R [ */
136 Rp [0] = 0 ;
137 for (row = 0 ; row < n2 ; row++)
138 {
139 Rp [row+1] = Rp [row] + Ri [row] ;
140 P [row] = Rp [row] ;
141 }
142 /* Ri no longer needed for row counts */
143
144 /* all entries in C are strong */
145 for (col = 0 ; col < n2 ; col++)
146 {
147 p2 = Cp [col+1] ;
148 for (p = Cp [col] ; p < p2 ; p++)
149 {
150 /* place the column index in row = Ci [p] */
151 Ri [P [Ci [p]]++] = col ;
152 }
153 }
154
155 /* contents of P no longer needed ] */
156
157 #ifndef NDEBUG
158 DEBUG0 (("==================R: row form of strong entries in A:\n")) ;
159 UMF_dump_col_matrix ((double *) NULL,
160 #ifdef COMPLEX
161 (double *) NULL,
162 #endif
163 Ri, Rp, n2, n2, Rp [n2]) ;
164 #endif
165 ASSERT (AMD_valid (n2, n2, Rp, Ri) == AMD_OK) ;
166
167 /* ---------------------------------------------------------------------- */
168 /* for each weak diagonal, find a pair of strong off-diagonal entries */
169 /* ---------------------------------------------------------------------- */
170
171 for (row = 0 ; row < n2 ; row++)
172 {
173 P [row] = EMPTY ;
174 }
175
176 unmatched = 0 ;
177 best = EMPTY ;
178 jdiff = EMPTY ;
179 jdeg = EMPTY ;
180
181 for (deg = mindeg ; deg <= maxdeg ; deg++)
182 {
183 /* find the next weak diagonal of lowest degree */
184 DEBUGm2 (("---------------------------------- Deg: " ID "\n", deg)) ;
185 for (k = Head [deg] ; k != EMPTY ; k = Next [k])
186 {
187 DEBUGm2 (("k: " ID "\n", k)) ;
188 if (P [k] == EMPTY)
189 {
190 /* C (k,k) is a weak diagonal entry. Find an index j != k such
191 * that C (j,k) and C (k,j) are both strong, and also such
192 * that Degree [j] is minimized. In case of a tie, pick
193 * the smallest index j. C and R contain the pattern of
194 * strong entries only.
195 *
196 * Note that row k of R and column k of C are both sorted. */
197
198 DEBUGm4 (("===== Weak diagonal k = " ID "\n", k)) ;
199 DEBUG1 (("Column k of C:\n")) ;
200 for (p = Cp [k] ; p < Cp [k+1] ; p++)
201 {
202 DEBUG1 ((" " ID ": deg " ID "\n", Ci [p], Degree [Ci [p]]));
203 }
204 DEBUG1 (("Row k of R (strong entries only):\n")) ;
205 for (p = Rp [k] ; p < Rp [k+1] ; p++)
206 {
207 DEBUG1 ((" " ID ": deg " ID "\n", Ri [p], Degree [Ri [p]]));
208 }
209
210 /* no (C (k,j), C (j,k)) pair exists yet */
211 j_best = EMPTY ;
212 jdiff_best = Int_MAX ;
213 jdeg_best = Int_MAX ;
214
215 /* pointers into column k (including values) */
216 cp1 = Cp [k] ;
217 cp2 = Cp [k+1] ;
218 cp = cp1 ;
219
220 /* pointers into row k (strong entries only, no values) */
221 rp1 = Rp [k] ;
222 rp2 = Rp [k+1] ;
223 rp = rp1 ;
224
225 /* while entries searched in column k and row k */
226 while (TRUE)
227 {
228
229 if (cp >= cp2)
230 {
231 /* no more entries in this column */
232 break ;
233 }
234
235 /* get C (j,k), which is strong */
236 j = Ci [cp] ;
237
238 if (rp >= rp2)
239 {
240 /* no more entries in this column */
241 break ;
242 }
243
244 /* get R (k,j2), which is strong */
245 j2 = Ri [rp] ;
246
247 if (j < j2)
248 {
249 /* C (j,k) is strong, but R (k,j) is not strong */
250 cp++ ;
251 continue ;
252 }
253
254 if (j2 < j)
255 {
256 /* C (k,j2) is strong, but R (j2,k) is not strong */
257 rp++ ;
258 continue ;
259 }
260
261 /* j == j2: C (j,k) is strong and R (k,j) is strong */
262
263 best = FALSE ;
264
265 if (P [j] == EMPTY)
266 {
267 /* j has not yet been matched */
268 jdeg = Degree [j] ;
269 jdiff = SCALAR_ABS (k-j) ;
270
271 DEBUG1 (("Try candidate j " ID " deg " ID " diff " ID
272 "\n", j, jdeg, jdiff)) ;
273
274 if (j_best == EMPTY)
275 {
276 /* this is the first candidate seen */
277 DEBUG1 ((" first\n")) ;
278 best = TRUE ;
279 }
280 else
281 {
282 if (jdeg < jdeg_best)
283 {
284 /* the degree of j is best seen so far. */
285 DEBUG1 ((" least degree\n")) ;
286 best = TRUE ;
287 }
288 else if (jdeg == jdeg_best)
289 {
290 /* degree of j and j_best are the same */
291 /* tie break by nearest node number */
292 if (jdiff < jdiff_best)
293 {
294 DEBUG1 ((" tie degree, closer\n")) ;
295 best = TRUE ;
296 }
297 else if (jdiff == jdiff_best)
298 {
299 /* |j-k| = |j_best-k|. For any given k
300 * and j_best there is only one other j
301 * than can be just as close as j_best.
302 * Tie break by picking the smaller of
303 * j and j_best */
304 DEBUG1 ((" tie degree, as close\n"));
305 best = j < j_best ;
306 }
307 }
308 else
309 {
310 /* j has higher degree than best so far */
311 best = FALSE ;
312 }
313 }
314 }
315
316 if (best)
317 {
318 /* j is best match for k */
319 /* found a strong pair, A (j,k) and A (k,j) */
320 DEBUG1 ((" --- Found pair k: " ID " j: " ID
321 " jdeg: " ID " jdiff: " ID "\n",
322 k, j, jdeg, jdiff)) ;
323 ASSERT (jdiff != EMPTY) ;
324 ASSERT (jdeg != EMPTY) ;
325 j_best = j ;
326 jdeg_best = jdeg ;
327 jdiff_best = jdiff ;
328 }
329
330 /* get the next entries in column k and row k */
331 cp++ ;
332 rp++ ;
333 }
334
335 /* save the pair (j,k), if we found one */
336 if (j_best != EMPTY)
337 {
338 j = j_best ;
339 DEBUGm4 ((" --- best pair j: " ID " for k: " ID "\n", j, k)) ;
340 P [k] = j ;
341 P [j] = k ;
342 }
343 else
344 {
345 /* no match was found for k */
346 unmatched++ ;
347 }
348 }
349 }
350 }
351
352 /* ---------------------------------------------------------------------- */
353 /* finalize the row permutation, P */
354 /* ---------------------------------------------------------------------- */
355
356 for (k = 0 ; k < n2 ; k++)
357 {
358 if (P [k] == EMPTY)
359 {
360 P [k] = k ;
361 }
362 }
363 ASSERT (UMF_is_permutation (P, Rp, n2, n2)) ;
364
365 return (unmatched) ;
366 }
367
368
369 /* ========================================================================== */
370 /* === UMF_2by2 ============================================================= */
371 /* ========================================================================== */
372
UMF_2by2(Int n,const Int Ap[],const Int Ai[],const double Ax[],const double Az[],double tol,Int scale,Int Cperm1[],Int Rperm1[],Int InvRperm1[],Int n1,Int nempty,Int Degree[],Int P[],Int * p_nweak,Int * p_unmatched,Int Ri[],Int Rp[],double Rs[],Int Head[],Int Next[],Int Ci[],Int Cp[])373 GLOBAL void UMF_2by2
374 (
375 /* input, not modified: */
376 Int n, /* A is n-by-n */
377 const Int Ap [ ], /* size n+1 */
378 const Int Ai [ ], /* size nz = Ap [n] */
379 const double Ax [ ], /* size nz if present */
380 #ifdef COMPLEX
381 const double Az [ ], /* size nz if present */
382 #endif
383 double tol, /* tolerance for determining whether or not an
384 * entry is numerically acceptable. If tol <= 0
385 * then all numerical values ignored. */
386 Int scale, /* scaling to perform (none, sum, or max) */
387 Int Cperm1 [ ], /* singleton permutations */
388 #ifndef NDEBUG
389 Int Rperm1 [ ], /* not needed, since Rperm1 = Cperm1 for submatrix S */
390 #endif
391 Int InvRperm1 [ ], /* inverse of Rperm1 */
392 Int n1, /* number of singletons */
393 Int nempty, /* number of empty rows/cols */
394
395 /* input, contents undefined on output: */
396 Int Degree [ ], /* Degree [j] is the number of off-diagonal
397 * entries in row/column j of S+S', where
398 * where S = A (Cperm1 [n1..], Rperm1 [n1..]).
399 * Note that S is not used, nor formed. */
400
401 /* output: */
402 Int P [ ], /* P [k] = i means original row i is kth row in S(P,:)
403 * where S = A (Cperm1 [n1..], Rperm1 [n1..]) */
404 Int *p_nweak,
405 Int *p_unmatched,
406
407 /* workspace (not defined on input or output): */
408 Int Ri [ ], /* of size >= max (nz, n) */
409 Int Rp [ ], /* of size n+1 */
410 double Rs [ ], /* of size n if present. Rs = sum (abs (A),2) or
411 * max (abs (A),2), the sum or max of each row. Unused
412 * if scale is equal to UMFPACK_SCALE_NONE. */
413 Int Head [ ], /* of size n. Head pointers for bucket sort */
414 Int Next [ ], /* of size n. Next pointers for bucket sort */
415 Int Ci [ ], /* size nz */
416 Int Cp [ ] /* size n+1 */
417 )
418 {
419
420 /* ---------------------------------------------------------------------- */
421 /* local variables */
422 /* ---------------------------------------------------------------------- */
423
424 Entry aij ;
425 double cmax, value, rs, ctol, dvalue ;
426 Int k, p, row, col, do_values, do_sum, do_max, do_scale, nweak, weak,
427 p1, p2, dfound, unmatched, n2, oldrow, newrow, oldcol, newcol, pp ;
428 #ifdef COMPLEX
429 Int split = SPLIT (Az) ;
430 #endif
431 #ifndef NRECIPROCAL
432 Int do_recip = FALSE ;
433 #endif
434
435 #ifndef NDEBUG
436 /* UMF_debug += 99 ; */
437 DEBUGm3 (("\n ==================================UMF_2by2: tol %g\n", tol)) ;
438 ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
439 for (k = n1 ; k < n - nempty ; k++)
440 {
441 ASSERT (Cperm1 [k] == Rperm1 [k]) ;
442 }
443 #endif
444
445 /* ---------------------------------------------------------------------- */
446 /* determine scaling options */
447 /* ---------------------------------------------------------------------- */
448
449 /* use the values, but only if they are present */
450 /* ignore the values if tol <= 0 */
451 do_values = (tol > 0) && (Ax != (double *) NULL) ;
452 if (do_values && (Rs != (double *) NULL))
453 {
454 do_sum = (scale == UMFPACK_SCALE_SUM) ;
455 do_max = (scale == UMFPACK_SCALE_MAX) ;
456 }
457 else
458 {
459 /* no scaling */
460 do_sum = FALSE ;
461 do_max = FALSE ;
462 }
463 do_scale = do_max || do_sum ;
464 DEBUGm3 (("do_values " ID " do_sum " ID " do_max " ID " do_scale " ID "\n",
465 do_values, do_sum, do_max, do_scale)) ;
466
467 /* ---------------------------------------------------------------------- */
468 /* compute the row scaling, if requested */
469 /* ---------------------------------------------------------------------- */
470
471 /* see also umf_kernel_init */
472
473 if (do_scale)
474 {
475 #ifndef NRECIPROCAL
476 double rsmin ;
477 #endif
478 for (row = 0 ; row < n ; row++)
479 {
480 Rs [row] = 0.0 ;
481 }
482 for (col = 0 ; col < n ; col++)
483 {
484 p2 = Ap [col+1] ;
485 for (p = Ap [col] ; p < p2 ; p++)
486 {
487 row = Ai [p] ;
488 ASSIGN (aij, Ax, Az, p, split) ;
489 APPROX_ABS (value, aij) ;
490 rs = Rs [row] ;
491 if (!SCALAR_IS_NAN (rs))
492 {
493 if (SCALAR_IS_NAN (value))
494 {
495 /* if any entry in a row is NaN, then the scale factor
496 * for the row is NaN. It will be set to 1 later. */
497 Rs [row] = value ;
498 }
499 else if (do_max)
500 {
501 Rs [row] = MAX (rs, value) ;
502 }
503 else
504 {
505 Rs [row] += value ;
506 }
507 }
508 }
509 }
510 #ifndef NRECIPROCAL
511 rsmin = Rs [0] ;
512 if (SCALAR_IS_ZERO (rsmin) || SCALAR_IS_NAN (rsmin))
513 {
514 rsmin = 1.0 ;
515 }
516 #endif
517 for (row = 0 ; row < n ; row++)
518 {
519 /* do not scale an empty row, or a row with a NaN */
520 rs = Rs [row] ;
521 if (SCALAR_IS_ZERO (rs) || SCALAR_IS_NAN (rs))
522 {
523 Rs [row] = 1.0 ;
524 }
525 #ifndef NRECIPROCAL
526 rsmin = MIN (rsmin, Rs [row]) ;
527 #endif
528 }
529
530 #ifndef NRECIPROCAL
531 /* multiply by the reciprocal if Rs is not too small */
532 do_recip = (rsmin >= RECIPROCAL_TOLERANCE) ;
533 if (do_recip)
534 {
535 /* invert the scale factors */
536 for (row = 0 ; row < n ; row++)
537 {
538 Rs [row] = 1.0 / Rs [row] ;
539 }
540 }
541 #endif
542 }
543
544 /* ---------------------------------------------------------------------- */
545 /* compute the max in each column and find diagonal */
546 /* ---------------------------------------------------------------------- */
547
548 nweak = 0 ;
549
550 #ifndef NDEBUG
551 for (k = 0 ; k < n ; k++)
552 {
553 ASSERT (Rperm1 [k] >= 0 && Rperm1 [k] < n) ;
554 ASSERT (InvRperm1 [Rperm1 [k]] == k) ;
555 }
556 #endif
557
558 n2 = n - n1 - nempty ;
559
560 /* use Ri to count the number of strong entries in each row */
561 for (row = 0 ; row < n2 ; row++)
562 {
563 Ri [row] = 0 ;
564 }
565
566 pp = 0 ;
567 ctol = 0 ;
568 dvalue = 1 ;
569
570 /* construct C = pruned submatrix, strong values only, column form */
571
572 for (k = n1 ; k < n - nempty ; k++)
573 {
574 oldcol = Cperm1 [k] ;
575 newcol = k - n1 ;
576 Next [newcol] = EMPTY ;
577 DEBUGm1 (("Column " ID " newcol " ID " oldcol " ID "\n", k, newcol, oldcol)) ;
578
579 Cp [newcol] = pp ;
580
581 dfound = FALSE ;
582 p1 = Ap [oldcol] ;
583 p2 = Ap [oldcol+1] ;
584 if (do_values)
585 {
586 cmax = 0 ;
587 dvalue = 0 ;
588
589 if (!do_scale)
590 {
591 /* no scaling */
592 for (p = p1 ; p < p2 ; p++)
593 {
594 oldrow = Ai [p] ;
595 ASSERT (oldrow >= 0 && oldrow < n) ;
596 newrow = InvRperm1 [oldrow] - n1 ;
597 ASSERT (newrow >= -n1 && newrow < n2) ;
598 if (newrow < 0) continue ;
599 ASSIGN (aij, Ax, Az, p, split) ;
600 APPROX_ABS (value, aij) ;
601 /* if either cmax or value is NaN, define cmax as NaN */
602 if (!SCALAR_IS_NAN (cmax))
603 {
604 if (SCALAR_IS_NAN (value))
605 {
606 cmax = value ;
607 }
608 else
609 {
610 cmax = MAX (cmax, value) ;
611 }
612 }
613 if (oldrow == oldcol)
614 {
615 /* we found the diagonal entry in this column */
616 dvalue = value ;
617 dfound = TRUE ;
618 ASSERT (newrow == newcol) ;
619 }
620 }
621 }
622 #ifndef NRECIPROCAL
623 else if (do_recip)
624 {
625 /* multiply by the reciprocal */
626 for (p = p1 ; p < p2 ; p++)
627 {
628 oldrow = Ai [p] ;
629 ASSERT (oldrow >= 0 && oldrow < n) ;
630 newrow = InvRperm1 [oldrow] - n1 ;
631 ASSERT (newrow >= -n1 && newrow < n2) ;
632 if (newrow < 0) continue ;
633 ASSIGN (aij, Ax, Az, p, split) ;
634 APPROX_ABS (value, aij) ;
635 value *= Rs [oldrow] ;
636 /* if either cmax or value is NaN, define cmax as NaN */
637 if (!SCALAR_IS_NAN (cmax))
638 {
639 if (SCALAR_IS_NAN (value))
640 {
641 cmax = value ;
642 }
643 else
644 {
645 cmax = MAX (cmax, value) ;
646 }
647 }
648 if (oldrow == oldcol)
649 {
650 /* we found the diagonal entry in this column */
651 dvalue = value ;
652 dfound = TRUE ;
653 ASSERT (newrow == newcol) ;
654 }
655 }
656 }
657 #endif
658 else
659 {
660 /* divide instead */
661 for (p = p1 ; p < p2 ; p++)
662 {
663 oldrow = Ai [p] ;
664 ASSERT (oldrow >= 0 && oldrow < n) ;
665 newrow = InvRperm1 [oldrow] - n1 ;
666 ASSERT (newrow >= -n1 && newrow < n2) ;
667 if (newrow < 0) continue ;
668 ASSIGN (aij, Ax, Az, p, split) ;
669 APPROX_ABS (value, aij) ;
670 value /= Rs [oldrow] ;
671 /* if either cmax or value is NaN, define cmax as NaN */
672 if (!SCALAR_IS_NAN (cmax))
673 {
674 if (SCALAR_IS_NAN (value))
675 {
676 cmax = value ;
677 }
678 else
679 {
680 cmax = MAX (cmax, value) ;
681 }
682 }
683 if (oldrow == oldcol)
684 {
685 /* we found the diagonal entry in this column */
686 dvalue = value ;
687 dfound = TRUE ;
688 ASSERT (newrow == newcol) ;
689 }
690 }
691 }
692
693 ctol = tol * cmax ;
694 DEBUGm1 ((" cmax col " ID " %g ctol %g\n", oldcol, cmax, ctol)) ;
695 }
696 else
697 {
698 for (p = p1 ; p < p2 ; p++)
699 {
700 oldrow = Ai [p] ;
701 ASSERT (oldrow >= 0 && oldrow < n) ;
702 newrow = InvRperm1 [oldrow] - n1 ;
703 ASSERT (newrow >= -n1 && newrow < n2) ;
704 if (newrow < 0) continue ;
705 Ci [pp++] = newrow ;
706 if (oldrow == oldcol)
707 {
708 /* we found the diagonal entry in this column */
709 ASSERT (newrow == newcol) ;
710 dfound = TRUE ;
711 }
712 /* count the entries in each column */
713 Ri [newrow]++ ;
714 }
715 }
716
717 /* ------------------------------------------------------------------ */
718 /* flag the weak diagonals */
719 /* ------------------------------------------------------------------ */
720
721 if (!dfound)
722 {
723 /* no diagonal entry present */
724 weak = TRUE ;
725 }
726 else
727 {
728 /* diagonal entry is present, check its value */
729 weak = (do_values) ? WEAK (dvalue, ctol) : FALSE ;
730 }
731 if (weak)
732 {
733 /* flag this column as weak */
734 DEBUG0 (("Weak!\n")) ;
735 Next [newcol] = IS_WEAK ;
736 nweak++ ;
737 }
738
739 /* ------------------------------------------------------------------ */
740 /* count entries in each row that are not numerically weak */
741 /* ------------------------------------------------------------------ */
742
743 if (do_values)
744 {
745 if (!do_scale)
746 {
747 /* no scaling */
748 for (p = p1 ; p < p2 ; p++)
749 {
750 oldrow = Ai [p] ;
751 newrow = InvRperm1 [oldrow] - n1 ;
752 if (newrow < 0) continue ;
753 ASSIGN (aij, Ax, Az, p, split) ;
754 APPROX_ABS (value, aij) ;
755 weak = WEAK (value, ctol) ;
756 if (!weak)
757 {
758 DEBUG0 ((" strong: row " ID ": %g\n", oldrow, value)) ;
759 Ci [pp++] = newrow ;
760 Ri [newrow]++ ;
761 }
762 }
763 }
764 #ifndef NRECIPROCAL
765 else if (do_recip)
766 {
767 /* multiply by the reciprocal */
768 for (p = p1 ; p < p2 ; p++)
769 {
770 oldrow = Ai [p] ;
771 newrow = InvRperm1 [oldrow] - n1 ;
772 if (newrow < 0) continue ;
773 ASSIGN (aij, Ax, Az, p, split) ;
774 APPROX_ABS (value, aij) ;
775 value *= Rs [oldrow] ;
776 weak = WEAK (value, ctol) ;
777 if (!weak)
778 {
779 DEBUG0 ((" strong: row " ID ": %g\n", oldrow, value)) ;
780 Ci [pp++] = newrow ;
781 Ri [newrow]++ ;
782 }
783 }
784 }
785 #endif
786 else
787 {
788 /* divide instead */
789 for (p = p1 ; p < p2 ; p++)
790 {
791 oldrow = Ai [p] ;
792 newrow = InvRperm1 [oldrow] - n1 ;
793 if (newrow < 0) continue ;
794 ASSIGN (aij, Ax, Az, p, split) ;
795 APPROX_ABS (value, aij) ;
796 value /= Rs [oldrow] ;
797 weak = WEAK (value, ctol) ;
798 if (!weak)
799 {
800 DEBUG0 ((" strong: row " ID ": %g\n", oldrow, value)) ;
801 Ci [pp++] = newrow ;
802 Ri [newrow]++ ;
803 }
804 }
805 }
806 }
807 }
808 Cp [n2] = pp ;
809 ASSERT (AMD_valid (n2, n2, Cp, Ci) == AMD_OK) ;
810
811 if (nweak == 0)
812 {
813 /* nothing to do, quick return */
814 DEBUGm2 (("\n =============================UMF_2by2: quick return\n")) ;
815 for (k = 0 ; k < n ; k++)
816 {
817 P [k] = k ;
818 }
819 *p_nweak = 0 ;
820 *p_unmatched = 0 ;
821 return ;
822 }
823
824 #ifndef NDEBUG
825 for (k = 0 ; k < n2 ; k++)
826 {
827 P [k] = EMPTY ;
828 }
829 for (k = 0 ; k < n2 ; k++)
830 {
831 ASSERT (Degree [k] >= 0 && Degree [k] < n2) ;
832 }
833 #endif
834
835 /* ---------------------------------------------------------------------- */
836 /* find the 2-by-2 permutation */
837 /* ---------------------------------------------------------------------- */
838
839 /* The matrix S is now mapped to the index range 0 to n2-1. We have
840 * S = A (Rperm [n1 .. n-nempty-1], Cperm [n1 .. n-nempty-1]), and then
841 * C = pattern of strong entries in S. A weak diagonal k in S is marked
842 * with Next [k] = IS_WEAK. */
843
844 unmatched = two_by_two (n2, Cp, Ci, Degree, Next, Ri, P, Rp, Head) ;
845
846 /* ---------------------------------------------------------------------- */
847
848 *p_nweak = nweak ;
849 *p_unmatched = unmatched ;
850
851 #ifndef NDEBUG
852 DEBUGm4 (("UMF_2by2: weak " ID " unmatched " ID "\n", nweak, unmatched)) ;
853 for (row = 0 ; row < n ; row++)
854 {
855 DEBUGm2 (("P [" ID "] = " ID "\n", row, P [row])) ;
856 }
857 DEBUGm2 (("\n =============================UMF_2by2: done\n\n")) ;
858 #endif
859 }
860