1 /* ========================================================================== */
2 /* === UMF_kernel =========================================================== */
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 Primary factorization routine. Called by UMFPACK_numeric.
12 Returns:
13 UMFPACK_OK if successful,
14 UMFPACK_ERROR_out_of_memory if out of memory, or
15 UMFPACK_ERROR_different_pattern if pattern of matrix (Ap and/or Ai)
16 has changed since the call to UMFPACK_*symbolic.
17 */
18
19 #include "umf_internal.h"
20 #include "umf_kernel.h"
21 #include "umf_kernel_init.h"
22 #include "umf_init_front.h"
23 #include "umf_start_front.h"
24 #include "umf_assemble.h"
25 #include "umf_scale_column.h"
26 #include "umf_local_search.h"
27 #include "umf_create_element.h"
28 #include "umf_extend_front.h"
29 #include "umf_blas3_update.h"
30 #include "umf_store_lu.h"
31 #include "umf_kernel_wrapup.h"
32
33 /* perform an action, and return if out of memory */
34 #define DO(action) { if (! (action)) { return (UMFPACK_ERROR_out_of_memory) ; }}
35
UMF_kernel(const Int Ap[],const Int Ai[],const double Ax[],const double Az[],NumericType * Numeric,WorkType * Work,SymbolicType * Symbolic)36 GLOBAL Int UMF_kernel
37 (
38 const Int Ap [ ],
39 const Int Ai [ ],
40 const double Ax [ ],
41 #ifdef COMPLEX
42 const double Az [ ],
43 #endif
44 NumericType *Numeric,
45 WorkType *Work,
46 SymbolicType *Symbolic
47 )
48 {
49
50 /* ---------------------------------------------------------------------- */
51 /* local variables */
52 /* ---------------------------------------------------------------------- */
53
54 Int j, f1, f2, chain, nchains, *Chain_start, status, fixQ, evaporate,
55 *Front_npivcol, jmax, nb, drop ;
56
57 /* ---------------------------------------------------------------------- */
58 /* initialize memory space and load the matrix. Optionally scale. */
59 /* ---------------------------------------------------------------------- */
60
61 if (!UMF_kernel_init (Ap, Ai, Ax,
62 #ifdef COMPLEX
63 Az,
64 #endif
65 Numeric, Work, Symbolic))
66 {
67 /* UMF_kernel_init is guaranteed to succeed, since UMFPACK_numeric */
68 /* either allocates enough space or if not, UMF_kernel does not get */
69 /* called. So running out of memory here is a fatal error, and means */
70 /* that the user changed Ap and/or Ai since the call to */
71 /* UMFPACK_*symbolic. */
72 DEBUGm4 (("kernel init failed\n")) ;
73 return (UMFPACK_ERROR_different_pattern) ;
74 }
75
76 /* ---------------------------------------------------------------------- */
77 /* get the symbolic factorization */
78 /* ---------------------------------------------------------------------- */
79
80 nchains = Symbolic->nchains ;
81 Chain_start = Symbolic->Chain_start ;
82 Front_npivcol = Symbolic->Front_npivcol ;
83 nb = Symbolic->nb ;
84 fixQ = Symbolic->fixQ ;
85 drop = Numeric->droptol > 0.0 ;
86
87 #ifndef NDEBUG
88 for (chain = 0 ; chain < nchains ; chain++)
89 {
90 Int i ;
91 f1 = Chain_start [chain] ;
92 f2 = Chain_start [chain+1] - 1 ;
93 DEBUG1 (("\nCHain: "ID" start "ID" end "ID"\n", chain, f1, f2)) ;
94 for (i = f1 ; i <= f2 ; i++)
95 {
96 DEBUG1 (("Front "ID", npivcol "ID"\n", i, Front_npivcol [i])) ;
97 }
98 }
99 #endif
100
101 /* ---------------------------------------------------------------------- */
102 /* factorize each chain of frontal matrices */
103 /* ---------------------------------------------------------------------- */
104
105 for (chain = 0 ; chain < nchains ; chain++)
106 {
107 f1 = Chain_start [chain] ;
108 f2 = Chain_start [chain+1] - 1 ;
109
110 /* ------------------------------------------------------------------ */
111 /* get the initial frontal matrix size for this chain */
112 /* ------------------------------------------------------------------ */
113
114 DO (UMF_start_front (chain, Numeric, Work, Symbolic)) ;
115
116 /* ------------------------------------------------------------------ */
117 /* factorize each front in the chain */
118 /* ------------------------------------------------------------------ */
119
120 for (Work->frontid = f1 ; Work->frontid <= f2 ; Work->frontid++)
121 {
122
123 /* -------------------------------------------------------------- */
124 /* Initialize the pivot column candidate set */
125 /* -------------------------------------------------------------- */
126
127 Work->ncand = Front_npivcol [Work->frontid] ;
128 Work->lo = Work->nextcand ;
129 Work->hi = Work->nextcand + Work->ncand - 1 ;
130 jmax = MIN (MAX_CANDIDATES, Work->ncand) ;
131 DEBUGm1 ((">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Starting front "
132 ID", npivcol: "ID"\n", Work->frontid, Work->ncand)) ;
133 if (fixQ)
134 {
135 /* do not modify the column order */
136 jmax = 1 ;
137 }
138 DEBUGm1 (("Initial candidates: ")) ;
139 for (j = 0 ; j < jmax ; j++)
140 {
141 DEBUGm1 ((" "ID, Work->nextcand)) ;
142 ASSERT (Work->nextcand <= Work->hi) ;
143 Work->Candidates [j] = Work->nextcand++ ;
144 }
145 Work->nCandidates = jmax ;
146 DEBUGm1 (("\n")) ;
147
148 /* -------------------------------------------------------------- */
149 /* Assemble and factorize the current frontal matrix */
150 /* -------------------------------------------------------------- */
151
152 while (Work->ncand > 0)
153 {
154
155 /* ---------------------------------------------------------- */
156 /* get the pivot row and column */
157 /* ---------------------------------------------------------- */
158
159 status = UMF_local_search (Numeric, Work, Symbolic) ;
160 if (status == UMFPACK_ERROR_different_pattern)
161 {
162 /* :: pattern change detected in umf_local_search :: */
163 /* input matrix has changed since umfpack_*symbolic */
164 DEBUGm4 (("local search failed\n")) ;
165 return (UMFPACK_ERROR_different_pattern) ;
166 }
167 if (status == UMFPACK_WARNING_singular_matrix)
168 {
169 /* no pivot found, discard and try again */
170 continue ;
171 }
172
173 /* ---------------------------------------------------------- */
174 /* update if front not extended or too many zeros in L,U */
175 /* ---------------------------------------------------------- */
176
177 if (Work->do_update)
178 {
179 UMF_blas3_update (Work) ;
180 if (drop)
181 {
182 DO (UMF_store_lu_drop (Numeric, Work)) ;
183 }
184 else
185 {
186 DO (UMF_store_lu (Numeric, Work)) ;
187 }
188 }
189
190 /* ---------------------------------------------------------- */
191 /* extend the frontal matrix, or start a new one */
192 /* ---------------------------------------------------------- */
193
194 if (Work->do_extend)
195 {
196 /* extend the current front */
197 DO (UMF_extend_front (Numeric, Work)) ;
198 }
199 else
200 {
201 /* finish the current front (if any) and start a new one */
202 DO (UMF_create_element (Numeric, Work, Symbolic)) ;
203 DO (UMF_init_front (Numeric, Work)) ;
204 }
205
206 /* ---------------------------------------------------------- */
207 /* Numerical & symbolic assembly into current frontal matrix */
208 /* ---------------------------------------------------------- */
209
210 if (fixQ)
211 {
212 UMF_assemble_fixq (Numeric, Work) ;
213 }
214 else
215 {
216 UMF_assemble (Numeric, Work) ;
217 }
218
219 /* ---------------------------------------------------------- */
220 /* scale the pivot column */
221 /* ---------------------------------------------------------- */
222
223 UMF_scale_column (Numeric, Work) ;
224
225 /* ---------------------------------------------------------- */
226 /* Numerical update if enough pivots accumulated */
227 /* ---------------------------------------------------------- */
228
229 evaporate = Work->fnrows == 0 || Work->fncols == 0 ;
230 if (Work->fnpiv >= nb || evaporate)
231 {
232 UMF_blas3_update (Work) ;
233 if (drop)
234 {
235 DO (UMF_store_lu_drop (Numeric, Work)) ;
236 }
237 else
238 {
239 DO (UMF_store_lu (Numeric, Work)) ;
240 }
241
242 }
243
244 Work->pivrow_in_front = FALSE ;
245 Work->pivcol_in_front = FALSE ;
246
247 /* ---------------------------------------------------------- */
248 /* If front is empty, evaporate it */
249 /* ---------------------------------------------------------- */
250
251 if (evaporate)
252 {
253 /* This does not create an element, just evaporates it.
254 * It ensures that a front is not 0-by-c or r-by-0. No
255 * memory is allocated, so it is guaranteed to succeed. */
256 (void) UMF_create_element (Numeric, Work, Symbolic) ;
257 Work->fnrows = 0 ;
258 Work->fncols = 0 ;
259 }
260 }
261 }
262
263 /* ------------------------------------------------------------------
264 * Wrapup the current frontal matrix. This is the last in a chain
265 * in the column elimination tree. The next frontal matrix
266 * cannot overlap with the current one, which will be its sibling
267 * in the column etree.
268 * ------------------------------------------------------------------ */
269
270 UMF_blas3_update (Work) ;
271 if (drop)
272 {
273 DO (UMF_store_lu_drop (Numeric, Work)) ;
274 }
275 else
276 {
277 DO (UMF_store_lu (Numeric, Work)) ;
278 }
279 Work->fnrows_new = Work->fnrows ;
280 Work->fncols_new = Work->fncols ;
281 DO (UMF_create_element (Numeric, Work, Symbolic)) ;
282
283 /* ------------------------------------------------------------------ */
284 /* current front is now empty */
285 /* ------------------------------------------------------------------ */
286
287 Work->fnrows = 0 ;
288 Work->fncols = 0 ;
289 }
290
291 /* ---------------------------------------------------------------------- */
292 /* end the last Lchain and Uchain and finalize the LU factors */
293 /* ---------------------------------------------------------------------- */
294
295 UMF_kernel_wrapup (Numeric, Symbolic, Work) ;
296
297 /* note that the matrix may be singular (this is OK) */
298 return (UMFPACK_OK) ;
299 }
300