1 // =============================================================================
2 // === spqr_fsize ==============================================================
3 // =============================================================================
4
5 // Compute the number of rows in front F, initialize its staircase, and
6 // create its Fmap
7
8 #include "spqr.hpp"
9
spqr_fsize(Long f,Long * Super,Long * Rp,Long * Rj,Long * Sleft,Long * Child,Long * Childp,Long * Cm,Long * Fmap,Long * Stair)10 Long spqr_fsize // returns # of rows of F
11 (
12 // inputs, not modified
13 Long f,
14 Long *Super, // size nf, from QRsym
15 Long *Rp, // size nf, from QRsym
16 Long *Rj, // size rjsize, from QRsym
17 Long *Sleft, // size n+2, from QRsym
18 Long *Child, // size nf, from QRsym
19 Long *Childp, // size nf+1, from QRsym
20 Long *Cm, // size nf
21
22 // outputs, not defined on input
23 Long *Fmap, // size n
24 Long *Stair // size fn
25 )
26 {
27 Long col1, col2, p1, p2, fp, fn, fm, col, p, j, c, pc, cm, ci, t, fpc ;
28
29 // -------------------------------------------------------------------------
30 // get the front F
31 // -------------------------------------------------------------------------
32
33 // pivotal columns Super [f] ... Super [f+1]-1
34 col1 = Super [f] ; // front F has columns col1:col2-1
35 col2 = Super [f+1] ;
36 p1 = Rp [f] ; // Rj [p1:p2-1] = columns in F
37 p2 = Rp [f+1] ;
38 fp = col2 - col1 ; // first fp columns are pivotal
39 fn = p2 - p1 ; // exact number of columns of F
40 ASSERT (fp > 0) ; // all fronts have at least one pivot column
41 ASSERT (fp <= fn) ;
42
43 PR (("\n---- Front: %ld\n", f)) ;
44 PR (("Get Fsize: col1 %ld col2 %ld fp %ld fn %ld\n", col1, col2, fp, fn)) ;
45
46 // -------------------------------------------------------------------------
47 // create the Fmap for front F
48 // -------------------------------------------------------------------------
49
50 for (p = p1, j = 0 ; p < p2 ; p++, j++)
51 {
52 col = Rj [p] ; // global column col is jth col of F
53 Fmap [col] = j ;
54 PR (("Fmap (%ld): %ld \n", col, j)) ;
55 }
56
57 // -------------------------------------------------------------------------
58 // initialize the staircase for front F
59 // -------------------------------------------------------------------------
60
61 // initialize the staircase with original rows of S
62 for (j = 0 ; j < fp ; j++)
63 {
64 // global column j+col1 is the jth pivot column of front F
65 col = j + col1 ;
66 Stair [j] = Sleft [col+1] - Sleft [col] ;
67 PR (("init rows, j: %ld count %ld\n", j, Stair [j])) ;
68 }
69
70 // contribution blocks from children will be added here
71 for ( ; j < fn ; j++)
72 {
73 Stair [j] = 0 ;
74 }
75
76 // -------------------------------------------------------------------------
77 // construct the staircase for each child
78 // -------------------------------------------------------------------------
79
80 for (p = Childp [f] ; p < Childp [f+1] ; p++)
81 {
82 c = Child [p] ; // get the child c of front F
83 PR (("child %ld\n", c)) ;
84 pc = Rp [c] ; // get the pattern of child R
85 cm = Cm [c] ; // # of rows in child C
86 // fnc = (Rp [c+1] - pc) ; // total # cols in child F
87 fpc = Super [c+1] - Super [c] ; // # of pivot cols in child
88 // cn = (fnc - fpc) ; // # of cols in child C
89 // ASSERT (cm >= 0 && cm <= cn) ;
90 pc += fpc ; // pointer to column indices in C
91 // ASSERT (pc + cn == Rp [c+1]) ;
92 // PR ((" cm %ld cn %ld\n", cm, cn)) ;
93
94 // add the child rows to the staircase
95 for (ci = 0 ; ci < cm ; ci++)
96 {
97 col = Rj [pc + ci] ; // leftmost col of this row of C
98 j = Fmap [col] ; // global col is jth col of F
99 PR ((" child col %ld j %ld\n", col, j)) ;
100 ASSERT (j >= 0 && j < fn) ;
101 Stair [j]++ ; // add this row to jth staircase
102 }
103 }
104
105 // -------------------------------------------------------------------------
106 // replace Stair with cumsum ([0 Stair]), and find # rows of F
107 // -------------------------------------------------------------------------
108
109 fm = 0 ;
110 for (j = 0 ; j < fn ; j++)
111 {
112 t = fm ;
113 fm += Stair [j] ;
114 Stair [j] = t ;
115 }
116 PR (("fm %ld %ld\n", fm, Stair [fn-1])) ;
117
118 // Note that the above code differs from the symbolic analysis. The
119 // nonzero values in column j of F will reside in row 0 of F to row
120 // Stair [j+1]-1 of F. Once the rows of S and the children of F are
121 // assembled, this will change to Stair [j]-1.
122 //
123 // When a row is assembled into F whose leftmost column is j, it is placed
124 // as row Stair [j], which is then incremented to accommodate the next row
125 // with leftmost column j. At that point, Stair [0:fn-1] will then "equal"
126 // Stair [0:fn-1] in the symbolic analysis (except that the latter is an
127 // upper bound if rank detection has found failed pivot columns).
128
129 return (fm) ;
130 }
131