1 /* ========================================================================== */
2 /* === UMF_extend_front ===================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* Copyright (c) 2005-2012 by Timothy A. Davis, http://www.suitesparse.com. */
7 /* All Rights Reserved. See ../Doc/License.txt for License. */
8 /* -------------------------------------------------------------------------- */
9
10 /* Called by kernel. */
11
12 #include "umf_internal.h"
13 #include "umf_extend_front.h"
14 #include "umf_grow_front.h"
15
16 /* ========================================================================== */
17 /* === zero_front =========================================================== */
18 /* ========================================================================== */
19
zero_front(Entry * Flblock,Entry * Fublock,Entry * Fcblock,Int fnrows,Int fncols,Int fnr_curr,Int fnc_curr,Int fnpiv,Int fnrows_extended,Int fncols_extended)20 PRIVATE void zero_front (
21 Entry *Flblock, Entry *Fublock, Entry *Fcblock,
22 Int fnrows, Int fncols, Int fnr_curr, Int fnc_curr,
23 Int fnpiv, Int fnrows_extended, Int fncols_extended)
24 {
25 Int j, i ;
26 Entry *F, *Fj, *Fi ;
27
28 Fj = Fcblock + fnrows ;
29 for (j = 0 ; j < fncols ; j++)
30 {
31 /* zero the new rows in the contribution block: */
32 F = Fj ;
33 Fj += fnr_curr ;
34 #pragma ivdep
35 for (i = fnrows ; i < fnrows_extended ; i++)
36 {
37 /* CLEAR (Fcblock [i + j*fnr_curr]) ; */
38 CLEAR_AND_INCREMENT (F) ;
39 }
40 }
41
42 Fj -= fnrows ;
43 for (j = fncols ; j < fncols_extended ; j++)
44 {
45 /* zero the new columns in the contribution block: */
46 F = Fj ;
47 Fj += fnr_curr ;
48 #pragma ivdep
49 for (i = 0 ; i < fnrows_extended ; i++)
50 {
51 /* CLEAR (Fcblock [i + j*fnr_curr]) ; */
52 CLEAR_AND_INCREMENT (F) ;
53 }
54 }
55
56 Fj = Flblock + fnrows ;
57 for (j = 0 ; j < fnpiv ; j++)
58 {
59 /* zero the new rows in L block: */
60 F = Fj ;
61 Fj += fnr_curr ;
62 #pragma ivdep
63 for (i = fnrows ; i < fnrows_extended ; i++)
64 {
65 /* CLEAR (Flblock [i + j*fnr_curr]) ; */
66 CLEAR_AND_INCREMENT (F) ;
67 }
68 }
69
70 Fi = Fublock + fncols ;
71 for (i = 0 ; i < fnpiv ; i++)
72 {
73 /* zero the new columns in U block: */
74 F = Fi ;
75 Fi += fnc_curr ;
76 #pragma ivdep
77 for (j = fncols ; j < fncols_extended ; j++)
78 {
79 /* CLEAR (Fublock [i*fnc_curr + j]) ; */
80 CLEAR_AND_INCREMENT (F) ;
81 }
82 }
83
84 }
85
86 /* ========================================================================== */
87 /* === UMF_extend_front ===================================================== */
88 /* ========================================================================== */
89
UMF_extend_front(NumericType * Numeric,WorkType * Work)90 GLOBAL Int UMF_extend_front
91 (
92 NumericType *Numeric,
93 WorkType *Work
94 )
95 {
96 /* ---------------------------------------------------------------------- */
97 /* local variables */
98 /* ---------------------------------------------------------------------- */
99
100 Int j, i, *Frows, row, col, *Wrow, fnr2, fnc2, *Frpos, *Fcpos, *Fcols,
101 fnrows_extended, rrdeg, ccdeg, fncols_extended, fnr_curr, fnc_curr,
102 fnrows, fncols, pos, fnpiv, *Wm ;
103 Entry *Wx, *Wy, *Fu, *Fl ;
104
105 /* ---------------------------------------------------------------------- */
106 /* get current frontal matrix and check for frontal growth */
107 /* ---------------------------------------------------------------------- */
108
109 fnpiv = Work->fnpiv ;
110
111 #ifndef NDEBUG
112 DEBUG2 (("EXTEND FRONT\n")) ;
113 DEBUG2 (("Work->fnpiv "ID"\n", fnpiv)) ;
114 ASSERT (Work->Flblock == Work->Flublock + Work->nb*Work->nb) ;
115 ASSERT (Work->Fublock == Work->Flblock + Work->fnr_curr*Work->nb) ;
116 ASSERT (Work->Fcblock == Work->Fublock + Work->nb*Work->fnc_curr) ;
117 DEBUG7 (("C block: ")) ;
118 UMF_dump_dense (Work->Fcblock, Work->fnr_curr, Work->fnrows, Work->fncols) ;
119 DEBUG7 (("L block: ")) ;
120 UMF_dump_dense (Work->Flblock, Work->fnr_curr, Work->fnrows, fnpiv);
121 DEBUG7 (("U' block: ")) ;
122 UMF_dump_dense (Work->Fublock, Work->fnc_curr, Work->fncols, fnpiv) ;
123 DEBUG7 (("LU block: ")) ;
124 UMF_dump_dense (Work->Flublock, Work->nb, fnpiv, fnpiv) ;
125 #endif
126
127 if (Work->do_grow)
128 {
129 fnr2 = UMF_FRONTAL_GROWTH * Work->fnrows_new + 2 ;
130 fnc2 = UMF_FRONTAL_GROWTH * Work->fncols_new + 2 ;
131 if (!UMF_grow_front (Numeric, fnr2, fnc2, Work, 1))
132 {
133 DEBUGm4 (("out of memory: extend front\n")) ;
134 return (FALSE) ;
135 }
136 }
137
138 fnr_curr = Work->fnr_curr ;
139 fnc_curr = Work->fnc_curr ;
140 ASSERT (Work->fnrows_new + 1 <= fnr_curr) ;
141 ASSERT (Work->fncols_new + 1 <= fnc_curr) ;
142 ASSERT (fnr_curr >= 0 && fnr_curr % 2 == 1) ;
143
144 /* ---------------------------------------------------------------------- */
145 /* get parameters */
146 /* ---------------------------------------------------------------------- */
147
148 Frows = Work->Frows ;
149 Frpos = Work->Frpos ;
150 Fcols = Work->Fcols ;
151 Fcpos = Work->Fcpos ;
152 fnrows = Work->fnrows ;
153 fncols = Work->fncols ;
154 rrdeg = Work->rrdeg ;
155 ccdeg = Work->ccdeg ;
156
157 /* scan starts at the first new column in Fcols */
158 /* also scan the pivot column if it was not in the front */
159 Work->fscan_col = fncols ;
160 Work->NewCols = Fcols ;
161
162 /* scan1 starts at the first new row in Frows */
163 /* also scan the pivot row if it was not in the front */
164 Work->fscan_row = fnrows ;
165 Work->NewRows = Frows ;
166
167 /* ---------------------------------------------------------------------- */
168 /* extend row pattern of the front with the new pivot column */
169 /* ---------------------------------------------------------------------- */
170
171 fnrows_extended = fnrows ;
172 fncols_extended = fncols ;
173
174 #ifndef NDEBUG
175 DEBUG2 (("Pivot col, before extension: "ID"\n", fnrows)) ;
176 for (i = 0 ; i < fnrows ; i++)
177 {
178 DEBUG2 ((" "ID": row "ID"\n", i, Frows [i])) ;
179 ASSERT (Frpos [Frows [i]] == i) ;
180 }
181 DEBUG2 (("Extending pivot column: pivcol_in_front: "ID"\n",
182 Work->pivcol_in_front)) ;
183 #endif
184
185 Fl = Work->Flblock + fnpiv * fnr_curr ;
186
187 if (Work->pivcol_in_front)
188 {
189 /* extended pattern and position already in Frows, Frpos. Values above
190 * the diagonal are already in LU block. Values on and below the
191 * diagonal are in Wy [0 .. fnrows_extended-1]. Copy into the L
192 * block. */
193 fnrows_extended += ccdeg ;
194 Wy = Work->Wy ;
195
196 for (i = 0 ; i < fnrows_extended ; i++)
197 {
198 Fl [i] = Wy [i] ;
199 #ifndef NDEBUG
200 row = Frows [i] ;
201 DEBUG2 ((" "ID": row "ID" ", i, row)) ;
202 EDEBUG2 (Fl [i]) ;
203 if (row == Work->pivrow) DEBUG2 ((" <- pivrow")) ;
204 DEBUG2 (("\n")) ;
205 if (i == fnrows - 1) DEBUG2 ((" :::::::\n")) ;
206 ASSERT (row >= 0 && row < Work->n_row) ;
207 ASSERT (Frpos [row] == i) ;
208 #endif
209 }
210
211 }
212 else
213 {
214 /* extended pattern,values is in (Wm,Wx), not yet in the front */
215 Entry *F ;
216 Fu = Work->Flublock + fnpiv * Work->nb ;
217 Wm = Work->Wm ;
218 Wx = Work->Wx ;
219 F = Fu ;
220 for (i = 0 ; i < fnpiv ; i++)
221 {
222 CLEAR_AND_INCREMENT (F) ;
223 }
224 F = Fl ;
225 for (i = 0 ; i < fnrows ; i++)
226 {
227 CLEAR_AND_INCREMENT (F) ;
228 }
229 for (i = 0 ; i < ccdeg ; i++)
230 {
231 row = Wm [i] ;
232 #ifndef NDEBUG
233 DEBUG2 ((" "ID": row "ID" (ext) ", fnrows_extended, row)) ;
234 EDEBUG2 (Wx [i]) ;
235 if (row == Work->pivrow) DEBUG2 ((" <- pivrow")) ;
236 DEBUG2 (("\n")) ;
237 ASSERT (row >= 0 && row < Work->n_row) ;
238 #endif
239 pos = Frpos [row] ;
240 if (pos < 0)
241 {
242 pos = fnrows_extended++ ;
243 Frows [pos] = row ;
244 Frpos [row] = pos ;
245 }
246 Fl [pos] = Wx [i] ;
247 }
248 }
249
250 ASSERT (fnrows_extended <= fnr_curr) ;
251
252 /* ---------------------------------------------------------------------- */
253 /* extend the column pattern of the front with the new pivot row */
254 /* ---------------------------------------------------------------------- */
255
256 #ifndef NDEBUG
257 DEBUG6 (("Pivot row, before extension: "ID"\n", fncols)) ;
258 for (j = 0 ; j < fncols ; j++)
259 {
260 DEBUG7 ((" "ID": col "ID"\n", j, Fcols [j])) ;
261 ASSERT (Fcpos [Fcols [j]] == j * fnr_curr) ;
262 }
263 DEBUG6 (("Extending pivot row:\n")) ;
264 #endif
265
266 if (Work->pivrow_in_front)
267 {
268 if (Work->pivcol_in_front)
269 {
270 ASSERT (Fcols == Work->Wrow) ;
271 for (j = fncols ; j < rrdeg ; j++)
272 {
273 #ifndef NDEBUG
274 col = Fcols [j] ;
275 DEBUG2 ((" "ID": col "ID" (ext)\n", j, col)) ;
276 ASSERT (col != Work->pivcol) ;
277 ASSERT (col >= 0 && col < Work->n_col) ;
278 ASSERT (Fcpos [col] < 0) ;
279 #endif
280 Fcpos [Fcols [j]] = j * fnr_curr ;
281 }
282 }
283 else
284 {
285 /* OUT-IN option: pivcol not in front, but pivrow is in front */
286 Wrow = Work->Wrow ;
287 ASSERT (IMPLIES (Work->pivcol_in_front, Wrow == Fcols)) ;
288 if (Wrow == Fcols)
289 {
290 /* Wrow and Fcols are equivalenced */
291 for (j = fncols ; j < rrdeg ; j++)
292 {
293 col = Wrow [j] ;
294 DEBUG2 ((" "ID": col "ID" (ext)\n", j, col)) ;
295 ASSERT (Fcpos [col] < 0) ;
296 /* Fcols [j] = col ; not needed */
297 Fcpos [col] = j * fnr_curr ;
298 }
299 }
300 else
301 {
302 for (j = fncols ; j < rrdeg ; j++)
303 {
304 col = Wrow [j] ;
305 DEBUG2 ((" "ID": col "ID" (ext)\n", j, col)) ;
306 ASSERT (Fcpos [col] < 0) ;
307 Fcols [j] = col ;
308 Fcpos [col] = j * fnr_curr ;
309 }
310 }
311 }
312 fncols_extended = rrdeg ;
313 }
314 else
315 {
316 ASSERT (Fcols != Work->Wrow) ;
317 Wrow = Work->Wrow ;
318 for (j = 0 ; j < rrdeg ; j++)
319 {
320 col = Wrow [j] ;
321 ASSERT (col >= 0 && col < Work->n_col) ;
322 if (Fcpos [col] < 0)
323 {
324 DEBUG2 ((" col:: "ID" (ext)\n", col)) ;
325 Fcols [fncols_extended] = col ;
326 Fcpos [col] = fncols_extended * fnr_curr ;
327 fncols_extended++ ;
328 }
329 }
330 }
331
332 /* ---------------------------------------------------------------------- */
333 /* pivot row and column have been extended */
334 /* ---------------------------------------------------------------------- */
335
336 #ifndef NDEBUG
337 ASSERT (fncols_extended <= fnc_curr) ;
338 ASSERT (fnrows_extended <= fnr_curr) ;
339
340 DEBUG6 (("Pivot col, after ext: "ID" "ID"\n", fnrows,fnrows_extended)) ;
341 for (i = 0 ; i < fnrows_extended ; i++)
342 {
343 row = Frows [i] ;
344 DEBUG7 ((" "ID": row "ID" pos "ID" old: %d", i, row, Frpos [row],
345 i < fnrows)) ;
346 if (row == Work->pivrow ) DEBUG7 ((" <-- pivrow")) ;
347 DEBUG7 (("\n")) ;
348 ASSERT (Frpos [Frows [i]] == i) ;
349 }
350
351 DEBUG6 (("Pivot row position: "ID"\n", Frpos [Work->pivrow])) ;
352 ASSERT (Frpos [Work->pivrow] >= 0) ;
353 ASSERT (Frpos [Work->pivrow] < fnrows_extended) ;
354
355 DEBUG6 (("Pivot row, after ext: "ID" "ID"\n", fncols,fncols_extended)) ;
356 for (j = 0 ; j < fncols_extended ; j++)
357 {
358 col = Fcols [j] ;
359 DEBUG7 ((" "ID": col "ID" pos "ID" old: %d", j, col, Fcpos [col],
360 j < fncols)) ;
361 if (col == Work->pivcol ) DEBUG7 ((" <-- pivcol")) ;
362 DEBUG7 (("\n")) ;
363 ASSERT (Fcpos [Fcols [j]] == j * fnr_curr) ;
364 }
365
366 DEBUG6 (("Pivot col position: "ID"\n", Fcpos [Work->pivcol])) ;
367 ASSERT (Fcpos [Work->pivcol] >= 0) ;
368 ASSERT (Fcpos [Work->pivcol] < fncols_extended * fnr_curr) ;
369
370 #endif
371
372 /* ---------------------------------------------------------------------- */
373 /* Zero the newly extended frontal matrix */
374 /* ---------------------------------------------------------------------- */
375
376 zero_front (Work->Flblock, Work->Fublock, Work->Fcblock,
377 fnrows, fncols, fnr_curr, fnc_curr,
378 fnpiv, fnrows_extended, fncols_extended) ;
379
380 /* ---------------------------------------------------------------------- */
381 /* finalize extended row and column pattern of the frontal matrix */
382 /* ---------------------------------------------------------------------- */
383
384 Work->fnrows = fnrows_extended ;
385 Work->fncols = fncols_extended ;
386
387 ASSERT (fnrows_extended == Work->fnrows_new + 1) ;
388 ASSERT (fncols_extended == Work->fncols_new + 1) ;
389
390 return (TRUE) ;
391
392 }
393