1 /* ========================================================================== */
2 /* === Modify/cholmod_rowdel ================================================ */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Modify Module.
7 * Copyright (C) 2005-2006, Timothy A. Davis and William W. Hager.
8 * http://www.suitesparse.com
9 * -------------------------------------------------------------------------- */
10
11 /* Deletes a row and column from an LDL' factorization. The row and column k
12 * is set to the kth row and column of the identity matrix. Optionally
13 * downdates the solution to Lx=b.
14 *
15 * workspace: Flag (nrow), Head (nrow+1), W (nrow*2), Iwork (2*nrow)
16 *
17 * Only real matrices are supported (exception: since only the pattern of R
18 * is used, it can have any valid xtype).
19 */
20
21 #ifndef NGPL
22 #ifndef NMODIFY
23
24 #include "cholmod_internal.h"
25 #include "cholmod_modify.h"
26
27
28 /* ========================================================================== */
29 /* === cholmod_rowdel ======================================================= */
30 /* ========================================================================== */
31
32 /* Sets the kth row and column of L to be the kth row and column of the identity
33 * matrix, and updates L(k+1:n,k+1:n) accordingly. To reduce the running time,
34 * the caller can optionally provide the nonzero pattern (or an upper bound) of
35 * kth row of L, as the sparse n-by-1 vector R. Provide R as NULL if you want
36 * CHOLMOD to determine this itself, which is easier for the caller, but takes
37 * a little more time.
38 */
39
CHOLMOD(rowdel)40 int CHOLMOD(rowdel)
41 (
42 /* ---- input ---- */
43 size_t k, /* row/column index to delete */
44 cholmod_sparse *R, /* NULL, or the nonzero pattern of kth row of L */
45 /* ---- in/out --- */
46 cholmod_factor *L, /* factor to modify */
47 /* --------------- */
48 cholmod_common *Common
49 )
50 {
51 double yk [2] ;
52 yk [0] = 0. ;
53 yk [1] = 0. ;
54 return (CHOLMOD(rowdel_mark) (k, R, yk, NULL, L, NULL, NULL, Common)) ;
55 }
56
57
58 /* ========================================================================== */
59 /* === cholmod_rowdel_solve ================================================= */
60 /* ========================================================================== */
61
62 /* Does the same as cholmod_rowdel, but also downdates the solution to Lx=b.
63 * When row/column k of A is "deleted" from the system A*y=b, this can induce
64 * a change to x, in addition to changes arising when L and b are modified.
65 * If this is the case, the kth entry of y is required as input (yk) */
66
CHOLMOD(rowdel_solve)67 int CHOLMOD(rowdel_solve)
68 (
69 /* ---- input ---- */
70 size_t k, /* row/column index to delete */
71 cholmod_sparse *R, /* NULL, or the nonzero pattern of kth row of L */
72 double yk [2], /* kth entry in the solution to A*y=b */
73 /* ---- in/out --- */
74 cholmod_factor *L, /* factor to modify */
75 cholmod_dense *X, /* solution to Lx=b (size n-by-1) */
76 cholmod_dense *DeltaB, /* change in b, zero on output */
77 /* --------------- */
78 cholmod_common *Common
79 )
80 {
81 return (CHOLMOD(rowdel_mark) (k, R, yk, NULL, L, X, DeltaB, Common)) ;
82 }
83
84
85 /* ========================================================================== */
86 /* === cholmod_rowdel_mark ================================================== */
87 /* ========================================================================== */
88
89 /* Does the same as cholmod_rowdel_solve, except only part of L is used in
90 * the update/downdate of the solution to Lx=b. This routine is an "expert"
91 * routine. It is meant for use in LPDASA only.
92 *
93 * if R == NULL then columns 0:k-1 of L are searched for row k. Otherwise, it
94 * searches columns in the set defined by the pattern of the first column of R.
95 * This is meant to be the pattern of row k of L (a superset of that pattern is
96 * OK too). R must be a permutation of a subset of 0:k-1.
97 */
98
CHOLMOD(rowdel_mark)99 int CHOLMOD(rowdel_mark)
100 (
101 /* ---- input ---- */
102 size_t kdel, /* row/column index to delete */
103 cholmod_sparse *R, /* NULL, or the nonzero pattern of kth row of L */
104 double yk [2], /* kth entry in the solution to A*y=b */
105 Int *colmark, /* Int array of size 1. See cholmod_updown.c */
106 /* ---- in/out --- */
107 cholmod_factor *L, /* factor to modify */
108 cholmod_dense *X, /* solution to Lx=b (size n-by-1) */
109 cholmod_dense *DeltaB, /* change in b, zero on output */
110 /* --------------- */
111 cholmod_common *Common
112 )
113 {
114 double dk, sqrt_dk, xk, dj, fl ;
115 double *Lx, *Cx, *W, *Xx, *Nx ;
116 Int *Li, *Lp, *Lnz, *Ci, *Rj, *Rp, *Iwork ;
117 cholmod_sparse *C, Cmatrix ;
118 Int j, p, pend, kk, lnz, n, Cp [2], do_solve, do_update, left, k,
119 right, middle, i, klast, given_row, rnz ;
120 size_t s ;
121 int ok = TRUE ;
122
123 /* ---------------------------------------------------------------------- */
124 /* check inputs */
125 /* ---------------------------------------------------------------------- */
126
127 RETURN_IF_NULL_COMMON (FALSE) ;
128 RETURN_IF_NULL (L, FALSE) ;
129 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_REAL, FALSE) ;
130 n = L->n ;
131 k = kdel ;
132 if (kdel >= L->n || k < 0)
133 {
134 ERROR (CHOLMOD_INVALID, "k invalid") ;
135 return (FALSE) ;
136 }
137 if (R == NULL)
138 {
139 Rj = NULL ;
140 rnz = EMPTY ;
141 }
142 else
143 {
144 RETURN_IF_XTYPE_INVALID (R, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
145 if (R->ncol != 1 || R->nrow != L->n)
146 {
147 ERROR (CHOLMOD_INVALID, "R invalid") ;
148 return (FALSE) ;
149 }
150 Rj = R->i ;
151 Rp = R->p ;
152 rnz = Rp [1] ;
153 }
154 do_solve = (X != NULL) && (DeltaB != NULL) ;
155 if (do_solve)
156 {
157 RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
158 RETURN_IF_XTYPE_INVALID (DeltaB, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
159 Xx = X->x ;
160 Nx = DeltaB->x ;
161 if (X->nrow != L->n || X->ncol != 1 || DeltaB->nrow != L->n ||
162 DeltaB->ncol != 1 || Xx == NULL || Nx == NULL)
163 {
164 ERROR (CHOLMOD_INVALID, "X and/or DeltaB invalid") ;
165 return (FALSE) ;
166 }
167 }
168 else
169 {
170 Xx = NULL ;
171 Nx = NULL ;
172 }
173 Common->status = CHOLMOD_OK ;
174
175 /* ---------------------------------------------------------------------- */
176 /* allocate workspace */
177 /* ---------------------------------------------------------------------- */
178
179 /* s = 2*n */
180 s = CHOLMOD(mult_size_t) (n, 2, &ok) ;
181 if (!ok)
182 {
183 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
184 return (FALSE) ;
185 }
186
187 CHOLMOD(allocate_work) (n, s, s, Common) ;
188 if (Common->status < CHOLMOD_OK)
189 {
190 return (FALSE) ;
191 }
192 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, Common)) ;
193
194 /* ---------------------------------------------------------------------- */
195 /* convert to simplicial numeric LDL' factor, if not already */
196 /* ---------------------------------------------------------------------- */
197
198 if (L->xtype == CHOLMOD_PATTERN || L->is_super || L->is_ll)
199 {
200 /* can only update/downdate a simplicial LDL' factorization */
201 CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, FALSE, FALSE, L,
202 Common) ;
203 if (Common->status < CHOLMOD_OK)
204 {
205 /* out of memory, L is returned unchanged */
206 return (FALSE) ;
207 }
208 }
209
210 /* ---------------------------------------------------------------------- */
211 /* get inputs */
212 /* ---------------------------------------------------------------------- */
213
214 /* inputs, not modified on output: */
215 Lp = L->p ; /* size n+1 */
216
217 /* outputs, contents defined on input for incremental case only: */
218 Lnz = L->nz ; /* size n */
219 Li = L->i ; /* size L->nzmax. Can change in size. */
220 Lx = L->x ; /* size L->nzmax. Can change in size. */
221
222 ASSERT (L->nz != NULL) ;
223
224 /* ---------------------------------------------------------------------- */
225 /* get workspace */
226 /* ---------------------------------------------------------------------- */
227
228 W = Common->Xwork ; /* size n, used only in cholmod_updown */
229 Cx = W + n ; /* use 2nd column of Xwork for C (size n) */
230 Iwork = Common->Iwork ;
231 Ci = Iwork + n ; /* size n (i/i/l) */
232 /* NOTE: cholmod_updown uses Iwork [0..n-1] (i/i/l) as Stack */
233
234 /* ---------------------------------------------------------------------- */
235 /* prune row k from all columns of L */
236 /* ---------------------------------------------------------------------- */
237
238 given_row = (rnz >= 0) ;
239 klast = given_row ? rnz : k ;
240 PRINT2 (("given_row "ID"\n", given_row)) ;
241
242 for (kk = 0 ; kk < klast ; kk++)
243 {
244 /* either search j = 0:k-1 or j = Rj [0:rnz-1] */
245 j = given_row ? (Rj [kk]) : (kk) ;
246
247 if (j < 0 || j >= k)
248 {
249 ERROR (CHOLMOD_INVALID, "R invalid") ;
250 return (FALSE) ;
251 }
252
253 PRINT2 (("Prune col j = "ID":\n", j)) ;
254
255 lnz = Lnz [j] ;
256 dj = Lx [Lp [j]] ;
257 ASSERT (Lnz [j] > 0 && Li [Lp [j]] == j) ;
258
259 if (lnz > 1)
260 {
261 left = Lp [j] ;
262 pend = left + lnz ;
263 right = pend - 1 ;
264
265 i = Li [right] ;
266
267 if (i < k)
268 {
269 /* row k is not in column j */
270 continue ;
271 }
272 else if (i == k)
273 {
274 /* k is the last row index in this column (quick delete) */
275 if (do_solve)
276 {
277 Xx [j] -= yk [0] * dj * Lx [right] ;
278 }
279 Lx [right] = 0 ;
280 }
281 else
282 {
283 /* binary search for row k in column j */
284 PRINT2 (("\nBinary search: lnz "ID" k = "ID"\n", lnz, k)) ;
285 while (left < right)
286 {
287 middle = (left + right) / 2 ;
288 PRINT2 (("left "ID" right "ID" middle "ID": ["ID" "ID""
289 ""ID"]\n", left, right, middle,
290 Li [left], Li [middle], Li [right])) ;
291 if (k > Li [middle])
292 {
293 left = middle + 1 ;
294 }
295 else
296 {
297 right = middle ;
298 }
299 }
300 ASSERT (left >= Lp [j] && left < pend) ;
301
302 #ifndef NDEBUG
303 /* brute force, linear-time search */
304 {
305 Int p3 = Lp [j] ;
306 i = EMPTY ;
307 PRINT2 (("Brute force:\n")) ;
308 for ( ; p3 < pend ; p3++)
309 {
310 i = Li [p3] ;
311 PRINT2 (("p "ID" ["ID"]\n", p3, i)) ;
312 if (i >= k)
313 {
314 break ;
315 }
316 }
317 if (i == k)
318 {
319 ASSERT (k == Li [p3]) ;
320 ASSERT (p3 == left) ;
321 }
322 }
323 #endif
324
325 if (k == Li [left])
326 {
327 if (do_solve)
328 {
329 Xx [j] -= yk [0] * dj * Lx [left] ;
330 }
331 /* found row k in column j. Prune it from the column.*/
332 Lx [left] = 0 ;
333 }
334 }
335 }
336 }
337
338 #ifndef NDEBUG
339 /* ensure that row k has been deleted from the matrix L */
340 for (j = 0 ; j < k ; j++)
341 {
342 Int lasti ;
343 lasti = EMPTY ;
344 p = Lp [j] ;
345 pend = p + Lnz [j] ;
346 /* look for row k in column j */
347 PRINT1 (("Pruned column "ID"\n", j)) ;
348 for ( ; p < pend ; p++)
349 {
350 i = Li [p] ;
351 PRINT2 ((" "ID"", i)) ;
352 PRINT2 ((" %g\n", Lx [p])) ;
353 ASSERT (IMPLIES (i == k, Lx [p] == 0)) ;
354 ASSERT (i > lasti) ;
355 lasti = i ;
356 }
357 PRINT1 (("\n")) ;
358 }
359 #endif
360
361 /* ---------------------------------------------------------------------- */
362 /* set diagonal and clear column k of L */
363 /* ---------------------------------------------------------------------- */
364
365 lnz = Lnz [k] - 1 ;
366 ASSERT (Lnz [k] > 0) ;
367
368 /* ---------------------------------------------------------------------- */
369 /* update/downdate */
370 /* ---------------------------------------------------------------------- */
371
372 /* update or downdate L (k+1:n, k+1:n) with the vector
373 * C = L (:,k) * sqrt (abs (D [k]))
374 * Do a numeric update if D[k] > 0, numeric downdate otherwise.
375 */
376
377 PRINT1 (("rowdel downdate lnz = "ID"\n", lnz)) ;
378
379 /* store the new unit diagonal */
380 p = Lp [k] ;
381 pend = p + lnz + 1 ;
382 dk = Lx [p] ;
383 Lx [p++] = 1 ;
384 PRINT2 (("D [k = "ID"] = %g\n", k, dk)) ;
385 ok = TRUE ;
386 fl = 0 ;
387
388 if (lnz > 0)
389 {
390 /* compute DeltaB for updown (in DeltaB) */
391 if (do_solve)
392 {
393 xk = Xx [k] - yk [0] * dk ;
394 for ( ; p < pend ; p++)
395 {
396 Nx [Li [p]] += Lx [p] * xk ;
397 }
398 }
399
400 do_update = IS_GT_ZERO (dk) ;
401 if (!do_update)
402 {
403 dk = -dk ;
404 }
405 sqrt_dk = sqrt (dk) ;
406 p = Lp [k] + 1 ;
407 for (kk = 0 ; kk < lnz ; kk++, p++)
408 {
409 Ci [kk] = Li [p] ;
410 Cx [kk] = Lx [p] * sqrt_dk ;
411 Lx [p] = 0 ; /* clear column k */
412 }
413 fl = lnz + 1 ;
414
415 /* create a n-by-1 sparse matrix to hold the single column */
416 C = &Cmatrix ;
417 C->nrow = n ;
418 C->ncol = 1 ;
419 C->nzmax = lnz ;
420 C->sorted = TRUE ;
421 C->packed = TRUE ;
422 C->p = Cp ;
423 C->i = Ci ;
424 C->x = Cx ;
425 C->nz = NULL ;
426 C->itype = L->itype ;
427 C->xtype = L->xtype ;
428 C->dtype = L->dtype ;
429 C->z = NULL ;
430 C->stype = 0 ;
431
432 Cp [0] = 0 ;
433 Cp [1] = lnz ;
434
435 /* numeric update if dk > 0, and with Lx=b change */
436 /* workspace: Flag (nrow), Head (nrow+1), W (nrow), Iwork (2*nrow) */
437 ok = CHOLMOD(updown_mark) (do_update ? (1) : (0), C, colmark,
438 L, X, DeltaB, Common) ;
439
440 /* clear workspace */
441 for (kk = 0 ; kk < lnz ; kk++)
442 {
443 Cx [kk] = 0 ;
444 }
445 }
446
447 Common->modfl += fl ;
448
449 if (do_solve)
450 {
451 /* kth equation becomes identity, so X(k) is now Y(k) */
452 Xx [k] = yk [0] ;
453 }
454
455 DEBUG (CHOLMOD(dump_factor) (L, "LDL factorization, L:", Common)) ;
456 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, Common)) ;
457 return (ok) ;
458 }
459 #endif
460 #endif
461