Mercurial > octave-nkf

comparison liboctave/UMFPACK/UMFPACK/Source/umf_create_element.c @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[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_create_element =================================================== */
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 Factorization of a frontal matrix is complete. Create a new element for
13 later assembly into a subsequent frontal matrix. Returns TRUE if
14 successful, FALSE if out of memory.
15 */
16
17 #include "umf_internal.h"
18 #include "umf_mem_alloc_element.h"
19 #include "umf_mem_alloc_tail_block.h"
20 #include "umf_mem_free_tail_block.h"
21 #include "umf_get_memory.h"
22
23 /* ========================================================================== */
24 /* === copy_column ========================================================== */
25 /* ========================================================================== */
26
27 PRIVATE void copy_column (Int len, Entry *X, Entry *Y)
28 {
29 Int i ;
30 #pragma ivdep
31 for (i = 0 ; i < len ; i++)
32 {
33 Y [i] = X [i] ;
34 }
35 }
36
37 /* ========================================================================== */
38 /* === UMF_create_element =================================================== */
39 /* ========================================================================== */
40
41 GLOBAL Int UMF_create_element
42 (
43 NumericType *Numeric,
44 WorkType *Work,
45 SymbolicType *Symbolic
46 )
47 {
48 /* ---------------------------------------------------------------------- */
49 /* local variables */
50 /* ---------------------------------------------------------------------- */
51
52 Int j, col, row, *Fcols, *Frows, fnrows, fncols, *Cols, len, needunits, t1,
53 t2, size, e, i, *E, *Fcpos, *Frpos, *Rows, eloc, fnr_curr, f,
54 got_memory, *Row_tuples, *Row_degree, *Row_tlen, *Col_tuples, max_mark,
55 *Col_degree, *Col_tlen, nn, n_row, n_col, r2, c2, do_Fcpos ;
56 Entry *C, *Fcol ;
57 Element *ep ;
58 Unit *p, *Memory ;
59 Tuple *tp, *tp1, *tp2, tuple, *tpend ;
60 #ifndef NDEBUG
61 DEBUG2 (("FRONTAL WRAPUP\n")) ;
62 UMF_dump_current_front (Numeric, Work, TRUE) ;
63 #endif
64
65 /* ---------------------------------------------------------------------- */
66 /* get parameters */
67 /* ---------------------------------------------------------------------- */
68
69 ASSERT (Work->fnpiv == 0) ;
70 ASSERT (Work->fnzeros == 0) ;
71 Row_degree = Numeric->Rperm ;
72 Row_tuples = Numeric->Uip ;
73 Row_tlen = Numeric->Uilen ;
74 Col_degree = Numeric->Cperm ;
75 Col_tuples = Numeric->Lip ;
76 Col_tlen = Numeric->Lilen ;
77 n_row = Work->n_row ;
78 n_col = Work->n_col ;
79 nn = MAX (n_row, n_col) ;
80 Fcols = Work->Fcols ;
81 Frows = Work->Frows ;
82 Fcpos = Work->Fcpos ;
83 Frpos = Work->Frpos ;
84 Memory = Numeric->Memory ;
85 fncols = Work->fncols ;
86 fnrows = Work->fnrows ;
87
88 tp = (Tuple *) NULL ;
89 tp1 = (Tuple *) NULL ;
90 tp2 = (Tuple *) NULL ;
91
92 /* ---------------------------------------------------------------------- */
93 /* add the current frontal matrix to the degrees of each column */
94 /* ---------------------------------------------------------------------- */
95
96 if (!Symbolic->fixQ)
97 {
98 /* but only if the column ordering is not fixed */
99 #pragma ivdep
100 for (j = 0 ; j < fncols ; j++)
101 {
102 /* add the current frontal matrix to the degree */
103 ASSERT (Fcols [j] >= 0 && Fcols [j] < n_col) ;
104 Col_degree [Fcols [j]] += fnrows ;
105 }
106 }
107
108 /* ---------------------------------------------------------------------- */
109 /* add the current frontal matrix to the degrees of each row */
110 /* ---------------------------------------------------------------------- */
111
112 #pragma ivdep
113 for (i = 0 ; i < fnrows ; i++)
114 {
115 /* add the current frontal matrix to the degree */
116 ASSERT (Frows [i] >= 0 && Frows [i] < n_row) ;
117 Row_degree [Frows [i]] += fncols ;
118 }
119
120 /* ---------------------------------------------------------------------- */
121 /* Reset the external degree counters */
122 /* ---------------------------------------------------------------------- */
123
124 E = Work->E ;
125 max_mark = MAX_MARK (nn) ;
126
127 if (!Work->pivcol_in_front)
128 {
129 /* clear the external column degrees. no more Usons of current front */
130 Work->cdeg0 += (nn + 1) ;
131 if (Work->cdeg0 >= max_mark)
132 {
133 /* guard against integer overflow. This is very rare */
134 DEBUG1 (("Integer overflow, cdeg\n")) ;
135 Work->cdeg0 = 1 ;
136 #pragma ivdep
137 for (e = 1 ; e <= Work->nel ; e++)
138 {
139 if (E [e])
140 {
141 ep = (Element *) (Memory + E [e]) ;
142 ep->cdeg = 0 ;
143 }
144 }
145 }
146 }
147
148 if (!Work->pivrow_in_front)
149 {
150 /* clear the external row degrees. no more Lsons of current front */
151 Work->rdeg0 += (nn + 1) ;
152 if (Work->rdeg0 >= max_mark)
153 {
154 /* guard against integer overflow. This is very rare */
155 DEBUG1 (("Integer overflow, rdeg\n")) ;
156 Work->rdeg0 = 1 ;
157 #pragma ivdep
158 for (e = 1 ; e <= Work->nel ; e++)
159 {
160 if (E [e])
161 {
162 ep = (Element *) (Memory + E [e]) ;
163 ep->rdeg = 0 ;
164 }
165 }
166 }
167 }
168
169 /* ---------------------------------------------------------------------- */
170 /* clear row/col offsets */
171 /* ---------------------------------------------------------------------- */
172
173 if (!Work->pivrow_in_front)
174 {
175 #pragma ivdep
176 for (j = 0 ; j < fncols ; j++)
177 {
178 Fcpos [Fcols [j]] = EMPTY ;
179 }
180 }
181
182 if (!Work->pivcol_in_front)
183 {
184 #pragma ivdep
185 for (i = 0 ; i < fnrows ; i++)
186 {
187 Frpos [Frows [i]] = EMPTY ;
188 }
189 }
190
191 if (fncols <= 0 || fnrows <= 0)
192 {
193 /* no element to create */
194 DEBUG2 (("Element evaporation\n")) ;
195 Work->prior_element = EMPTY ;
196 return (TRUE) ;
197 }
198
199 /* ---------------------------------------------------------------------- */
200 /* create element for later assembly */
201 /* ---------------------------------------------------------------------- */
202
203 #ifndef NDEBUG
204 UMF_allocfail = FALSE ;
205 if (UMF_gprob > 0)
206 {
207 double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ;
208 DEBUG4 (("Check random %e %e\n", rrr, UMF_gprob)) ;
209 UMF_allocfail = rrr < UMF_gprob ;
210 if (UMF_allocfail) DEBUGm2 (("Random garbage collection (create)\n"));
211 }
212 #endif
213
214 needunits = 0 ;
215 got_memory = FALSE ;
216 eloc = UMF_mem_alloc_element (Numeric, fnrows, fncols, &Rows, &Cols, &C,
217 &needunits, &ep) ;
218
219 /* if UMF_get_memory needs to be called */
220 if (Work->do_grow)
221 {
222 /* full compaction of current frontal matrix, since UMF_grow_front will
223 * be called next anyway. */
224 r2 = fnrows ;
225 c2 = fncols ;
226 do_Fcpos = FALSE ;
227 }
228 else
229 {
230 /* partial compaction. */
231 r2 = MAX (fnrows, Work->fnrows_new + 1) ;
232 c2 = MAX (fncols, Work->fncols_new + 1) ;
233 /* recompute Fcpos if pivot row is in the front */
234 do_Fcpos = Work->pivrow_in_front ;
235 }
236
237 if (!eloc)
238 {
239 /* Do garbage collection, realloc, and try again. */
240 /* Compact the current front if it needs to grow anyway. */
241 /* Note that there are no pivot rows or columns in the current front */
242 DEBUGm3 (("get_memory from umf_create_element, 1\n")) ;
243 if (!UMF_get_memory (Numeric, Work, needunits, r2, c2, do_Fcpos))
244 {
245 /* :: out of memory in umf_create_element (1) :: */
246 DEBUGm4 (("out of memory: create element (1)\n")) ;
247 return (FALSE) ; /* out of memory */
248 }
249 got_memory = TRUE ;
250 Memory = Numeric->Memory ;
251 eloc = UMF_mem_alloc_element (Numeric, fnrows, fncols, &Rows, &Cols, &C,
252 &needunits, &ep) ;
253 ASSERT (eloc >= 0) ;
254 if (!eloc)
255 {
256 /* :: out of memory in umf_create_element (2) :: */
257 DEBUGm4 (("out of memory: create element (2)\n")) ;
258 return (FALSE) ; /* out of memory */
259 }
260 }
261
262 e = ++(Work->nel) ; /* get the name of this new frontal matrix */
263 Work->prior_element = e ;
264 DEBUG8 (("wrapup e "ID" nel "ID"\n", e, Work->nel)) ;
265
266 ASSERT (e > 0 && e < Work->elen) ;
267 ASSERT (E [e] == 0) ;
268 E [e] = eloc ;
269
270 if (Work->pivcol_in_front)
271 {
272 /* the new element is a Uson of the next frontal matrix */
273 ep->cdeg = Work->cdeg0 ;
274 }
275
276 if (Work->pivrow_in_front)
277 {
278 /* the new element is an Lson of the next frontal matrix */
279 ep->rdeg = Work->rdeg0 ;
280 }
281
282 /* ---------------------------------------------------------------------- */
283 /* copy frontal matrix into the new element */
284 /* ---------------------------------------------------------------------- */
285
286 #pragma ivdep
287 for (i = 0 ; i < fnrows ; i++)
288 {
289 Rows [i] = Frows [i] ;
290 }
291 #pragma ivdep
292 for (i = 0 ; i < fncols ; i++)
293 {
294 Cols [i] = Fcols [i] ;
295 }
296 Fcol = Work->Fcblock ;
297 DEBUG0 (("copy front "ID" by "ID"\n", fnrows, fncols)) ;
298 fnr_curr = Work->fnr_curr ;
299 ASSERT (fnr_curr >= 0 && fnr_curr % 2 == 1) ;
300 for (j = 0 ; j < fncols ; j++)
301 {
302 copy_column (fnrows, Fcol, C) ;
303 #if 0
304 #ifdef USE_NO_BLAS
305 copy_column (fnrows, Fcol, C) ;
306 #else
307 could also use BLAS-COPY (fnrows, Fcol, C) here, but it is typically
308 not as fast as the inlined copy_column subroutine, above.
309 #endif
310 for (i = 0 ; i < fnrows ; i++)
311 {
312 C [i] = Fcol [i] ;
313 }
314 #endif
315 Fcol += fnr_curr ;
316 C += fnrows ;
317 }
318
319 DEBUG8 (("element copied\n")) ;
320
321 /* ---------------------------------------------------------------------- */
322 /* add tuples for the new element */
323 /* ---------------------------------------------------------------------- */
324
325 tuple.e = e ;
326
327 if (got_memory)
328 {
329
330 /* ------------------------------------------------------------------ */
331 /* UMF_get_memory ensures enough space exists for each new tuple */
332 /* ------------------------------------------------------------------ */
333
334 /* place (e,f) in the element list of each column */
335 for (tuple.f = 0 ; tuple.f < fncols ; tuple.f++)
336 {
337 col = Fcols [tuple.f] ;
338 ASSERT (col >= 0 && col < n_col) ;
339 ASSERT (NON_PIVOTAL_COL (col)) ;
340 ASSERT (Col_tuples [col]) ;
341 tp = ((Tuple *) (Memory + Col_tuples [col])) + Col_tlen [col]++ ;
342 *tp = tuple ;
343 }
344
345 /* ------------------------------------------------------------------ */
346
347 /* place (e,f) in the element list of each row */
348 for (tuple.f = 0 ; tuple.f < fnrows ; tuple.f++)
349 {
350 row = Frows [tuple.f] ;
351 ASSERT (row >= 0 && row < n_row) ;
352 ASSERT (NON_PIVOTAL_ROW (row)) ;
353 ASSERT (Row_tuples [row]) ;
354 tp = ((Tuple *) (Memory + Row_tuples [row])) + Row_tlen [row]++ ;
355 *tp = tuple ;
356 }
357
358 }
359 else
360 {
361
362 /* ------------------------------------------------------------------ */
363 /* place (e,f) in the element list of each column */
364 /* ------------------------------------------------------------------ */
365
366 /* might not have enough space for each tuple */
367
368 for (tuple.f = 0 ; tuple.f < fncols ; tuple.f++)
369 {
370 col = Fcols [tuple.f] ;
371 ASSERT (col >= 0 && col < n_col) ;
372 ASSERT (NON_PIVOTAL_COL (col)) ;
373 t1 = Col_tuples [col] ;
374 DEBUG1 (("Placing on col:"ID" , tuples at "ID"\n",
375 col, Col_tuples [col])) ;
376
377 size = 0 ;
378 len = 0 ;
379
380 if (t1)
381 {
382 p = Memory + t1 ;
383 tp = (Tuple *) p ;
384 size = GET_BLOCK_SIZE (p) ;
385 len = Col_tlen [col] ;
386 tp2 = tp + len ;
387 }
388
389 needunits = UNITS (Tuple, len + 1) ;
390 DEBUG1 (("len: "ID" size: "ID" needunits: "ID"\n",
391 len, size, needunits));
392
393 if (needunits > size && t1)
394 {
395 /* prune the tuples */
396 tp1 = tp ;
397 tp2 = tp ;
398 tpend = tp + len ;
399 for ( ; tp < tpend ; tp++)
400 {
401 e = tp->e ;
402 ASSERT (e > 0 && e <= Work->nel) ;
403 if (!E [e]) continue ; /* element already deallocated */
404 f = tp->f ;
405 p = Memory + E [e] ;
406 ep = (Element *) p ;
407 p += UNITS (Element, 1) ;
408 Cols = (Int *) p ;
409 ;
410 if (Cols [f] == EMPTY) continue ; /* already assembled */
411 ASSERT (col == Cols [f]) ;
412 *tp2++ = *tp ; /* leave the tuple in the list */
413 }
414 len = tp2 - tp1 ;
415 Col_tlen [col] = len ;
416 needunits = UNITS (Tuple, len + 1) ;
417 }
418
419 if (needunits > size)
420 {
421 /* no room exists - reallocate elsewhere */
422 DEBUG1 (("REALLOCATE Col: "ID", size "ID" to "ID"\n",
423 col, size, 2*needunits)) ;
424
425 #ifndef NDEBUG
426 UMF_allocfail = FALSE ;
427 if (UMF_gprob > 0) /* a double relop, but ignore NaN case */
428 {
429 double rrr = ((double) (rand ( ))) /
430 (((double) RAND_MAX) + 1) ;
431 DEBUG1 (("Check random %e %e\n", rrr, UMF_gprob)) ;
432 UMF_allocfail = rrr < UMF_gprob ;
433 if (UMF_allocfail) DEBUGm2 (("Random gar. (col tuple)\n")) ;
434 }
435 #endif
436
437 needunits = MIN (2*needunits, (Int) UNITS (Tuple, nn)) ;
438 t2 = UMF_mem_alloc_tail_block (Numeric, needunits) ;
439 if (!t2)
440 {
441 /* :: get memory in umf_create_element (1) :: */
442 /* get memory, reconstruct all tuple lists, and return */
443 /* Compact the current front if it needs to grow anyway. */
444 /* Note: no pivot rows or columns in the current front */
445 DEBUGm4 (("get_memory from umf_create_element, 1\n")) ;
446 return (UMF_get_memory (Numeric, Work, 0, r2, c2,do_Fcpos));
447 }
448 Col_tuples [col] = t2 ;
449 tp2 = (Tuple *) (Memory + t2) ;
450 if (t1)
451 {
452 for (i = 0 ; i < len ; i++)
453 {
454 *tp2++ = *tp1++ ;
455 }
456 UMF_mem_free_tail_block (Numeric, t1) ;
457 }
458 }
459
460 /* place the new (e,f) tuple in the element list of the column */
461 Col_tlen [col]++ ;
462 *tp2 = tuple ;
463 }
464
465 /* ------------------------------------------------------------------ */
466 /* place (e,f) in the element list of each row */
467 /* ------------------------------------------------------------------ */
468
469 for (tuple.f = 0 ; tuple.f < fnrows ; tuple.f++)
470 {
471 row = Frows [tuple.f] ;
472 ASSERT (row >= 0 && row < n_row) ;
473 ASSERT (NON_PIVOTAL_ROW (row)) ;
474 t1 = Row_tuples [row] ;
475 DEBUG1 (("Placing on row:"ID" , tuples at "ID"\n",
476 row, Row_tuples [row])) ;
477
478 size = 0 ;
479 len = 0 ;
480 if (t1)
481 {
482 p = Memory + t1 ;
483 tp = (Tuple *) p ;
484 size = GET_BLOCK_SIZE (p) ;
485 len = Row_tlen [row] ;
486 tp2 = tp + len ;
487 }
488
489 needunits = UNITS (Tuple, len + 1) ;
490 DEBUG1 (("len: "ID" size: "ID" needunits: "ID"\n",
491 len, size, needunits)) ;
492
493 if (needunits > size && t1)
494 {
495 /* prune the tuples */
496 tp1 = tp ;
497 tp2 = tp ;
498 tpend = tp + len ;
499 for ( ; tp < tpend ; tp++)
500 {
501 e = tp->e ;
502 ASSERT (e > 0 && e <= Work->nel) ;
503 if (!E [e])
504 {
505 continue ; /* element already deallocated */
506 }
507 f = tp->f ;
508 p = Memory + E [e] ;
509 ep = (Element *) p ;
510 p += UNITS (Element, 1) ;
511 Cols = (Int *) p ;
512 Rows = Cols + (ep->ncols) ;
513 if (Rows [f] == EMPTY) continue ; /* already assembled */
514 ASSERT (row == Rows [f]) ;
515 *tp2++ = *tp ; /* leave the tuple in the list */
516 }
517 len = tp2 - tp1 ;
518 Row_tlen [row] = len ;
519 needunits = UNITS (Tuple, len + 1) ;
520 }
521
522 if (needunits > size)
523 {
524 /* no room exists - reallocate elsewhere */
525 DEBUG1 (("REALLOCATE Row: "ID", size "ID" to "ID"\n",
526 row, size, 2*needunits)) ;
527
528 #ifndef NDEBUG
529 UMF_allocfail = FALSE ;
530 if (UMF_gprob > 0) /* a double relop, but ignore NaN case */
531 {
532 double rrr = ((double) (rand ( ))) /
533 (((double) RAND_MAX) + 1) ;
534 DEBUG1 (("Check random %e %e\n", rrr, UMF_gprob)) ;
535 UMF_allocfail = rrr < UMF_gprob ;
536 if (UMF_allocfail) DEBUGm2 (("Random gar. (row tuple)\n")) ;
537 }
538 #endif
539
540 needunits = MIN (2*needunits, (Int) UNITS (Tuple, nn)) ;
541 t2 = UMF_mem_alloc_tail_block (Numeric, needunits) ;
542 if (!t2)
543 {
544 /* :: get memory in umf_create_element (2) :: */
545 /* get memory, reconstruct all tuple lists, and return */
546 /* Compact the current front if it needs to grow anyway. */
547 /* Note: no pivot rows or columns in the current front */
548 DEBUGm4 (("get_memory from umf_create_element, 2\n")) ;
549 return (UMF_get_memory (Numeric, Work, 0, r2, c2,do_Fcpos));
550 }
551 Row_tuples [row] = t2 ;
552 tp2 = (Tuple *) (Memory + t2) ;
553 if (t1)
554 {
555 for (i = 0 ; i < len ; i++)
556 {
557 *tp2++ = *tp1++ ;
558 }
559 UMF_mem_free_tail_block (Numeric, t1) ;
560 }
561 }
562
563 /* place the new (e,f) tuple in the element list of the row */
564 Row_tlen [row]++ ;
565 *tp2 = tuple ;
566 }
567
568 }
569
570 /* ---------------------------------------------------------------------- */
571
572 #ifndef NDEBUG
573 DEBUG1 (("Done extending\nFINAL: element row pattern: len="ID"\n", fncols));
574 for (j = 0 ; j < fncols ; j++) DEBUG1 ((""ID"\n", Fcols [j])) ;
575 DEBUG1 (("FINAL: element col pattern: len="ID"\n", fnrows)) ;
576 for (j = 0 ; j < fnrows ; j++) DEBUG1 ((""ID"\n", Frows [j])) ;
577 for (j = 0 ; j < fncols ; j++)
578 {
579 col = Fcols [j] ;
580 ASSERT (col >= 0 && col < n_col) ;
581 UMF_dump_rowcol (1, Numeric, Work, col, !Symbolic->fixQ) ;
582 }
583 for (j = 0 ; j < fnrows ; j++)
584 {
585 row = Frows [j] ;
586 ASSERT (row >= 0 && row < n_row) ;
587 UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
588 }
589 if (n_row < 1000 && n_col < 1000)
590 {
591 UMF_dump_memory (Numeric) ;
592 }
593 DEBUG1 (("New element, after filling with stuff: "ID"\n", e)) ;
594 UMF_dump_element (Numeric, Work, e, TRUE) ;
595 if (nn < 1000)
596 {
597 DEBUG4 (("Matrix dump, after New element: "ID"\n", e)) ;
598 UMF_dump_matrix (Numeric, Work, TRUE) ;
599 }
600 DEBUG3 (("FRONTAL WRAPUP DONE\n")) ;
601 #endif
602
603 return (TRUE) ;
604 }