1 // =============================================================================
2 // === spqr_debug ==============================================================
3 // =============================================================================
4
5 #include "spqr.hpp"
6
7 // =============================================================================
8 // This file is for debugging only.
9 // =============================================================================
10
11
12 #ifndef NDEBUG
13
14 #ifndef NPRINT
15
16 // =============================================================================
17 // === spqrDebug_print =========================================================
18 // =============================================================================
19
spqrDebug_print(double x)20 void spqrDebug_print
21 (
22 double x
23 )
24 {
25 PR ((" %10.4g", x)) ;
26 }
27
spqrDebug_print(Complex x)28 void spqrDebug_print
29 (
30 Complex x
31 )
32 {
33 PR ((" (%10.4g + 1i*(%10.4g))", x.real ( ), x.imag ( ))) ;
34 }
35
36 // =============================================================================
37 // === spqrDebug_dumpdense =====================================================
38 // =============================================================================
39
spqrDebug_dumpdense(Entry * A,Long m,Long n,Long lda,cholmod_common * cc)40 template <typename Entry> void spqrDebug_dumpdense
41 (
42 Entry *A,
43 Long m,
44 Long n,
45 Long lda,
46 cholmod_common *cc
47 )
48 {
49 Long i, j ;
50 if (cc == NULL) return ;
51 PR (("Dense: m %ld n %ld lda %ld p %p\n", m, n, lda, A)) ;
52 if (m < 0 || n < 0 || lda < m || A == NULL)
53 {
54 PR (("bad dense matrix!\n")) ;
55 ASSERT (0) ;
56 return ;
57 }
58
59 for (j = 0 ; j < n ; j++)
60 {
61 PR ((" --- column %ld of %ld\n", j, n)) ;
62 for (i = 0 ; i < m ; i++)
63 {
64 if (i == j) PR ((" [ diag: ")) ;
65 else PR ((" row %4ld ", i)) ;
66 spqrDebug_print (A [i + j*lda]) ;
67 if (i == j) PR ((" ]\n")) ;
68 else PR (("\n")) ;
69 }
70 PR (("\n")) ;
71 }
72
73 #if 0
74 for (i = 0 ; i < m ; i++)
75 {
76 for (j = 0 ; j < n ; j++)
77 {
78 // if (A [i+j*lda] != (Entry) 0)
79 // {
80 // PR (("X")) ;
81 // }
82 // else
83 // {
84 // PR ((".")) ;
85 // }
86 spqrDebug_print (A [i + j*lda]) ;
87 }
88 PR (("\n")) ;
89 }
90 #endif
91
92 }
93
94 template void spqrDebug_dumpdense <double>
95 (
96 double *A,
97 Long m,
98 Long n,
99 Long lda,
100 cholmod_common *cc
101 ) ;
102
103 template void spqrDebug_dumpdense <Complex>
104 (
105 Complex *A,
106 Long m,
107 Long n,
108 Long lda,
109 cholmod_common *cc
110 ) ;
111
112 // =============================================================================
113 // === spqrDebug_dumpsparse ====================================================
114 // =============================================================================
115
spqrDebug_dumpsparse(Long * Ap,Long * Ai,Entry * Ax,Long m,Long n,cholmod_common * cc)116 template <typename Entry> void spqrDebug_dumpsparse
117 (
118 Long *Ap,
119 Long *Ai,
120 Entry *Ax,
121 Long m,
122 Long n,
123 cholmod_common *cc
124 )
125 {
126 Long p, i, j ;
127 if (cc == NULL) return ;
128 PR (("\nSparse: m %ld n %ld nz %ld Ap %p Ai %p Ax %p\n",
129 m, n, Ap [n], Ap, Ai,Ax)) ;
130 if (m < 0 || n < 0 || Ax == NULL || Ap == NULL || Ai == NULL) return ;
131 for (j = 0 ; j < n ; j++)
132 {
133 PR ((" column %ld\n", j)) ;
134 for (p = Ap [j] ; p < Ap [j+1] ; p++)
135 {
136 i = Ai [p] ;
137 PR ((" %ld :", i)) ;
138 spqrDebug_print (Ax [p]) ;
139 PR (("\n")) ;
140 ASSERT (i >= 0 && i < m) ;
141 }
142 }
143 }
144
145 template void spqrDebug_dumpsparse <double>
146 (
147 Long *Ap,
148 Long *Ai,
149 double *Ax,
150 Long m,
151 Long n,
152 cholmod_common *cc
153 ) ;
154
155 template void spqrDebug_dumpsparse <Complex>
156 (
157 Long *Ap,
158 Long *Ai,
159 Complex *Ax,
160 Long m,
161 Long n,
162 cholmod_common *cc
163 ) ;
164
165 #endif
166
167 // =============================================================================
168 // === spqrDebug_listcount =====================================================
169 // =============================================================================
170
171 #ifdef DEBUG_EXPENSIVE
172
173 // returns # of times x is in the List [0..len-1]
spqrDebug_listcount(Long x,Long * List,Long len,Long what,cholmod_common * cc)174 Long spqrDebug_listcount
175 (
176 Long x, Long *List, Long len, Long what,
177 cholmod_common *cc
178 )
179 {
180 Long k, nfound = 0 ;
181 if (cc == NULL) return (EMPTY) ;
182 if (what == 0)
183 {
184 k = 0 ;
185 PR (("\nQfill, j %ld len %ld\n", x, len)) ;
186 }
187 if (what == 1)
188 {
189 k = 0 ;
190 PR (("\nQrows, i %ld len %ld\n", x, len)) ;
191 }
192 for (k = 0 ; k < len ; k++)
193 {
194 if (List [k] == x) nfound++ ;
195 PR ((" %ld ( %ld ) %ld\n", x, List [k], nfound)) ;
196 }
197 PR (("total found %ld\n\n", nfound)) ;
198 return (nfound) ;
199 }
200 #endif
201
202 // =============================================================================
203 // === spqrDebug_rhsize ========================================================
204 // =============================================================================
205
206 // Count the number of entries in the R+H block for a single front.
207
spqrDebug_rhsize(Long m,Long n,Long npiv,Long * Stair,cholmod_common * cc)208 Long spqrDebug_rhsize // returns # of entries in R+H
209 (
210 // input, not modified
211 Long m, // # of rows in F
212 Long n, // # of columns in F
213 Long npiv, // number of pivotal columns in F
214 Long *Stair, // size n; column j is dead if Stair [j] == 0.
215 // Only the first npiv columns can be dead.
216 cholmod_common *cc
217 )
218 {
219 Long k, h, t, rm, rhsize = 0 ;
220
221 ASSERT (m >= 0 && n >= 0 && npiv <= n && npiv >= 0) ;
222
223 if (cc == NULL) return (EMPTY) ;
224 if (m <= 0 || n <= 0) return (0) ; // nothing to do
225
226 PR (("Try RHSIZE: m %ld n %ld npiv %ld\n", m, n, npiv)) ;
227
228 // -------------------------------------------------------------------------
229 // count the squeezed part of R+H
230 // -------------------------------------------------------------------------
231
232 rm = 0 ; // number of rows in R (:,0:k)
233 for (k = 0 ; k < npiv ; k++)
234 {
235 // get the staircase
236 t = Stair [k] ; // F (0:t-1,k) contains R and H
237 if (t == 0)
238 {
239 t = rm ; // dead col, R (0:rm-1,k) only, no H
240 }
241 else if (rm < m)
242 {
243 rm++ ; // col k not dead; one more row of R
244 }
245 PR ((" for RHSIZE, k %ld Stair %ld t %ld (piv)\n", k, Stair[k], t)) ;
246 // pack R (0:rm-1,k) and H (rm:t-1,k)
247 rhsize += t ;
248 }
249
250 // -------------------------------------------------------------------------
251 // count the rectangular part of R and trapezoidal part of H
252 // -------------------------------------------------------------------------
253
254 h = rm ; // the column of H starts in row h
255 for ( ; k < n ; k++)
256 {
257 // get the staircase
258 t = Stair [k] ;
259 // pack R (0:rm-1,k)
260 rhsize += rm ;
261 h = MIN (h+1, m) ; // one more row of C to skip over
262 // pack H (h:t-1,k)
263 PR ((" for RHSIZE, k %ld Stair %ld t %ld\n", k, Stair[k], t)) ;
264 rhsize += (t-h) ;
265 }
266
267 PR ((" RHSIZE: m %ld n %ld npiv %ld is %ld\n", m, n, npiv, rhsize)) ;
268 return (rhsize) ; // return # of entries in R+H
269 }
270
271
272 // =============================================================================
273 // === spqrDebug_dump_Parent ===================================================
274 // =============================================================================
275
spqrDebug_dump_Parent(Long n,Long * Parent,const char * filename)276 void spqrDebug_dump_Parent (Long n, Long *Parent, const char *filename)
277 {
278 FILE *pfile = fopen (filename, "w") ;
279 if (Parent == NULL)
280 {
281 fprintf (pfile, "0\n") ;
282 }
283 else
284 {
285 for (Long f = 0 ; f < n ; f++)
286 {
287 fprintf (pfile, "%ld\n", 1+Parent [f]) ;
288 }
289 }
290 fclose (pfile) ;
291 }
292 #endif
293