1 // =============================================================================
2 // === spqr_hpinv ==============================================================
3 // =============================================================================
4
5 // Finalizes the row permutation that is implicit in the pattern of H. This
6 // must be done sequentially, after all threads have finished factorizing the
7 // matrix and finding all the dead columns. Also determines QRnum->maxfm.
8
9 #include "spqr.hpp"
10
spqr_hpinv(spqr_symbolic * QRsym,spqr_numeric<Entry> * QRnum,Long * W)11 template <typename Entry> void spqr_hpinv
12 (
13 // input
14 spqr_symbolic *QRsym,
15 // input/output
16 spqr_numeric <Entry> *QRnum,
17 // workspace
18 Long *W // size QRnum->m
19 )
20 {
21 Long *Hi, *Hii, *Hip, *HPinv, *Hr, *Super, *Rp, *Hm, *Sleft, *PLinv ;
22 Long nf, m, n, f, rm, i, row1, row2, fm, fn, fp, cm, cn, maxfm ;
23
24 // -------------------------------------------------------------------------
25 // get the contents of the QRsym and QRnum objects
26 // -------------------------------------------------------------------------
27
28 // this function must not be called if Householder vectors weren't kept
29 ASSERT (QRnum->keepH) ;
30
31 nf = QRsym->nf ;
32 m = QRsym->m ;
33 n = QRsym->n ;
34 Hr = QRnum->Hr ;
35 Hm = QRnum->Hm ;
36 Hii = QRnum->Hii ;
37 Hip = QRsym->Hip ;
38 HPinv = QRnum->HPinv ;
39 Super = QRsym->Super ;
40 Rp = QRsym->Rp ;
41 Sleft = QRsym->Sleft ;
42 PLinv = QRsym->PLinv ;
43 maxfm = 0 ;
44
45 #ifndef NDEBUG
46 for (f = 0 ; f < nf ; f++)
47 {
48 Long j ;
49 rm = 0 ;
50 for (j = Super [f] ; j < Super [f+1] ; j++)
51 {
52 if (!(QRnum->Rdead [j]))
53 {
54 rm++ ; // column j is not dead
55 }
56 }
57 ASSERT (Hr [f] == rm) ; // # rows in R block
58 }
59 for (i = 0 ; i < m ; i++)
60 {
61 W [i] = EMPTY ;
62 PR (("For S, PLinv row perm (%ld) = %ld\n", i, PLinv [i])) ;
63 }
64 #endif
65
66 // -------------------------------------------------------------------------
67 // extract the inverse permutation for R1
68 // -------------------------------------------------------------------------
69
70 row1 = 0 ; // number of squeezed rows of R
71 row2 = m ; // for ordering empty rows of R
72
73 // order the empty rows of S
74 ASSERT (Sleft [n+1] == m) ;
75 PR (("Sleft [%ld] = %ld, m = %ld\n", n, Sleft [n], m)) ;
76 for (i = Sleft [n] ; i < m ; i++)
77 {
78 // row i of S is empty, and so it does not appear as a row in the Hi
79 // pattern of any front
80 PR (("empty row %ld\n", i)) ;
81 W [i] = (--row2) ;
82 }
83
84 // order the non-empty rows of S that appear in frontal matrices
85 for (f = 0 ; f < nf ; f++) // traverse in natural order
86 {
87 Hi = &Hii [Hip [f]] ; // list of row indices of H
88
89 // order the pivotal rows of F
90 rm = Hr [f] ; // number of rows in R block
91 for (i = 0 ; i < rm ; i++)
92 {
93 // row1 is a row of R; it is row Hi [i] of H and S
94 ASSERT (Hi [i] >= 0 && Hi [i] < m) ;
95 ASSERT (W [Hi [i]] == EMPTY) ;
96 W [Hi [i]] = row1++ ;
97 }
98
99 // order the non-pivotal rows of F which are not in the C block
100 fp = Super [f+1] - Super [f] ;
101 fn = Rp [f+1] - Rp [f] ;
102 fm = Hm [f] ;
103 maxfm = MAX (maxfm, fm) ;
104 cn = fn - fp ;
105 cm = MIN (fm - rm, cn) ;
106 for (i = fm-1 ; i >= rm + cm ; i--)
107 {
108 // row2 is an empty row of R (does not appear in the rank-size R);
109 // it is row Hi [i] of H and S
110 ASSERT (Hi [i] >= 0 && Hi [i] < m) ;
111 ASSERT (W [Hi [i]] == EMPTY) ;
112 W [Hi [i]] = (--row2) ;
113 }
114 }
115 ASSERT (row1 == QRnum->rank) ;
116 ASSERT (row1 == row2) ;
117
118 QRnum->maxfm = maxfm ;
119
120 #ifndef NDEBUG
121 for (i = 0 ; i < m ; i++)
122 {
123 PR (("H row perm W (%ld) = %ld\n", i, W [i])) ;
124 ASSERT (W [i] >= 0 && W [i] < m) ;
125 }
126 #endif
127
128 // -------------------------------------------------------------------------
129 // combine the permutations
130 // -------------------------------------------------------------------------
131
132 // At this point, W [i] = k if row i of S is row k of R, and
133 // PLinv [i] = k if row i of A is row k of S. Combine the two permutations
134 // into a single one, HPinv [row i of A] = row k of R
135
136 for (i = 0 ; i < m ; i++)
137 {
138 HPinv [i] = W [PLinv [i]] ;
139 PR (("Combined H row perm HPinv (%ld) = %ld\n", i, HPinv [i])) ;
140 }
141
142 // -------------------------------------------------------------------------
143 // revise the pattern of the frontal matrices
144 // -------------------------------------------------------------------------
145
146 for (f = 0 ; f < nf ; f++)
147 {
148 Hi = &Hii [Hip [f]] ; // list of row indices of H
149 fm = Hm [f] ;
150 for (i = 0 ; i < fm ; i++)
151 {
152 ASSERT (Hi [i] >= 0 && Hi [i] < m) ;
153 Hi [i] = W [Hi [i]] ;
154 }
155 // Now Hi [0..fm-1] contains a list of row indices for front F which
156 // correspond to rows of R. The row indices are also sorted.
157 #ifndef NDEBUG
158 for (i = 1 ; i < fm ; i++)
159 {
160 ASSERT (Hi [i] > Hi [i-1]) ;
161 }
162 #endif
163 }
164 }
165
166
167 // =============================================================================
168
169 template void spqr_hpinv <double>
170 (
171 // input
172 spqr_symbolic *QRsym,
173 // input/output
174 spqr_numeric <double> *QRnum,
175 // workspace
176 Long *W // size QRnum->m
177 ) ;
178
179 // =============================================================================
180
181 template void spqr_hpinv <Complex>
182 (
183 // input
184 spqr_symbolic *QRsym,
185 // input/output
186 spqr_numeric <Complex> *QRnum,
187 // workspace
188 Long *W // size QRnum->m
189 ) ;
190