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