← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/umfpack/umf_row_search.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/pysparse/umfpack/umf_row_search.c
1
 
/* ========================================================================== */
2
 
/* === UMF_row_search ======================================================= */
3
 
/* ========================================================================== */
4
 
 
5
 
/* -------------------------------------------------------------------------- */
6
 
/* UMFPACK Version 4.1 (Apr. 30, 2003), Copyright (c) 2003 by Timothy A.      */
7
 
/* Davis.  All Rights Reserved.  See ../README for License.                   */
8
 
/* email: davis@cise.ufl.edu    CISE Department, Univ. of Florida.            */
9
 
/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
10
 
/* -------------------------------------------------------------------------- */
11
 
 
12
 
/*
13
 
    Find two candidate pivot rows in a column: the best one in the front,
14
 
    and the best one not in the front.  Return the two pivot row patterns and
15
 
    their exact degrees.  Called by UMF_local_search.
16
 
 
17
 
    Returns UMFPACK_OK if successful, or UMFPACK_WARNING_singular_matrix or
18
 
    UMFPACK_ERROR_different_pattern if not.
19
 
 
20
 
*/
21
 
 
22
 
#include "umf_internal.h"
23
 
#include "umf_row_search.h"
24
 
 
25
 
GLOBAL Int UMF_row_search
26
 
(
27
 
    NumericType *Numeric,
28
 
    WorkType *Work,
29
 
    SymbolicType *Symbolic,
30
 
    Int cdeg0,                  /* length of column in Front */
31
 
    Int cdeg1,                  /* length of column outside Front */
32
 
    const Int Pattern [ ],      /* pattern of column, Pattern [0..cdeg1 -1] */
33
 
    const Int Pos [ ],          /* Pos [Pattern [0..cdeg1 -1]] = 0..cdeg1 -1 */
34
 
    Int pivrow [2],             /* pivrow [IN] and pivrow [OUT] */
35
 
    Int rdeg [2],               /* rdeg [IN] and rdeg [OUT] */
36
 
    Int W_i [ ],                /* pattern of pivrow [IN], */
37
 
                                /* either Fcols or Woi */
38
 
    Int W_o [ ],                /* pattern of pivrow [OUT], */
39
 
                                /* either Wio or Woo */
40
 
    Int prior_pivrow [2],       /* the two other rows just scanned, if any */
41
 
    const Entry Wxy [ ],        /* numerical values Wxy [0..cdeg1-1],
42
 
                                   either Wx or Wy */
43
 
 
44
 
    Int pivcol,                 /* the candidate column being searched */
45
 
    Int freebie [ ]
46
 
)
47
 
{
48
 
 
49
 
    /* ---------------------------------------------------------------------- */
50
 
    /* local variables */
51
 
    /* ---------------------------------------------------------------------- */
52
 
 
53
 
    double maxval, toler, toler2, value, pivot [2] ;
54
 
    Int i, row, deg, col, *Frpos, fnrows, *E, j, ncols, *Cols, *Rows,
55
 
        e, f, Wrpflag, *Fcpos, fncols, tpi, max_rdeg, nans_in_col, was_offdiag,
56
 
        diag_row, prefer_diagonal, *Wrp, found, *Diagonal_map ;
57
 
    Tuple *tp, *tpend, *tp1, *tp2 ;
58
 
    Unit *Memory, *p ;
59
 
    Element *ep ;
60
 
    Int *Row_tuples, *Row_degree, *Row_tlen ;
61
 
 
62
 
#ifndef NDEBUG
63
 
    Int *Col_degree ;
64
 
    DEBUG2 (("Row_search:\n")) ;
65
 
    for (i = 0 ; i < cdeg1 ; i++)
66
 
    {
67
 
        row = Pattern [i] ;
68
 
        DEBUG4 (("   row: "ID"\n", row)) ;
69
 
        ASSERT (row >= 0 && row < Numeric->n_row) ;
70
 
        ASSERT (i == Pos [row]) ;
71
 
    }
72
 
    /* If row is not in Pattern [0..cdeg1-1], then Pos [row] == EMPTY */
73
 
    if (UMF_debug > 0 || Numeric->n_row < 1000)
74
 
    {
75
 
        Int cnt = cdeg1 ;
76
 
        DEBUG4 (("Scan all rows:\n")) ;
77
 
        for (row = 0 ; row < Numeric->n_row ; row++)
78
 
        {
79
 
            if (Pos [row] < 0)
80
 
            {
81
 
                cnt++ ;
82
 
            }
83
 
            else
84
 
            {
85
 
                DEBUG4 (("   row: "ID" pos "ID"\n", row, Pos [row])) ;
86
 
            }
87
 
        }
88
 
        ASSERT (cnt == Numeric->n_row) ;
89
 
    }
90
 
    Col_degree = Numeric->Cperm ;   /* for NON_PIVOTAL_COL macro only */
91
 
    ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;
92
 
    ASSERT (NON_PIVOTAL_COL (pivcol)) ;
93
 
#endif
94
 
 
95
 
    pivot [IN] = 0. ;
96
 
    pivot [OUT] = 0. ;
97
 
 
98
 
    /* ---------------------------------------------------------------------- */
99
 
    /* get parameters */
100
 
    /* ---------------------------------------------------------------------- */
101
 
 
102
 
    Row_degree = Numeric->Rperm ;
103
 
    Row_tuples = Numeric->Uip ;
104
 
    Row_tlen   = Numeric->Uilen ;
105
 
    Wrp = Work->Wrp ;
106
 
    Frpos = Work->Frpos ;
107
 
    E = Work->E ;
108
 
    Memory = Numeric->Memory ;
109
 
    fnrows = Work->fnrows ;
110
 
 
111
 
    prefer_diagonal = Symbolic->prefer_diagonal ;
112
 
    Diagonal_map = Work->Diagonal_map ;
113
 
 
114
 
    if (Diagonal_map)
115
 
    {
116
 
        diag_row = Diagonal_map [pivcol] ;
117
 
        was_offdiag = diag_row < 0 ;
118
 
        if (was_offdiag)
119
 
        {
120
 
            /* the "diagonal" entry in this column was permuted here by an
121
 
             * earlier pivot choice.  The tighter off-diagonal tolerance will
122
 
             * be used instead of the symmetric tolerance. */
123
 
            diag_row = FLIP (diag_row) ;
124
 
        }
125
 
        ASSERT (diag_row >= 0 && diag_row < Numeric->n_row) ;
126
 
    }
127
 
    else
128
 
    {
129
 
        diag_row = EMPTY ;      /* unused */
130
 
        was_offdiag = EMPTY ;   /* unused */
131
 
    }
132
 
 
133
 
    /* pivot row degree cannot exceed max_rdeg */
134
 
    max_rdeg = Work->fncols_max ;
135
 
 
136
 
    /* ---------------------------------------------------------------------- */
137
 
    /* scan pivot column for candidate rows */
138
 
    /* ---------------------------------------------------------------------- */
139
 
 
140
 
    maxval = 0.0 ;
141
 
    nans_in_col = FALSE ;
142
 
 
143
 
    for (i = 0 ; i < cdeg1 ; i++)
144
 
    {
145
 
        APPROX_ABS (value, Wxy [i]) ;
146
 
        if (SCALAR_IS_NAN (value))
147
 
        {
148
 
            nans_in_col = TRUE ;
149
 
            maxval = value ;
150
 
            break ;
151
 
        }
152
 
        /* This test can now ignore the NaN case: */
153
 
        maxval = MAX (maxval, value) ;
154
 
    }
155
 
 
156
 
    /* if maxval is zero, the matrix is numerically singular */
157
 
 
158
 
    toler = Numeric->relpt * maxval ;
159
 
    toler2 = Numeric->relpt2 * maxval ;
160
 
    toler2 = was_offdiag ? toler : toler2 ;
161
 
 
162
 
    DEBUG5 (("Row_search begins [ maxval %g toler %g %g\n",
163
 
        maxval, toler, toler2)) ;
164
 
    if (SCALAR_IS_NAN (toler) || SCALAR_IS_NAN (toler2))
165
 
    {
166
 
        nans_in_col = TRUE ;
167
 
    }
168
 
 
169
 
    if (!nans_in_col)
170
 
    {
171
 
 
172
 
        /* look for the diagonal entry, if it exists */
173
 
        found = FALSE ;
174
 
        ASSERT (!SCALAR_IS_NAN (toler)) ;
175
 
 
176
 
        if (prefer_diagonal)
177
 
        {
178
 
            ASSERT (diag_row != EMPTY) ;
179
 
            i = Pos [diag_row] ;
180
 
            if (i >= 0)
181
 
            {
182
 
                double a ;
183
 
                ASSERT (i < cdeg1) ;
184
 
                ASSERT (diag_row == Pattern [i]) ;
185
 
 
186
 
                APPROX_ABS (a, Wxy [i]) ;
187
 
 
188
 
                ASSERT (!SCALAR_IS_NAN (a)) ;
189
 
                ASSERT (!SCALAR_IS_NAN (toler2)) ;
190
 
 
191
 
                if (SCALAR_IS_NONZERO (a) && a >= toler2)
192
 
                {
193
 
                    /* found it! */
194
 
                    DEBUG3 (("Symmetric pivot: "ID" "ID"\n", pivcol, diag_row));
195
 
                    found = TRUE ;
196
 
                    if (Frpos [diag_row] >= 0 && Frpos [diag_row] < fnrows)
197
 
                    {
198
 
                        pivrow [IN] = diag_row ;
199
 
                        pivrow [OUT] = EMPTY ;
200
 
                    }
201
 
                    else
202
 
                    {
203
 
                        pivrow [IN] = EMPTY ;
204
 
                        pivrow [OUT] = diag_row ;
205
 
                    }
206
 
                }
207
 
            }
208
 
        }
209
 
 
210
 
        /* either no diagonal found, or we didn't look for it */
211
 
        if (!found)
212
 
        {
213
 
            if (cdeg0 > 0)
214
 
            {
215
 
 
216
 
                /* this is a column in the front */
217
 
                for (i = 0 ; i < cdeg0 ; i++)
218
 
                {
219
 
                    double a ;
220
 
                    APPROX_ABS (a, Wxy [i]) ;
221
 
                    ASSERT (!SCALAR_IS_NAN (a)) ;
222
 
                    ASSERT (!SCALAR_IS_NAN (toler)) ;
223
 
                    if (SCALAR_IS_NONZERO (a) && a >= toler)
224
 
                    {
225
 
                        row = Pattern [i] ;
226
 
                        deg = Row_degree [row] ;
227
 
#ifndef NDEBUG
228
 
                        DEBUG6 ((ID" candidate row "ID" deg "ID" absval %g\n",
229
 
                            i, row, deg, a)) ;
230
 
                        UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
231
 
#endif
232
 
                        ASSERT (Frpos [row] >= 0 && Frpos [row] < fnrows) ;
233
 
                        ASSERT (Frpos [row] == i) ;
234
 
                        /* row is in the current front */
235
 
                        DEBUG4 ((" in front\n")) ;
236
 
                        if (deg < rdeg [IN]
237
 
                            /* break ties by picking the largest entry: */
238
 
                               || (deg == rdeg [IN] && a > pivot [IN])
239
 
                            /* break ties by picking the diagonal entry: */
240
 
                            /* || (deg == rdeg [IN] && row == diag_row) */
241
 
                           )
242
 
                        {
243
 
                            /* best row in front, so far */
244
 
                            pivrow [IN] = row ;
245
 
                            rdeg [IN] = deg ;
246
 
                            pivot [IN] = a ;
247
 
                        }
248
 
                    }
249
 
                }
250
 
                for ( ; i < cdeg1 ; i++)
251
 
                {
252
 
                    double a ;
253
 
                    APPROX_ABS (a, Wxy [i]) ;
254
 
                    ASSERT (!SCALAR_IS_NAN (a)) ;
255
 
                    ASSERT (!SCALAR_IS_NAN (toler)) ;
256
 
                    if (SCALAR_IS_NONZERO (a) && a >= toler)
257
 
                    {
258
 
                        row = Pattern [i] ;
259
 
                        deg = Row_degree [row] ;
260
 
#ifndef NDEBUG
261
 
                        DEBUG6 ((ID" candidate row "ID" deg "ID" absval %g\n",
262
 
                            i, row, deg, a)) ;
263
 
                        UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
264
 
#endif
265
 
                        ASSERT (Frpos [row] == i) ;
266
 
                        /* row is not in the current front */
267
 
                        DEBUG4 ((" NOT in front\n")) ;
268
 
                        if (deg < rdeg [OUT]
269
 
                            /* break ties by picking the largest entry: */
270
 
                               || (deg == rdeg [OUT] && a > pivot [OUT])
271
 
                            /* break ties by picking the diagonal entry: */
272
 
                            /* || (deg == rdeg [OUT] && row == diag_row) */
273
 
                           )
274
 
                        {
275
 
                            /* best row not in front, so far */
276
 
                            pivrow [OUT] = row ;
277
 
                            rdeg [OUT] = deg ;
278
 
                            pivot [OUT] = a ;
279
 
                        }
280
 
                    }
281
 
                }
282
 
 
283
 
            }
284
 
            else
285
 
            {
286
 
 
287
 
                /* this column is not in the front */
288
 
                for (i = 0 ; i < cdeg1 ; i++)
289
 
                {
290
 
                    double a ;
291
 
                    APPROX_ABS (a, Wxy [i]) ;
292
 
                    ASSERT (!SCALAR_IS_NAN (a)) ;
293
 
                    ASSERT (!SCALAR_IS_NAN (toler)) ;
294
 
                    if (SCALAR_IS_NONZERO (a) && a >= toler)
295
 
                    {
296
 
                        row = Pattern [i] ;
297
 
                        deg = Row_degree [row] ;
298
 
#ifndef NDEBUG
299
 
                        DEBUG6 ((ID" candidate row "ID" deg "ID" absval %g\n",
300
 
                            i, row, deg, a)) ;
301
 
                        UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
302
 
#endif
303
 
                        if (Frpos [row] >= 0 && Frpos [row] < fnrows)
304
 
                        {
305
 
                            /* row is in the current front */
306
 
                            DEBUG4 ((" in front\n")) ;
307
 
                            if (deg < rdeg [IN]
308
 
                            /* break ties by picking the largest entry: */
309
 
                               || (deg == rdeg [IN] && a > pivot [IN])
310
 
                            /* break ties by picking the diagonal entry: */
311
 
                            /* || (deg == rdeg [IN] && row == diag_row) */
312
 
                               )
313
 
                            {
314
 
                                /* best row in front, so far */
315
 
                                pivrow [IN] = row ;
316
 
                                rdeg [IN] = deg ;
317
 
                                pivot [IN] = a ;
318
 
                            }
319
 
                        }
320
 
                        else
321
 
                        {
322
 
                            /* row is not in the current front */
323
 
                            DEBUG4 ((" NOT in front\n")) ;
324
 
                            if (deg < rdeg [OUT]
325
 
                            /* break ties by picking the largest entry: */
326
 
                               || (deg == rdeg[OUT] && a > pivot [OUT])
327
 
                            /* break ties by picking the diagonal entry: */
328
 
                            /* || (deg == rdeg[OUT] && row == diag_row) */
329
 
                               )
330
 
                            {
331
 
                                /* best row not in front, so far */
332
 
                                pivrow [OUT] = row ;
333
 
                                rdeg [OUT] = deg ;
334
 
                                pivot [OUT] = a ;
335
 
                            }
336
 
                        }
337
 
                    }
338
 
                }
339
 
            }
340
 
        }
341
 
    }
342
 
 
343
 
    /* ---------------------------------------------------------------------- */
344
 
    /* NaN handling */
345
 
    /* ---------------------------------------------------------------------- */
346
 
 
347
 
    /* if cdeg1 > 0 then we must have found a pivot row ... unless NaN's */
348
 
    /* exist.  Try with no numerical tests if no pivot found. */
349
 
 
350
 
    if (cdeg1 > 0 && pivrow [IN] == EMPTY && pivrow [OUT] == EMPTY)
351
 
    {
352
 
        /* cleanup for the NaN case */
353
 
        DEBUG0 (("Found a NaN in pivot column!\n")) ;
354
 
 
355
 
        /* grab the first entry in the pivot column, ignoring degree, */
356
 
        /* numerical stability, and symmetric preference */
357
 
        row = Pattern [0] ;
358
 
        deg = Row_degree [row] ;
359
 
        if (Frpos [row] >= 0 && Frpos [row] < fnrows)
360
 
        {
361
 
            /* row is in the current front */
362
 
            DEBUG4 ((" in front\n")) ;
363
 
            pivrow [IN] = row ;
364
 
            rdeg [IN] = deg ;
365
 
        }
366
 
        else
367
 
        {
368
 
            /* row is not in the current front */
369
 
            DEBUG4 ((" NOT in front\n")) ;
370
 
            pivrow [OUT] = row ;
371
 
            rdeg [OUT] = deg ;
372
 
        }
373
 
 
374
 
        /* We are now guaranteed to have a pivot, no matter how broken */
375
 
        /* (non-IEEE compliant) the underlying numerical operators are. */
376
 
        /* This is particularly a problem for Microsoft compilers (they do */
377
 
        /* not handle NaN's properly). Now try to find a sparser pivot, if */
378
 
        /* possible. */
379
 
 
380
 
        for (i = 1 ; i < cdeg1 ; i++)
381
 
        {
382
 
            row = Pattern [i] ;
383
 
            deg = Row_degree [row] ;
384
 
 
385
 
            if (Frpos [row] >= 0 && Frpos [row] < fnrows)
386
 
            {
387
 
                /* row is in the current front */
388
 
                DEBUG4 ((" in front\n")) ;
389
 
                if (deg < rdeg [IN] || (deg == rdeg [IN] && row == diag_row))
390
 
                {
391
 
                    /* best row in front, so far */
392
 
                    pivrow [IN] = row ;
393
 
                    rdeg [IN] = deg ;
394
 
                }
395
 
            }
396
 
            else
397
 
            {
398
 
                /* row is not in the current front */
399
 
                DEBUG4 ((" NOT in front\n")) ;
400
 
                if (deg < rdeg [OUT] || (deg == rdeg [OUT] && row == diag_row))
401
 
                {
402
 
                    /* best row not in front, so far */
403
 
                    pivrow [OUT] = row ;
404
 
                    rdeg [OUT] = deg ;
405
 
                }
406
 
            }
407
 
        }
408
 
    }
409
 
 
410
 
    /* We found a pivot if there are entries (even zero ones) in pivot col */
411
 
    ASSERT (IMPLIES (cdeg1 > 0, pivrow[IN] != EMPTY || pivrow[OUT] != EMPTY)) ;
412
 
 
413
 
    /* If there are no entries in the pivot column, then no pivot is found */
414
 
    ASSERT (IMPLIES (cdeg1 == 0, pivrow[IN] == EMPTY && pivrow[OUT] == EMPTY)) ;
415
 
 
416
 
    /* ---------------------------------------------------------------------- */
417
 
    /* check for singular matrix */
418
 
    /* ---------------------------------------------------------------------- */
419
 
 
420
 
    if (cdeg1  == 0)
421
 
    {
422
 
        if (fnrows > 0)
423
 
        {
424
 
            /*
425
 
                Get the pivrow [OUT][IN] from the current front.
426
 
                The frontal matrix looks like this:
427
 
 
428
 
                        pivcol[OUT]
429
 
                        |
430
 
                        v
431
 
                x x x x 0   <- so grab this row as the pivrow [OUT][IN].
432
 
                x x x x 0
433
 
                x x x x 0
434
 
                0 0 0 0 0
435
 
 
436
 
                The current frontal matrix has some rows in it.  The degree
437
 
                of the pivcol[OUT] is zero.  The column is empty, and the
438
 
                current front does not contribute to it.
439
 
 
440
 
            */
441
 
            pivrow [IN] = Work->Frows [0] ;
442
 
            DEBUGm4 (("Got zero pivrow[OUT][IN] "ID" from current front\n",
443
 
                pivrow [IN])) ;
444
 
        }
445
 
        else
446
 
        {
447
 
 
448
 
            /*
449
 
                Get a pivot row from the row-merge tree, use as
450
 
                pivrow [OUT][OUT].   pivrow [IN] remains EMPTY.
451
 
                This can only happen if the current front is 0-by-0.
452
 
            */
453
 
 
454
 
            Int *Front_leftmostdesc, *Front_1strow, *Front_new1strow, row1,
455
 
                row2, fleftmost, nfr, n_row, frontid ;
456
 
 
457
 
            ASSERT (Work->fncols == 0) ;
458
 
 
459
 
            Front_leftmostdesc = Symbolic->Front_leftmostdesc ;
460
 
            Front_1strow = Symbolic->Front_1strow ;
461
 
            Front_new1strow = Work->Front_new1strow ;
462
 
            nfr = Symbolic->nfr ;
463
 
            n_row = Numeric->n_row ;
464
 
            frontid = Work->frontid ;
465
 
 
466
 
            DEBUGm4 (("Note: pivcol: "ID" is empty front "ID"\n",
467
 
                pivcol, frontid)) ;
468
 
#ifndef NDEBUG
469
 
            DEBUG1 (("Calling dump rowmerge\n")) ;
470
 
            UMF_dump_rowmerge (Numeric, Symbolic, Work) ;
471
 
#endif
472
 
 
473
 
            /* Row-merge set is the non-pivotal rows in the range */
474
 
            /* Front_new1strow [Front_leftmostdesc [frontid]] to */
475
 
            /* Front_1strow [frontid+1] - 1. */
476
 
            /* If this is empty, then use the empty rows, in the range */
477
 
            /* Front_new1strow [nfr] to n_row-1. */
478
 
            /* If this too is empty, then pivrow [OUT] will be empty. */
479
 
            /* In both cases, update Front_new1strow [...]. */
480
 
 
481
 
            fleftmost = Front_leftmostdesc [frontid] ;
482
 
            row1 = Front_new1strow [fleftmost] ;
483
 
            row2 = Front_1strow [frontid+1] - 1 ;
484
 
            DEBUG1 (("Leftmost: "ID" Rows ["ID" to "ID"] srch ["ID" to "ID"]\n",
485
 
                fleftmost, Front_1strow [frontid], row2, row1, row2)) ;
486
 
 
487
 
            /* look in the range row1 ... row2 */
488
 
            for (row = row1 ; row <= row2 ; row++)
489
 
            {
490
 
                DEBUG3 (("   Row: "ID"\n", row)) ;
491
 
                if (NON_PIVOTAL_ROW (row))
492
 
                {
493
 
                    /* found it */
494
 
                    DEBUG3 (("   Row: "ID" found\n", row)) ;
495
 
                    ASSERT (Frpos [row] == EMPTY) ;
496
 
                    pivrow [OUT] = row ;
497
 
                    DEBUGm4 (("got row merge pivrow %d\n", pivrow [OUT])) ;
498
 
                    break ;
499
 
                }
500
 
            }
501
 
            Front_new1strow [fleftmost] = row ;
502
 
 
503
 
            if (pivrow [OUT] == EMPTY)
504
 
            {
505
 
                /* not found, look in empty row set in "dummy" front */
506
 
                row1 = Front_new1strow [nfr] ;
507
 
                row2 = n_row-1 ;
508
 
                DEBUG3 (("Empty: "ID" Rows ["ID" to "ID"] srch["ID" to "ID"]\n",
509
 
                    nfr, Front_1strow [nfr], row2, row1, row2)) ;
510
 
 
511
 
                /* look in the range row1 ... row2 */
512
 
                for (row = row1 ; row <= row2 ; row++)
513
 
                {
514
 
                    DEBUG3 (("   Empty Row: "ID"\n", row)) ;
515
 
                    if (NON_PIVOTAL_ROW (row))
516
 
                    {
517
 
                        /* found it */
518
 
                        DEBUG3 (("   Empty Row: "ID" found\n", row)) ;
519
 
                        ASSERT (Frpos [row] == EMPTY) ;
520
 
                        pivrow [OUT] = row ;
521
 
                        DEBUGm4 (("got dummy row pivrow %d\n", pivrow [OUT])) ;
522
 
                        break ;
523
 
                    }
524
 
                }
525
 
                Front_new1strow [nfr] = row ;
526
 
            }
527
 
 
528
 
            if (pivrow [OUT] == EMPTY)
529
 
            {
530
 
                /* Row-merge set is empty.  We can just discard */
531
 
                /* the candidate pivot column. */
532
 
                DEBUG0 (("Note: row-merge set empty\n")) ;
533
 
                DEBUGm4 (("got no pivrow \n")) ;
534
 
                return (UMFPACK_WARNING_singular_matrix) ;
535
 
            }
536
 
        }
537
 
    }
538
 
 
539
 
    /* ---------------------------------------------------------------------- */
540
 
    /* construct the candidate row in the front, if any */
541
 
    /* ---------------------------------------------------------------------- */
542
 
 
543
 
#ifndef NDEBUG
544
 
    /* check Wrp */
545
 
    ASSERT (Work->Wrpflag > 0) ;
546
 
    if (UMF_debug > 0 || Work->n_col < 1000)
547
 
    {
548
 
        for (i = 0 ; i < Work->n_col ; i++)
549
 
        {
550
 
            ASSERT (Wrp [i] < Work->Wrpflag) ;
551
 
        }
552
 
    }
553
 
#endif
554
 
 
555
 
#ifndef NDEBUG
556
 
    DEBUG4 (("pivrow [IN]: "ID"\n", pivrow [IN])) ;
557
 
    UMF_dump_rowcol (0, Numeric, Work, pivrow [IN], TRUE) ;
558
 
#endif
559
 
 
560
 
    if (pivrow [IN] != EMPTY)
561
 
    {
562
 
 
563
 
        /* the row merge candidate row is not pivrow [IN] */
564
 
        freebie [IN] = (pivrow [IN] == prior_pivrow [IN]) && (cdeg1  > 0) ;
565
 
        ASSERT (cdeg1  >= 0) ;
566
 
 
567
 
        if (!freebie [IN])
568
 
        {
569
 
            /* include current front in the degree of this row */
570
 
 
571
 
            Fcpos = Work->Fcpos ;
572
 
            fncols = Work->fncols ;
573
 
 
574
 
            Wrpflag = Work->Wrpflag ;
575
 
 
576
 
            /* -------------------------------------------------------------- */
577
 
            /* construct the pattern of the IN row */
578
 
            /* -------------------------------------------------------------- */
579
 
 
580
 
#ifndef NDEBUG
581
 
            /* check Fcols */
582
 
            DEBUG5 (("ROW ASSEMBLE: rdeg "ID"\nREDUCE ROW "ID"\n",
583
 
                fncols, pivrow [IN])) ;
584
 
            for (j = 0 ; j < fncols ; j++)
585
 
            {
586
 
                col = Work->Fcols [j] ;
587
 
                ASSERT (col >= 0 && col < Work->n_col) ;
588
 
                ASSERT (Fcpos [col] >= 0) ;
589
 
            }
590
 
            if (UMF_debug > 0 || Work->n_col < 1000)
591
 
            {
592
 
                Int cnt = fncols ;
593
 
                for (col = 0 ; col < Work->n_col ; col++)
594
 
                {
595
 
                    if (Fcpos [col] < 0) cnt++ ;
596
 
                }
597
 
                ASSERT (cnt == Work->n_col) ;
598
 
            }
599
 
#endif
600
 
 
601
 
            rdeg [IN] = fncols ;
602
 
 
603
 
            ASSERT (pivrow [IN] >= 0 && pivrow [IN] < Work->n_row) ;
604
 
            ASSERT (NON_PIVOTAL_ROW (pivrow [IN])) ;
605
 
 
606
 
            /* add the pivot column itself */
607
 
            ASSERT (Wrp [pivcol] != Wrpflag) ;
608
 
            if (Fcpos [pivcol] < 0)
609
 
            {
610
 
                DEBUG3 (("Adding pivot col to pivrow [IN] pattern\n")) ;
611
 
                if (rdeg [IN] >= max_rdeg)
612
 
                {
613
 
                    /* :: pattern change (in) :: */
614
 
                    return (UMFPACK_ERROR_different_pattern) ;
615
 
                }
616
 
                Wrp [pivcol] = Wrpflag ;
617
 
                W_i [rdeg [IN]++] = pivcol ;
618
 
            }
619
 
 
620
 
            tpi = Row_tuples [pivrow [IN]] ;
621
 
            if (tpi)
622
 
            {
623
 
                tp = (Tuple *) (Memory + tpi) ;
624
 
                tp1 = tp ;
625
 
                tp2 = tp ;
626
 
                tpend = tp + Row_tlen [pivrow [IN]] ;
627
 
                for ( ; tp < tpend ; tp++)
628
 
                {
629
 
                    e = tp->e ;
630
 
                    ASSERT (e > 0 && e <= Work->nel) ;
631
 
                    if (!E [e])
632
 
                    {
633
 
                        continue ;      /* element already deallocated */
634
 
                    }
635
 
                    f = tp->f ;
636
 
                    p = Memory + E [e] ;
637
 
                    ep = (Element *) p ;
638
 
                    p += UNITS (Element, 1) ;
639
 
                    Cols = (Int *) p ;
640
 
                    ncols = ep->ncols ;
641
 
                    Rows = Cols + ncols ;
642
 
                    if (Rows [f] == EMPTY)
643
 
                    {
644
 
                        continue ;      /* row already assembled */
645
 
                    }
646
 
                    ASSERT (pivrow [IN] == Rows [f]) ;
647
 
 
648
 
                    for (j = 0 ; j < ncols ; j++)
649
 
                    {
650
 
                        col = Cols [j] ;
651
 
                        ASSERT (col >= EMPTY && col < Work->n_col) ;
652
 
                        if ((col >= 0) && (Wrp [col] != Wrpflag)
653
 
                            && Fcpos [col] <0)
654
 
                        {
655
 
                            ASSERT (NON_PIVOTAL_COL (col)) ;
656
 
                            if (rdeg [IN] >= max_rdeg)
657
 
                            {
658
 
                                /* :: pattern change (rdeg in failure) :: */
659
 
                                DEBUGm4 (("rdeg [IN] >= max_rdeg failure\n")) ;
660
 
                                return (UMFPACK_ERROR_different_pattern) ;
661
 
                            }
662
 
                            Wrp [col] = Wrpflag ;
663
 
                            W_i [rdeg [IN]++] = col ;
664
 
                        }
665
 
                    }
666
 
 
667
 
                    *tp2++ = *tp ;      /* leave the tuple in the list */
668
 
                }
669
 
                Row_tlen [pivrow [IN]] = tp2 - tp1 ;
670
 
            }
671
 
 
672
 
#ifndef NDEBUG
673
 
            DEBUG4 (("Reduced IN row:\n")) ;
674
 
            for (j = 0 ; j < fncols ; j++)
675
 
            {
676
 
                DEBUG6 ((" "ID" "ID" "ID"\n",
677
 
                    j, Work->Fcols [j], Fcpos [Work->Fcols [j]])) ;
678
 
                ASSERT (Fcpos [Work->Fcols [j]] >= 0) ;
679
 
            }
680
 
            for (j = fncols ; j < rdeg [IN] ; j++)
681
 
            {
682
 
                DEBUG6 ((" "ID" "ID" "ID"\n", j, W_i [j], Wrp [W_i [j]]));
683
 
                ASSERT (W_i [j] >= 0 && W_i [j] < Work->n_col) ;
684
 
                ASSERT (Wrp [W_i [j]] == Wrpflag) ;
685
 
            }
686
 
            /* mark the end of the pattern in case we scan it by mistake */
687
 
            /* Note that this means W_i must be of size >= fncols_max + 1 */
688
 
            W_i [rdeg [IN]] = EMPTY ;
689
 
#endif
690
 
 
691
 
            /* rdeg [IN] is now the exact degree of the IN row */
692
 
 
693
 
            /* clear Work->Wrp. */
694
 
            Work->Wrpflag++ ;
695
 
            /* All Wrp [0..n_col] is now < Wrpflag */
696
 
        }
697
 
    }
698
 
 
699
 
#ifndef NDEBUG
700
 
    /* check Wrp */
701
 
    if (UMF_debug > 0 || Work->n_col < 1000)
702
 
    {
703
 
        for (i = 0 ; i < Work->n_col ; i++)
704
 
        {
705
 
            ASSERT (Wrp [i] < Work->Wrpflag) ;
706
 
        }
707
 
    }
708
 
#endif
709
 
 
710
 
    /* ---------------------------------------------------------------------- */
711
 
    /* construct the candidate row not in the front, if any */
712
 
    /* ---------------------------------------------------------------------- */
713
 
 
714
 
#ifndef NDEBUG
715
 
    DEBUG4 (("pivrow [OUT]: "ID"\n", pivrow [OUT])) ;
716
 
    UMF_dump_rowcol (0, Numeric, Work, pivrow [OUT], TRUE) ;
717
 
#endif
718
 
 
719
 
    /* If this is a candidate row from the row merge set, force it to be */
720
 
    /* scanned (ignore prior_pivrow [OUT]). */
721
 
 
722
 
    if (pivrow [OUT] != EMPTY)
723
 
    {
724
 
        freebie [OUT] = (pivrow [OUT] == prior_pivrow [OUT]) && cdeg1  > 0 ;
725
 
        ASSERT (cdeg1  >= 0) ;
726
 
 
727
 
        if (!freebie [OUT])
728
 
        {
729
 
 
730
 
            Wrpflag = Work->Wrpflag ;
731
 
 
732
 
            /* -------------------------------------------------------------- */
733
 
            /* construct the pattern of the row */
734
 
            /* -------------------------------------------------------------- */
735
 
 
736
 
            rdeg [OUT] = 0 ;
737
 
 
738
 
            ASSERT (pivrow [OUT] >= 0 && pivrow [OUT] < Work->n_row) ;
739
 
            ASSERT (NON_PIVOTAL_ROW (pivrow [OUT])) ;
740
 
 
741
 
            /* add the pivot column itself */
742
 
            ASSERT (Wrp [pivcol] != Wrpflag) ;
743
 
            DEBUG3 (("Adding pivot col to pivrow [OUT] pattern\n")) ;
744
 
            if (rdeg [OUT] >= max_rdeg)
745
 
            {
746
 
                /* :: pattern change (out) :: */
747
 
                return (UMFPACK_ERROR_different_pattern) ;
748
 
            }
749
 
            Wrp [pivcol] = Wrpflag ;
750
 
            W_o [rdeg [OUT]++] = pivcol ;
751
 
 
752
 
            tpi = Row_tuples [pivrow [OUT]] ;
753
 
            if (tpi)
754
 
            {
755
 
                tp = (Tuple *) (Memory + tpi) ;
756
 
                tp1 = tp ;
757
 
                tp2 = tp ;
758
 
                tpend = tp + Row_tlen [pivrow [OUT]] ;
759
 
                for ( ; tp < tpend ; tp++)
760
 
                {
761
 
                    e = tp->e ;
762
 
                    ASSERT (e > 0 && e <= Work->nel) ;
763
 
                    if (!E [e])
764
 
                    {
765
 
                        continue ;      /* element already deallocated */
766
 
                    }
767
 
                    f = tp->f ;
768
 
                    p = Memory + E [e] ;
769
 
                    ep = (Element *) p ;
770
 
                    p += UNITS (Element, 1) ;
771
 
                    Cols = (Int *) p ;
772
 
                    ncols = ep->ncols ;
773
 
                    Rows = Cols + ncols ;
774
 
                    if (Rows [f] == EMPTY)
775
 
                    {
776
 
                        continue ;      /* row already assembled */
777
 
                    }
778
 
                    ASSERT (pivrow [OUT] == Rows [f]) ;
779
 
 
780
 
                    for (j = 0 ; j < ncols ; j++)
781
 
                    {
782
 
                        col = Cols [j] ;
783
 
                        ASSERT (col >= EMPTY && col < Work->n_col) ;
784
 
                        if ((col >= 0) && (Wrp [col] != Wrpflag))
785
 
                        {
786
 
                            ASSERT (NON_PIVOTAL_COL (col)) ;
787
 
                            if (rdeg [OUT] >= max_rdeg)
788
 
                            {
789
 
                                /* :: pattern change (rdeg out failure) :: */
790
 
                                DEBUGm4 (("rdeg [OUT] failure\n")) ;
791
 
                                return (UMFPACK_ERROR_different_pattern) ;
792
 
                            }
793
 
                            Wrp [col] = Wrpflag ;
794
 
                            W_o [rdeg [OUT]++] = col ;
795
 
                        }
796
 
                    }
797
 
                    *tp2++ = *tp ;      /* leave the tuple in the list */
798
 
                }
799
 
                Row_tlen [pivrow [OUT]] = tp2 - tp1 ;
800
 
            }
801
 
 
802
 
#ifndef NDEBUG
803
 
            DEBUG4 (("Reduced row OUT:\n")) ;
804
 
            for (j = 0 ; j < rdeg [OUT] ; j++)
805
 
            {
806
 
                DEBUG6 ((" "ID" "ID" "ID"\n", j, W_o [j], Wrp [W_o [j]])) ;
807
 
                ASSERT (W_o [j] >= 0 && W_o [j] < Work->n_col) ;
808
 
                ASSERT (Wrp [W_o [j]] == Wrpflag) ;
809
 
            }
810
 
            /* mark the end of the pattern in case we scan it by mistake */
811
 
            /* Note that this means W_o must be of size >= fncols_max + 1 */
812
 
            W_o [rdeg [OUT]] = EMPTY ;
813
 
#endif
814
 
 
815
 
            /* rdeg [OUT] is now the exact degree of the row */
816
 
 
817
 
            /* clear Work->Wrp. */
818
 
            Work->Wrpflag++ ;
819
 
            /* All Wrp [0..n] is now < Wrpflag */
820
 
 
821
 
        }
822
 
 
823
 
    }
824
 
    DEBUG5 (("Row_search end ] \n")) ;
825
 
 
826
 
#ifndef NDEBUG
827
 
    /* check Wrp */
828
 
    if (UMF_debug > 0 || Work->n_col < 1000)
829
 
    {
830
 
        for (i = 0 ; i < Work->n_col ; i++)
831
 
        {
832
 
            ASSERT (Wrp [i] < Work->Wrpflag) ;
833
 
        }
834
 
    }
835
 
#endif
836
 
 
837
 
    return (UMFPACK_OK) ;
838
 
}