1 /* ========================================================================== */
2 /* === UMF_scale_column ===================================================== */
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 /*
11 Scale the current pivot column, move the pivot row and column into place,
12 and log the permutation.
13 */
14
15 #include "umf_internal.h"
16 #include "umf_scale_column.h"
17 #include "umf_mem_free_tail_block.h"
18 #include "umf_scale.h"
19
20 /* ========================================================================== */
21 /* === shift_pivot_row ====================================================== */
22 /* ========================================================================== */
23
24 /* Except for the BLAS, most of the time is typically spent in the following
25 * shift_pivot_row routine. It copies the pivot row into the U block, and
26 * then fills in the whole in the C block by shifting the last row of C into
27 * the row vacated by the pivot row.
28 */
29
shift_pivot_row(Entry * Fd,Entry * Fs,Entry * Fe,Int len,Int d)30 PRIVATE void shift_pivot_row (Entry *Fd, Entry *Fs, Entry *Fe, Int len, Int d)
31 {
32 Int j ;
33 for (j = 0 ; j < len ; j++)
34 {
35 Fd [j] = Fs [j*d] ;
36 Fs [j*d] = Fe [j*d] ;
37 }
38 }
39
40 /* ========================================================================== */
41 /* === UMF_scale_column ===================================================== */
42 /* ========================================================================== */
43
UMF_scale_column(NumericType * Numeric,WorkType * Work)44 GLOBAL void UMF_scale_column
45 (
46 NumericType *Numeric,
47 WorkType *Work
48 )
49 {
50 /* ---------------------------------------------------------------------- */
51 /* local variables */
52 /* ---------------------------------------------------------------------- */
53
54 Entry pivot_value ;
55 Entry *Fcol, *Flublock, *Flblock, *Fublock, *Fcblock ;
56 Int k, k1, fnr_curr, fnrows, fncols, *Frpos, *Fcpos, pivrow, pivcol,
57 *Frows, *Fcols, fnc_curr, fnpiv, *Row_tuples, nb,
58 *Col_tuples, *Rperm, *Cperm, fspos, col2, row2 ;
59 #ifndef NDEBUG
60 Int *Col_degree, *Row_degree ;
61 #endif
62
63 /* ---------------------------------------------------------------------- */
64 /* get parameters */
65 /* ---------------------------------------------------------------------- */
66
67 fnrows = Work->fnrows ;
68 fncols = Work->fncols ;
69 fnpiv = Work->fnpiv ;
70
71 /* ---------------------------------------------------------------------- */
72
73 Rperm = Numeric->Rperm ;
74 Cperm = Numeric->Cperm ;
75
76 /* ---------------------------------------------------------------------- */
77
78 Flublock = Work->Flublock ;
79 Flblock = Work->Flblock ;
80 Fublock = Work->Fublock ;
81 Fcblock = Work->Fcblock ;
82
83 fnr_curr = Work->fnr_curr ;
84 fnc_curr = Work->fnc_curr ;
85 Frpos = Work->Frpos ;
86 Fcpos = Work->Fcpos ;
87 Frows = Work->Frows ;
88 Fcols = Work->Fcols ;
89 pivrow = Work->pivrow ;
90 pivcol = Work->pivcol ;
91
92 ASSERT (pivrow >= 0 && pivrow < Work->n_row) ;
93 ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;
94
95 #ifndef NDEBUG
96 Col_degree = Numeric->Cperm ; /* for NON_PIVOTAL_COL macro */
97 Row_degree = Numeric->Rperm ; /* for NON_PIVOTAL_ROW macro */
98 #endif
99
100 Row_tuples = Numeric->Uip ;
101 Col_tuples = Numeric->Lip ;
102 nb = Work->nb ;
103
104 #ifndef NDEBUG
105 ASSERT (fnrows == Work->fnrows_new + 1) ;
106 ASSERT (fncols == Work->fncols_new + 1) ;
107 DEBUG1 (("SCALE COL: fnrows "ID" fncols "ID"\n", fnrows, fncols)) ;
108 DEBUG2 (("\nFrontal matrix, including all space:\n"
109 "fnr_curr "ID" fnc_curr "ID" nb "ID"\n"
110 "fnrows "ID" fncols "ID" fnpiv "ID"\n",
111 fnr_curr, fnc_curr, nb, fnrows, fncols, fnpiv)) ;
112 DEBUG2 (("\nJust the active part:\n")) ;
113 DEBUG7 (("C block: ")) ;
114 UMF_dump_dense (Fcblock, fnr_curr, fnrows, fncols) ;
115 DEBUG7 (("L block: ")) ;
116 UMF_dump_dense (Flblock, fnr_curr, fnrows, fnpiv);
117 DEBUG7 (("U' block: ")) ;
118 UMF_dump_dense (Fublock, fnc_curr, fncols, fnpiv) ;
119 DEBUG7 (("LU block: ")) ;
120 UMF_dump_dense (Flublock, nb, fnpiv, fnpiv) ;
121 #endif
122
123 /* ====================================================================== */
124 /* === Shift pivot row and column ======================================= */
125 /* ====================================================================== */
126
127 /* ---------------------------------------------------------------------- */
128 /* move pivot column into place */
129 /* ---------------------------------------------------------------------- */
130
131 /* Note that the pivot column is already in place. Just shift the last
132 * column into the position vacated by the pivot column. */
133
134 fspos = Fcpos [pivcol] ;
135
136 /* one less column in the contribution block */
137 fncols = --(Work->fncols) ;
138
139 if (fspos != fncols * fnr_curr)
140 {
141
142 Int fs = fspos / fnr_curr ;
143
144 DEBUG6 (("Shift pivot column in front\n")) ;
145 DEBUG6 (("fspos: "ID" flpos: "ID"\n", fspos, fncols * fnr_curr)) ;
146
147 /* ------------------------------------------------------------------ */
148 /* move Fe => Fs */
149 /* ------------------------------------------------------------------ */
150
151 /* column of the contribution block: */
152 {
153 /* Fs: current position of pivot column in contribution block */
154 /* Fe: position of last column in contribution block */
155 Int i ;
156 Entry *Fs, *Fe ;
157 Fs = Fcblock + fspos ;
158 Fe = Fcblock + fncols * fnr_curr ;
159 #pragma ivdep
160 for (i = 0 ; i < fnrows ; i++)
161 {
162 Fs [i] = Fe [i] ;
163 }
164 }
165
166 /* column of the U2 block */
167 {
168 /* Fs: current position of pivot column in U block */
169 /* Fe: last column in U block */
170 Int i ;
171 Entry *Fs, *Fe ;
172 Fs = Fublock + fs ;
173 Fe = Fublock + fncols ;
174 #pragma ivdep
175 for (i = 0 ; i < fnpiv ; i++)
176 {
177 Fs [i * fnc_curr] = Fe [i * fnc_curr] ;
178 }
179 }
180
181 /* move column Fe to Fs in the Fcols pattern */
182 col2 = Fcols [fncols] ;
183 Fcols [fs] = col2 ;
184 Fcpos [col2] = fspos ;
185 }
186
187 /* pivot column is no longer in the frontal matrix */
188 Fcpos [pivcol] = EMPTY ;
189
190 #ifndef NDEBUG
191 DEBUG2 (("\nFrontal matrix after col swap, including all space:\n"
192 "fnr_curr "ID" fnc_curr "ID" nb "ID"\n"
193 "fnrows "ID" fncols "ID" fnpiv "ID"\n",
194 fnr_curr, fnc_curr, nb,
195 fnrows, fncols, fnpiv)) ;
196 DEBUG2 (("\nJust the active part:\n")) ;
197 DEBUG7 (("C block: ")) ;
198 UMF_dump_dense (Fcblock, fnr_curr, fnrows, fncols) ;
199 DEBUG7 (("L block: ")) ;
200 UMF_dump_dense (Flblock, fnr_curr, fnrows, fnpiv+1);
201 DEBUG7 (("U' block: ")) ;
202 UMF_dump_dense (Fublock, fnc_curr, fncols, fnpiv) ;
203 DEBUG7 (("LU block: ")) ;
204 UMF_dump_dense (Flublock, nb, fnpiv, fnpiv+1) ;
205 #endif
206
207 /* ---------------------------------------------------------------------- */
208 /* move pivot row into place */
209 /* ---------------------------------------------------------------------- */
210
211 fspos = Frpos [pivrow] ;
212
213 /* one less row in the contribution block */
214 fnrows = --(Work->fnrows) ;
215
216 DEBUG6 (("Swap/shift pivot row in front:\n")) ;
217 DEBUG6 (("fspos: "ID" flpos: "ID"\n", fspos, fnrows)) ;
218
219 if (fspos == fnrows)
220 {
221
222 /* ------------------------------------------------------------------ */
223 /* move Fs => Fd */
224 /* ------------------------------------------------------------------ */
225
226 DEBUG6 (("row case 1\n")) ;
227
228 /* row of the contribution block: */
229 {
230 Int j ;
231 Entry *Fd, *Fs ;
232 Fd = Fublock + fnpiv * fnc_curr ;
233 Fs = Fcblock + fspos ;
234 #pragma ivdep
235 for (j = 0 ; j < fncols ; j++)
236 {
237 Fd [j] = Fs [j * fnr_curr] ;
238 }
239 }
240
241 /* row of the L2 block: */
242 if (Work->pivrow_in_front)
243 {
244 Int j ;
245 Entry *Fd, *Fs ;
246 Fd = Flublock + fnpiv ;
247 Fs = Flblock + fspos ;
248 #pragma ivdep
249 for (j = 0 ; j <= fnpiv ; j++)
250 {
251 Fd [j * nb] = Fs [j * fnr_curr] ;
252 }
253 }
254 else
255 {
256 Int j ;
257 Entry *Fd, *Fs ;
258 Fd = Flublock + fnpiv ;
259 Fs = Flblock + fspos ;
260 #pragma ivdep
261 for (j = 0 ; j < fnpiv ; j++)
262 {
263 ASSERT (IS_ZERO (Fs [j * fnr_curr])) ;
264 CLEAR (Fd [j * nb]) ;
265 }
266 Fd [fnpiv * nb] = Fs [fnpiv * fnr_curr] ;
267 }
268 }
269 else
270 {
271
272 /* ------------------------------------------------------------------ */
273 /* move Fs => Fd */
274 /* move Fe => Fs */
275 /* ------------------------------------------------------------------ */
276
277 DEBUG6 (("row case 2\n")) ;
278 /* this is the most common case, by far */
279
280 /* row of the contribution block: */
281 {
282 /* Fd: destination of pivot row on U block */
283 /* Fs: current position of pivot row in contribution block */
284 /* Fe: position of last row in contribution block */
285 Entry *Fd, *Fs, *Fe ;
286 Fd = Fublock + fnpiv * fnc_curr ;
287 Fs = Fcblock + fspos ;
288 Fe = Fcblock + fnrows ;
289 shift_pivot_row (Fd, Fs, Fe, fncols, fnr_curr) ;
290 }
291
292 /* row of the L2 block: */
293 if (Work->pivrow_in_front)
294 {
295 /* Fd: destination of pivot row in LU block */
296 /* Fs: current position of pivot row in L block */
297 /* Fe: last row in L block */
298 Int j ;
299 Entry *Fd, *Fs, *Fe ;
300 Fd = Flublock + fnpiv ;
301 Fs = Flblock + fspos ;
302 Fe = Flblock + fnrows ;
303 #pragma ivdep
304 for (j = 0 ; j <= fnpiv ; j++)
305 {
306 Fd [j * nb] = Fs [j * fnr_curr] ;
307 Fs [j * fnr_curr] = Fe [j * fnr_curr] ;
308 }
309 }
310 else
311 {
312 Int j ;
313 Entry *Fd, *Fs, *Fe ;
314 Fd = Flublock + fnpiv ;
315 Fs = Flblock + fspos ;
316 Fe = Flblock + fnrows ;
317 #pragma ivdep
318 for (j = 0 ; j < fnpiv ; j++)
319 {
320 ASSERT (IS_ZERO (Fs [j * fnr_curr])) ;
321 CLEAR (Fd [j * nb]) ;
322 Fs [j * fnr_curr] = Fe [j * fnr_curr] ;
323 }
324 Fd [fnpiv * nb] = Fs [fnpiv * fnr_curr] ;
325 Fs [fnpiv * fnr_curr] = Fe [fnpiv * fnr_curr] ;
326 }
327
328 /* move row Fe to Fs in the Frows pattern */
329 row2 = Frows [fnrows] ;
330 Frows [fspos] = row2 ;
331 Frpos [row2] = fspos ;
332
333 }
334 /* pivot row is no longer in the frontal matrix */
335 Frpos [pivrow] = EMPTY ;
336
337 #ifndef NDEBUG
338 DEBUG2 (("\nFrontal matrix after row swap, including all space:\n"
339 "fnr_curr "ID" fnc_curr "ID" nb "ID"\n"
340 "fnrows "ID" fncols "ID" fnpiv "ID"\n",
341 Work->fnr_curr, Work->fnc_curr, Work->nb,
342 Work->fnrows, Work->fncols, Work->fnpiv)) ;
343 DEBUG2 (("\nJust the active part:\n")) ;
344 DEBUG7 (("C block: ")) ;
345 UMF_dump_dense (Fcblock, fnr_curr, fnrows, fncols) ;
346 DEBUG7 (("L block: ")) ;
347 UMF_dump_dense (Flblock, fnr_curr, fnrows, fnpiv+1);
348 DEBUG7 (("U' block: ")) ;
349 UMF_dump_dense (Fublock, fnc_curr, fncols, fnpiv+1) ;
350 DEBUG7 (("LU block: ")) ;
351 UMF_dump_dense (Flublock, nb, fnpiv+1, fnpiv+1) ;
352 #endif
353
354 /* ---------------------------------------------------------------------- */
355 /* Frpos [row] >= 0 for each row in pivot column pattern. */
356 /* offset into pattern is given by: */
357 /* Frpos [row] == offset - 1 */
358 /* Frpos [pivrow] is EMPTY */
359
360 /* Fcpos [col] >= 0 for each col in pivot row pattern. */
361 /* Fcpos [col] == (offset - 1) * fnr_curr */
362 /* Fcpos [pivcol] is EMPTY */
363
364 /* Fcols [0..fncols-1] is the pivot row pattern (excl pivot cols) */
365 /* Frows [0..fnrows-1] is the pivot col pattern (excl pivot rows) */
366
367 /* ====================================================================== */
368 /* === scale pivot column =============================================== */
369 /* ====================================================================== */
370
371 /* pivot column (except for pivot entry itself) */
372 Fcol = Flblock + fnpiv * fnr_curr ;
373 /* fnpiv-th pivot in frontal matrix located in Flublock (fnpiv, fnpiv) */
374 pivot_value = Flublock [fnpiv + fnpiv * nb] ;
375
376 /* this is the kth global pivot */
377 k = Work->npiv + fnpiv ;
378
379 DEBUG4 (("Pivot value: ")) ;
380 EDEBUG4 (pivot_value) ;
381 DEBUG4 (("\n")) ;
382
383 UMF_scale (fnrows, pivot_value, Fcol) ;
384
385 /* ---------------------------------------------------------------------- */
386 /* deallocate the pivot row and pivot column tuples */
387 /* ---------------------------------------------------------------------- */
388
389 UMF_mem_free_tail_block (Numeric, Row_tuples [pivrow]) ;
390 UMF_mem_free_tail_block (Numeric, Col_tuples [pivcol]) ;
391
392 Row_tuples [pivrow] = 0 ;
393 Col_tuples [pivcol] = 0 ;
394
395 DEBUG5 (("number of pivots prior to this one: "ID"\n", k)) ;
396 ASSERT (NON_PIVOTAL_ROW (pivrow)) ;
397 ASSERT (NON_PIVOTAL_COL (pivcol)) ;
398
399 /* save row and column inverse permutation */
400 k1 = ONES_COMPLEMENT (k) ;
401 Rperm [pivrow] = k1 ; /* aliased with Row_degree */
402 Cperm [pivcol] = k1 ; /* aliased with Col_degree */
403
404 ASSERT (!NON_PIVOTAL_ROW (pivrow)) ;
405 ASSERT (!NON_PIVOTAL_COL (pivcol)) ;
406
407 /* ---------------------------------------------------------------------- */
408 /* Keep track of the pivot order. This is the kth pivot row and column. */
409 /* ---------------------------------------------------------------------- */
410
411 /* keep track of pivot rows and columns in the LU, L, and U blocks */
412 ASSERT (fnpiv < MAXNB) ;
413 Work->Pivrow [fnpiv] = pivrow ;
414 Work->Pivcol [fnpiv] = pivcol ;
415
416 /* ====================================================================== */
417 /* === one step in the factorization is done ============================ */
418 /* ====================================================================== */
419
420 /* One more step is done, except for pending updates to the U and C blocks
421 * of this frontal matrix. Those are saved up, and applied by
422 * UMF_blas3_update when enough pivots have accumulated. Also, the
423 * LU factors for these pending pivots have not yet been stored. */
424
425 Work->fnpiv++ ;
426
427 #ifndef NDEBUG
428 DEBUG7 (("Current frontal matrix: (after pivcol scale)\n")) ;
429 UMF_dump_current_front (Numeric, Work, TRUE) ;
430 #endif
431
432 }
433