Mercurial > octave-nkf
comparison liboctave/UMFPACK/UMFPACK/Source/umf_scale_column.c @ 5164:57077d0ddc8e
[project @ 2005-02-25 19:55:24 by jwe]
author | jwe |
---|---|
date | Fri, 25 Feb 2005 19:55:28 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
5163:9f3299378193 | 5164:57077d0ddc8e |
---|---|
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 | |
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 | |
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 } |