1 /* ========================================================================== */
2 /* === UMF_start_front ====================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Copyright (c) Timothy A. Davis, CISE, */
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 /* Allocate the initial frontal matrix working array for a single chain. The
12 * front does not have to be big enough, since if it's too small it will get
13 * reallocated. The size computed here is just an estimate. */
14
15 #include "umf_internal.h"
16 #include "umf_start_front.h"
17 #include "umf_grow_front.h"
18
UMF_start_front(Int chain,NumericType * Numeric,WorkType * Work,SymbolicType * Symbolic)19 GLOBAL Int UMF_start_front /* returns TRUE if successful, FALSE otherwise */
20 (
21 Int chain,
22 NumericType *Numeric,
23 WorkType *Work,
24 SymbolicType *Symbolic
25 )
26 {
27 double maxbytes ;
28 Int fnrows_max, fncols_max, fnr2, fnc2, fsize, fcurr_size, maxfrsize,
29 overflow, nb, f, cdeg ;
30
31 nb = Symbolic->nb ;
32 fnrows_max = Symbolic->Chain_maxrows [chain] ;
33 fncols_max = Symbolic->Chain_maxcols [chain] ;
34
35 DEBUGm2 (("Start Front for chain "ID". fnrows_max "ID" fncols_max "ID"\n",
36 chain, fnrows_max, fncols_max)) ;
37
38 Work->fnrows_max = fnrows_max ;
39 Work->fncols_max = fncols_max ;
40 Work->any_skip = FALSE ;
41
42 maxbytes = sizeof (Entry) *
43 (double) (fnrows_max + nb) * (double) (fncols_max + nb) ;
44 fcurr_size = Work->fcurr_size ;
45
46 if (Symbolic->prefer_diagonal)
47 {
48 /* Get a rough upper bound on the degree of the first pivot column in
49 * this front. Note that Col_degree is not maintained if diagonal
50 * pivoting is preferred. For most matrices, the first pivot column
51 * of the first frontal matrix of a new chain has only one tuple in
52 * it anyway, so this bound is exact in that case. */
53 Int col, tpi, e, *E, *Col_tuples, *Col_tlen, *Cols ;
54 Tuple *tp, *tpend ;
55 Unit *Memory, *p ;
56 Element *ep ;
57 E = Work->E ;
58 Memory = Numeric->Memory ;
59 Col_tuples = Numeric->Lip ;
60 Col_tlen = Numeric->Lilen ;
61 col = Work->nextcand ;
62 tpi = Col_tuples [col] ;
63 tp = (Tuple *) Memory + tpi ;
64 tpend = tp + Col_tlen [col] ;
65 cdeg = 0 ;
66 DEBUGm3 (("\n=============== start front: col "ID" tlen "ID"\n",
67 col, Col_tlen [col])) ;
68 for ( ; tp < tpend ; tp++)
69 {
70 DEBUG1 (("Tuple ("ID","ID")\n", tp->e, tp->f)) ;
71 e = tp->e ;
72 if (!E [e]) continue ;
73 f = tp->f ;
74 p = Memory + E [e] ;
75 ep = (Element *) p ;
76 p += UNITS (Element, 1) ;
77 Cols = (Int *) p ;
78 if (Cols [f] == EMPTY) continue ;
79 DEBUG1 ((" nrowsleft "ID"\n", ep->nrowsleft)) ;
80 cdeg += ep->nrowsleft ;
81 }
82 #ifndef NDEBUG
83 DEBUGm3 (("start front cdeg: "ID" col "ID"\n", cdeg, col)) ;
84 UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
85 #endif
86
87 /* cdeg is now the rough upper bound on the degree of the next pivot
88 * column. */
89
90 /* If AMD was called, we know the maximum number of nonzeros in any
91 * column of L. Use this as an upper bound for cdeg, but add 2 to
92 * account for a small amount of off-diagonal pivoting. */
93 if (Symbolic->amd_dmax > 0)
94 {
95 cdeg = MIN (cdeg, Symbolic->amd_dmax) ;
96 }
97
98 /* Increase it to account for larger columns later on.
99 * Also ensure that it's larger than zero. */
100 cdeg += 2 ;
101
102 /* cdeg cannot be larger than fnrows_max */
103 cdeg = MIN (cdeg, fnrows_max) ;
104
105 }
106 else
107 {
108 /* don't do the above cdeg computation */
109 cdeg = 0 ;
110 }
111
112 DEBUGm2 (("fnrows max "ID" fncols_max "ID"\n", fnrows_max, fncols_max)) ;
113
114 /* the current frontal matrix is empty */
115 ASSERT (Work->fnrows == 0 && Work->fncols == 0 && Work->fnpiv == 0) ;
116
117 /* maximum row dimension is always odd, to avoid bad cache effects */
118 ASSERT (fnrows_max >= 0) ;
119 ASSERT (fnrows_max % 2 == 1) ;
120
121 /* ----------------------------------------------------------------------
122 * allocate working array for current frontal matrix:
123 * minimum size: 1-by-1
124 * maximum size: fnrows_max-by-fncols_max
125 * desired size:
126 *
127 * if Numeric->front_alloc_init >= 0:
128 *
129 * for unsymmetric matrices:
130 * Numeric->front_alloc_init * (fnrows_max-by-fncols_max)
131 *
132 * for symmetric matrices (diagonal pivoting preference, actually):
133 * Numeric->front_alloc_init * (fnrows_max-by-fncols_max), or
134 * cdeg*cdeg, whichever is smaller.
135 *
136 * if Numeric->front_alloc_init < 0:
137 * allocate a front of size -Numeric->front_alloc_init.
138 *
139 * Allocate the whole thing if it's small (less than 2*nb^2). Make sure the
140 * leading dimension of the frontal matrix is odd.
141 *
142 * Also allocate the nb-by-nb LU block, the dr-by-nb L block, and the
143 * nb-by-dc U block.
144 * ---------------------------------------------------------------------- */
145
146 /* get the maximum front size, avoiding integer overflow */
147 overflow = INT_OVERFLOW (maxbytes) ;
148 if (overflow)
149 {
150 /* :: int overflow, max front size :: */
151 maxfrsize = Int_MAX / sizeof (Entry) ;
152 }
153 else
154 {
155 maxfrsize = (fnrows_max + nb) * (fncols_max + nb) ;
156 }
157 ASSERT (!INT_OVERFLOW ((double) maxfrsize * sizeof (Entry))) ;
158
159 if (Numeric->front_alloc_init < 0)
160 {
161 /* allocate a front of -Numeric->front_alloc_init entries */
162 fsize = -Numeric->front_alloc_init ;
163 fsize = MAX (1, fsize) ;
164 }
165 else
166 {
167 if (INT_OVERFLOW (Numeric->front_alloc_init * maxbytes))
168 {
169 /* :: int overflow, requested front size :: */
170 fsize = Int_MAX / sizeof (Entry) ;
171 }
172 else
173 {
174 fsize = Numeric->front_alloc_init * maxfrsize ;
175 }
176
177 if (cdeg > 0)
178 {
179 /* diagonal pivoting is in use. cdeg was computed above */
180 Int fsize2 ;
181
182 /* add the L and U blocks */
183 cdeg += nb ;
184
185 if (INT_OVERFLOW (((double) cdeg * (double) cdeg) * sizeof (Entry)))
186 {
187 /* :: int overflow, symmetric front size :: */
188 fsize2 = Int_MAX / sizeof (Entry) ;
189 }
190 else
191 {
192 fsize2 = MAX (cdeg * cdeg, fcurr_size) ;
193 }
194 fsize = MIN (fsize, fsize2) ;
195 }
196 }
197
198 fsize = MAX (fsize, 2*nb*nb) ;
199
200 /* fsize and maxfrsize are now safe from integer overflow. They both
201 * include the size of the pivot blocks. */
202 ASSERT (!INT_OVERFLOW ((double) fsize * sizeof (Entry))) ;
203
204 Work->fnrows_new = 0 ;
205 Work->fncols_new = 0 ;
206
207 /* desired size is fnr2-by-fnc2 (includes L and U blocks): */
208 DEBUGm2 ((" fsize "ID" fcurr_size "ID"\n", fsize, fcurr_size)) ;
209 DEBUGm2 ((" maxfrsize "ID" fnr_curr "ID" fnc_curr "ID"\n", maxfrsize,
210 Work->fnr_curr, Work->fnc_curr)) ;
211
212 if (fsize >= maxfrsize && !overflow)
213 {
214 /* max working array is small, allocate all of it */
215 fnr2 = fnrows_max + nb ;
216 fnc2 = fncols_max + nb ;
217 fsize = maxfrsize ;
218 DEBUGm1 ((" sufficient for ("ID"+"ID")-by-("ID"+"ID")\n",
219 fnrows_max, nb, fncols_max, nb)) ;
220 }
221 else
222 {
223 /* allocate a smaller working array */
224 if (fnrows_max <= fncols_max)
225 {
226 fnr2 = (Int) sqrt ((double) fsize) ;
227 /* make sure fnr2 is odd */
228 fnr2 = MAX (fnr2, 1) ;
229 if (fnr2 % 2 == 0) fnr2++ ;
230 fnr2 = MIN (fnr2, fnrows_max + nb) ;
231 fnc2 = fsize / fnr2 ;
232 }
233 else
234 {
235 fnc2 = (Int) sqrt ((double) fsize) ;
236 fnc2 = MIN (fnc2, fncols_max + nb) ;
237 fnr2 = fsize / fnc2 ;
238 /* make sure fnr2 is odd */
239 fnr2 = MAX (fnr2, 1) ;
240 if (fnr2 % 2 == 0)
241 {
242 fnr2++ ;
243 fnc2 = fsize / fnr2 ;
244 }
245 }
246 DEBUGm1 ((" smaller "ID"-by-"ID"\n", fnr2, fnc2)) ;
247 }
248 fnr2 = MIN (fnr2, fnrows_max + nb) ;
249 fnc2 = MIN (fnc2, fncols_max + nb) ;
250 ASSERT (fnr2 % 2 == 1) ;
251 ASSERT (fnr2 * fnc2 <= fsize) ;
252
253 fnr2 -= nb ;
254 fnc2 -= nb ;
255 ASSERT (fnr2 >= 0) ;
256 ASSERT (fnc2 >= 0) ;
257
258 if (fsize > fcurr_size)
259 {
260 DEBUGm1 ((" Grow front \n")) ;
261 Work->do_grow = TRUE ;
262 if (!UMF_grow_front (Numeric, fnr2, fnc2, Work, -1))
263 {
264 /* since the minimum front size is 1-by-1, it would be nearly
265 * impossible to run out of memory here. */
266 DEBUGm4 (("out of memory: start front\n")) ;
267 return (FALSE) ;
268 }
269 }
270 else
271 {
272 /* use the existing front */
273 DEBUGm1 ((" existing front ok\n")) ;
274 Work->fnr_curr = fnr2 ;
275 Work->fnc_curr = fnc2 ;
276 Work->Flblock = Work->Flublock + nb * nb ;
277 Work->Fublock = Work->Flblock + nb * fnr2 ;
278 Work->Fcblock = Work->Fublock + nb * fnc2 ;
279 }
280 ASSERT (Work->Flblock == Work->Flublock + Work->nb*Work->nb) ;
281 ASSERT (Work->Fublock == Work->Flblock + Work->fnr_curr*Work->nb) ;
282 ASSERT (Work->Fcblock == Work->Fublock + Work->nb*Work->fnc_curr) ;
283 return (TRUE) ;
284 }
285