1 /* ========================================================================== */
2 /* === Cholesky/t_cholmod_rowfac ============================================ */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Template routine for cholmod_rowfac. Supports any numeric xtype
10 * (real, complex, or zomplex).
11 *
12 * workspace: Iwork (n), Flag (n), Xwork (n if real, 2*n if complex)
13 */
14
15 #include "cholmod_template.h"
16
17 #ifdef MASK
TEMPLATE(cholmod_rowfac_mask)18 static int TEMPLATE (cholmod_rowfac_mask)
19 #else
20 static int TEMPLATE (cholmod_rowfac)
21 #endif
22 (
23 /* ---- input ---- */
24 cholmod_sparse *A, /* matrix to factorize */
25 cholmod_sparse *F, /* used for A*A' case only. F=A' or A(:,f)' */
26 double beta [2], /* factorize beta*I+A or beta*I+AA' (beta [0] only) */
27 size_t kstart, /* first row to factorize */
28 size_t kend, /* last row to factorize is kend-1 */
29 #ifdef MASK
30 /* These inputs are used for cholmod_rowfac_mask only */
31 Int *mask, /* size A->nrow. if mask[i] >= maskmark
32 then W(i) is set to zero */
33 Int maskmark,
34 Int *RLinkUp, /* size A->nrow. link list of rows to compute */
35 #endif
36 /* ---- in/out --- */
37 cholmod_factor *L,
38 /* --------------- */
39 cholmod_common *Common
40 )
41 {
42 double yx [2], lx [2], fx [2], dk [1], di [1], fl = 0 ;
43 #ifdef ZOMPLEX
44 double yz [1], lz [1], fz [1] ;
45 #endif
46 double *Ax, *Az, *Lx, *Lz, *Wx, *Wz, *Fx, *Fz ;
47 Int *Ap, *Anz, *Ai, *Lp, *Lnz, *Li, *Lnext, *Flag, *Stack, *Fp, *Fi, *Fnz,
48 *Iwork ;
49 Int i, p, k, t, pf, pfend, top, s, mark, pend, n, lnz, is_ll, multadds,
50 use_dbound, packed, stype, Fpacked, sorted, nzmax, len, parent ;
51 #ifndef REAL
52 Int dk_imaginary ;
53 #endif
54
55 /* ---------------------------------------------------------------------- */
56 /* get inputs */
57 /* ---------------------------------------------------------------------- */
58
59 PRINT1 (("\nin cholmod_rowfac, kstart %d kend %d stype %d\n",
60 kstart, kend, A->stype)) ;
61 DEBUG (CHOLMOD(dump_factor) (L, "Initial L", Common)) ;
62
63 n = A->nrow ;
64 stype = A->stype ;
65
66 if (stype > 0)
67 {
68 /* symmetric upper case: F is not needed. It may be NULL */
69 Fp = NULL ;
70 Fi = NULL ;
71 Fx = NULL ;
72 Fz = NULL ;
73 Fnz = NULL ;
74 Fpacked = TRUE ;
75 }
76 else
77 {
78 /* unsymmetric case: F is required. */
79 Fp = F->p ;
80 Fi = F->i ;
81 Fx = F->x ;
82 Fz = F->z ;
83 Fnz = F->nz ;
84 Fpacked = F->packed ;
85 }
86
87 Ap = A->p ; /* size A->ncol+1, column pointers of A */
88 Ai = A->i ; /* size nz = Ap [A->ncol], row indices of A */
89 Ax = A->x ; /* size nz, numeric values of A */
90 Az = A->z ;
91 Anz = A->nz ;
92 packed = A->packed ;
93 sorted = A->sorted ;
94
95 use_dbound = IS_GT_ZERO (Common->dbound) ;
96
97 /* get the current factors L (and D for LDL'); allocate space if needed */
98 is_ll = L->is_ll ;
99 if (L->xtype == CHOLMOD_PATTERN)
100 {
101 /* ------------------------------------------------------------------ */
102 /* L is symbolic only; allocate and initialize L (and D for LDL') */
103 /* ------------------------------------------------------------------ */
104
105 /* workspace: none */
106 CHOLMOD(change_factor) (A->xtype, is_ll, FALSE, FALSE, TRUE, L, Common);
107 if (Common->status < CHOLMOD_OK)
108 {
109 /* out of memory */
110 return (FALSE) ;
111 }
112 ASSERT (L->minor == (size_t) n) ;
113 }
114 else if (kstart == 0 && kend == (size_t) n)
115 {
116 /* ------------------------------------------------------------------ */
117 /* refactorization; reset L->nz and L->minor to restart factorization */
118 /* ------------------------------------------------------------------ */
119
120 L->minor = n ;
121 Lnz = L->nz ;
122 for (k = 0 ; k < n ; k++)
123 {
124 Lnz [k] = 1 ;
125 }
126 }
127
128 ASSERT (is_ll == L->is_ll) ;
129 ASSERT (L->xtype != CHOLMOD_PATTERN) ;
130 DEBUG (CHOLMOD(dump_factor) (L, "L ready", Common)) ;
131 DEBUG (CHOLMOD(dump_sparse) (A, "A ready", Common)) ;
132 DEBUG (if (stype == 0) CHOLMOD(dump_sparse) (F, "F ready", Common)) ;
133
134 /* inputs, can be modified on output: */
135 Lp = L->p ; /* size n+1 */
136 ASSERT (Lp != NULL) ;
137
138 /* outputs, contents defined on input for incremental case only: */
139 Lnz = L->nz ; /* size n */
140 Lnext = L->next ; /* size n+2 */
141 Li = L->i ; /* size L->nzmax, can change in size */
142 Lx = L->x ; /* size L->nzmax or 2*L->nzmax, can change in size */
143 Lz = L->z ; /* size L->nzmax for zomplex case, can change in size */
144 nzmax = L->nzmax ;
145 ASSERT (Lnz != NULL && Li != NULL && Lx != NULL) ;
146
147 /* ---------------------------------------------------------------------- */
148 /* get workspace */
149 /* ---------------------------------------------------------------------- */
150
151 Iwork = Common->Iwork ;
152 Stack = Iwork ; /* size n (i/i/l) */
153 Flag = Common->Flag ; /* size n, Flag [i] < mark must hold */
154 Wx = Common->Xwork ; /* size n if real, 2*n if complex or
155 * zomplex. Xwork [i] == 0 must hold. */
156 Wz = Wx + n ; /* size n for zomplex case only */
157 mark = Common->mark ;
158 ASSERT ((Int) Common->xworksize >= (L->xtype == CHOLMOD_REAL ? 1:2)*n) ;
159
160 /* ---------------------------------------------------------------------- */
161 /* compute LDL' or LL' factorization by rows */
162 /* ---------------------------------------------------------------------- */
163
164 #ifdef MASK
165 #define NEXT(k) k = RLinkUp [k]
166 #else
167 #define NEXT(k) k++
168 #endif
169
170 for (k = kstart ; k < ((Int) kend) ; NEXT(k))
171 {
172 PRINT1 (("\n===============K "ID" Lnz [k] "ID"\n", k, Lnz [k])) ;
173
174 /* ------------------------------------------------------------------ */
175 /* compute pattern of kth row of L and scatter kth input column */
176 /* ------------------------------------------------------------------ */
177
178 /* column k of L is currently empty */
179 ASSERT (Lnz [k] == 1) ;
180 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, Common)) ;
181
182 top = n ; /* Stack is empty */
183 Flag [k] = mark ; /* do not include diagonal entry in Stack */
184
185 /* use Li [Lp [i]+1] for etree */
186 #define PARENT(i) (Lnz [i] > 1) ? (Li [Lp [i] + 1]) : EMPTY
187
188 if (stype > 0)
189 {
190 /* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
191 p = Ap [k] ;
192 pend = (packed) ? (Ap [k+1]) : (p + Anz [k]) ;
193 /* W [i] = Ax [i] ; scatter column of A */
194 #define SCATTER ASSIGN(Wx,Wz,i, Ax,Az,p)
195 SUBTREE ;
196 #undef SCATTER
197 }
198 else
199 {
200 /* scatter kth col of triu (beta*I+AA'), get pattern L(k,:) */
201 pf = Fp [k] ;
202 pfend = (Fpacked) ? (Fp [k+1]) : (pf + Fnz [k]) ;
203 for ( ; pf < pfend ; pf++)
204 {
205 /* get nonzero entry F (t,k) */
206 t = Fi [pf] ;
207 /* fk = Fx [pf] */
208 ASSIGN (fx, fz, 0, Fx, Fz, pf) ;
209 p = Ap [t] ;
210 pend = (packed) ? (Ap [t+1]) : (p + Anz [t]) ;
211 multadds = 0 ;
212 /* W [i] += Ax [p] * fx ; scatter column of A*A' */
213 #define SCATTER MULTADD (Wx,Wz,i, Ax,Az,p, fx,fz,0) ; multadds++ ;
214 SUBTREE ;
215 #undef SCATTER
216 #ifdef REAL
217 fl += 2 * ((double) multadds) ;
218 #else
219 fl += 8 * ((double) multadds) ;
220 #endif
221 }
222 }
223
224 #undef PARENT
225
226 /* ------------------------------------------------------------------ */
227 /* if mask is present, set the corresponding entries in W to zero */
228 /* ------------------------------------------------------------------ */
229
230 #ifdef MASK
231 /* remove the dead element of Wx */
232 if (mask != NULL)
233 {
234
235 #if 0
236 /* older version */
237 for (p = n; p > top;)
238 {
239 i = Stack [--p] ;
240 if ( mask [i] >= 0 )
241 {
242 CLEAR (Wx,Wz,i) ; /* set W(i) to zero */
243 }
244 }
245 #endif
246
247 for (s = top ; s < n ; s++)
248 {
249 i = Stack [s] ;
250 if (mask [i] >= maskmark)
251 {
252 CLEAR (Wx,Wz,i) ; /* set W(i) to zero */
253 }
254 }
255
256 }
257 #endif
258
259 /* nonzero pattern of kth row of L is now in Stack [top..n-1].
260 * Flag [Stack [top..n-1]] is equal to mark, but no longer needed */
261
262 /* mark = CHOLMOD(clear_flag) (Common) ; */
263 CHOLMOD_CLEAR_FLAG (Common) ;
264 mark = Common->mark ;
265
266 /* ------------------------------------------------------------------ */
267 /* compute kth row of L and store in column form */
268 /* ------------------------------------------------------------------ */
269
270 /* Solve L (0:k-1, 0:k-1) * y (0:k-1) = b (0:k-1) where
271 * b (0:k) = A (0:k,k) or A(0:k,:) * F(:,k) is in W and Stack.
272 *
273 * For LDL' factorization:
274 * L (k, 0:k-1) = y (0:k-1) ./ D (0:k-1)
275 * D (k) = b (k) - L (k, 0:k-1) * y (0:k-1)
276 *
277 * For LL' factorization:
278 * L (k, 0:k-1) = y (0:k-1)
279 * L (k,k) = sqrt (b (k) - L (k, 0:k-1) * L (0:k-1, k))
280 */
281
282 /* dk = W [k] + beta */
283 ADD_REAL (dk,0, Wx,k, beta,0) ;
284
285 #ifndef REAL
286 /* In the unsymmetric case, the imaginary part of W[k] must be real,
287 * since F is assumed to be the complex conjugate transpose of A. In
288 * the symmetric case, W[k] is the diagonal of A. If the imaginary part
289 * of W[k] is nonzero, then the Cholesky factorization cannot be
290 * computed; A is not positive definite */
291 dk_imaginary = (stype > 0) ? (IMAG_IS_NONZERO (Wx,Wz,k)) : FALSE ;
292 #endif
293
294 /* W [k] = 0.0 ; */
295 CLEAR (Wx,Wz,k) ;
296
297 for (s = top ; s < n ; s++)
298 {
299 /* get i for each nonzero entry L(k,i) */
300 i = Stack [s] ;
301
302 /* y = W [i] ; */
303 ASSIGN (yx,yz,0, Wx,Wz,i) ;
304
305 /* W [i] = 0.0 ; */
306 CLEAR (Wx,Wz,i) ;
307
308 lnz = Lnz [i] ;
309 p = Lp [i] ;
310 ASSERT (lnz > 0 && Li [p] == i) ;
311 pend = p + lnz ;
312
313 /* di = Lx [p] ; the diagonal entry L or D(i,i), which is real */
314 ASSIGN_REAL (di,0, Lx,p) ;
315
316 if (i >= (Int) L->minor || IS_ZERO (di [0]))
317 {
318 /* For the LL' factorization, L(i,i) is zero. For the LDL',
319 * D(i,i) is zero. Skip column i of L, and set L(k,i) = 0. */
320 CLEAR (lx,lz,0) ;
321 p = pend ;
322 }
323 else if (is_ll)
324 {
325 #ifdef REAL
326 fl += 2 * ((double) (pend - p - 1)) + 3 ;
327 #else
328 fl += 8 * ((double) (pend - p - 1)) + 6 ;
329 #endif
330 /* forward solve using L (i:(k-1),i) */
331 /* divide by L(i,i), which must be real and nonzero */
332 /* y /= di [0] */
333 DIV_REAL (yx,yz,0, yx,yz,0, di,0) ;
334 for (p++ ; p < pend ; p++)
335 {
336 /* W [Li [p]] -= Lx [p] * y ; */
337 MULTSUB (Wx,Wz,Li[p], Lx,Lz,p, yx,yz,0) ;
338 }
339 /* do not scale L; compute dot product for L(k,k) */
340 /* L(k,i) = conj(y) ; */
341 ASSIGN_CONJ (lx,lz,0, yx,yz,0) ;
342 /* d -= conj(y) * y ; */
343 LLDOT (dk,0, yx,yz,0) ;
344 }
345 else
346 {
347 #ifdef REAL
348 fl += 2 * ((double) (pend - p - 1)) + 3 ;
349 #else
350 fl += 8 * ((double) (pend - p - 1)) + 6 ;
351 #endif
352 /* forward solve using D (i,i) and L ((i+1):(k-1),i) */
353 for (p++ ; p < pend ; p++)
354 {
355 /* W [Li [p]] -= Lx [p] * y ; */
356 MULTSUB (Wx,Wz,Li[p], Lx,Lz,p, yx,yz,0) ;
357 }
358 /* Scale L (k,0:k-1) for LDL' factorization, compute D (k,k)*/
359 #ifdef REAL
360 /* L(k,i) = y/d */
361 lx [0] = yx [0] / di [0] ;
362 /* d -= L(k,i) * y */
363 dk [0] -= lx [0] * yx [0] ;
364 #else
365 /* L(k,i) = conj(y) ; */
366 ASSIGN_CONJ (lx,lz,0, yx,yz,0) ;
367 /* L(k,i) /= di ; */
368 DIV_REAL (lx,lz,0, lx,lz,0, di,0) ;
369 /* d -= conj(y) * y / di */
370 LDLDOT (dk,0, yx,yz,0, di,0) ;
371 #endif
372 }
373
374 /* determine if column i of L can hold the new L(k,i) entry */
375 if (p >= Lp [Lnext [i]])
376 {
377 /* column i needs to grow */
378 PRINT1 (("Factor Colrealloc "ID", old Lnz "ID"\n", i, Lnz [i]));
379 if (!CHOLMOD(reallocate_column) (i, lnz + 1, L, Common))
380 {
381 /* out of memory, L is now simplicial symbolic */
382 for (i = 0 ; i < n ; i++)
383 {
384 /* W [i] = 0 ; */
385 CLEAR (Wx,Wz,i) ;
386 }
387 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, n, Common)) ;
388 return (FALSE) ;
389 }
390 Li = L->i ; /* L->i, L->x, L->z may have moved */
391 Lx = L->x ;
392 Lz = L->z ;
393 p = Lp [i] + lnz ; /* contents of L->p changed */
394 ASSERT (p < Lp [Lnext [i]]) ;
395 }
396
397 /* store L (k,i) in the column form matrix of L */
398 Li [p] = k ;
399 /* Lx [p] = L(k,i) ; */
400 ASSIGN (Lx,Lz,p, lx,lz,0) ;
401 Lnz [i]++ ;
402 }
403
404 /* ------------------------------------------------------------------ */
405 /* ensure abs (d) >= dbound if dbound is given, and store it in L */
406 /* ------------------------------------------------------------------ */
407
408 p = Lp [k] ;
409 Li [p] = k ;
410
411 if (k >= (Int) L->minor)
412 {
413 /* the matrix is already not positive definite */
414 dk [0] = 0 ;
415 }
416 else if (use_dbound)
417 {
418 /* modify the diagonal to force LL' or LDL' to exist */
419 dk [0] = CHOLMOD(dbound) (is_ll ? fabs (dk [0]) : dk [0], Common) ;
420 }
421 else if ((is_ll ? (IS_LE_ZERO (dk [0])) : (IS_ZERO (dk [0])))
422 #ifndef REAL
423 || dk_imaginary
424 #endif
425 )
426 {
427 /* the matrix has just been found to be not positive definite */
428 dk [0] = 0 ;
429 L->minor = k ;
430 ERROR (CHOLMOD_NOT_POSDEF, "not positive definite") ;
431 }
432
433 if (is_ll)
434 {
435 /* this is counted as one flop, below */
436 dk [0] = sqrt (dk [0]) ;
437 }
438
439 /* Lx [p] = D(k,k) = d ; real part only */
440 ASSIGN_REAL (Lx,p, dk,0) ;
441 CLEAR_IMAG (Lx,Lz,p) ;
442 }
443
444 #undef NEXT
445
446 if (is_ll) fl += MAX ((Int) kend - (Int) kstart, 0) ; /* count sqrt's */
447 Common->rowfacfl = fl ;
448
449 DEBUG (CHOLMOD(dump_factor) (L, "final cholmod_rowfac", Common)) ;
450 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, n, Common)) ;
451 return (TRUE) ;
452 }
453 #undef PATTERN
454 #undef REAL
455 #undef COMPLEX
456 #undef ZOMPLEX
457