1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
|
//------------------------------------------------------------------------------
// CHOLMOD/Modify/cholmod_rowdel: delete row/column from an LDL' factorization
//------------------------------------------------------------------------------
// CHOLMOD/Modify Module. Copyright (C) 2005-2023, Timothy A. Davis,
// and William W. Hager. All Rights Reserved.
// SPDX-License-Identifier: GPL-2.0+
//------------------------------------------------------------------------------
// Deletes a row and column from an LDL' factorization. The row and column k
// is set to the kth row and column of the identity matrix. Optionally
// downdates the solution to Lx=b.
//
// workspace: Flag (nrow), Head (nrow+1), W (nrow*2), Iwork (2*nrow)
//
// Only real matrices (single or double) are supported (exception: since only
// the pattern of R is used, it can have any valid xtype).
#include "cholmod_internal.h"
#ifndef NGPL
#ifndef NMODIFY
//------------------------------------------------------------------------------
// t_cholmod_rowdel_worker
//------------------------------------------------------------------------------
#define DOUBLE
#define REAL
#include "t_cholmod_rowdel_worker.c"
#undef DOUBLE
#define SINGLE
#define REAL
#include "t_cholmod_rowdel_worker.c"
//------------------------------------------------------------------------------
// cholmod_rowdel
//------------------------------------------------------------------------------
// Sets the kth row and column of L to be the kth row and column of the identity
// matrix, and updates L(k+1:n,k+1:n) accordingly. To reduce the running time,
// the caller can optionally provide the nonzero pattern (or an upper bound) of
// kth row of L, as the sparse n-by-1 vector R. Provide R as NULL if you want
// CHOLMOD to determine this itself, which is easier for the caller, but takes
// a little more time.
int CHOLMOD(rowdel)
(
// input:
size_t k, // row/column index to delete
cholmod_sparse *R, // NULL, or the nonzero pattern of kth row of L
// input/output:
cholmod_factor *L, // factor to modify
cholmod_common *Common
)
{
double yk [2] ;
yk [0] = 0. ;
yk [1] = 0. ;
return (CHOLMOD(rowdel_mark) (k, R, yk, NULL, L, NULL, NULL, Common)) ;
}
//------------------------------------------------------------------------------
// cholmod_rowdel_solve
//------------------------------------------------------------------------------
// Does the same as cholmod_rowdel, but also downdates the solution to Lx=b.
// When row/column k of A is "deleted" from the system A*y=b, this can induce
// a change to x, in addition to changes arising when L and b are modified.
// If this is the case, the kth entry of y is required as input (yk).
int CHOLMOD(rowdel_solve)
(
// input:
size_t k, // row/column index to delete
cholmod_sparse *R, // NULL, or the nonzero pattern of kth row of L
double yk [2], // kth entry in the solution to A*y=b
// input/output:
cholmod_factor *L, // factor to modify
cholmod_dense *X, // solution to Lx=b (size n-by-1)
cholmod_dense *DeltaB, // change in b, zero on output
cholmod_common *Common
)
{
return (CHOLMOD(rowdel_mark) (k, R, yk, NULL, L, X, DeltaB, Common)) ;
}
//------------------------------------------------------------------------------
// cholmod_rowdel_mark
//------------------------------------------------------------------------------
// Does the same as cholmod_rowdel_solve, except only part of L is used in
// the update/downdate of the solution to Lx=b. This routine is an "expert"
// routine. It is meant for use in LPDASA only.
//
// if R == NULL then columns 0:k-1 of L are searched for row k. Otherwise, it
// searches columns in the set defined by the pattern of the first column of R.
// This is meant to be the pattern of row k of L (a superset of that pattern is
// OK too). R must be a permutation of a subset of 0:k-1.
int CHOLMOD(rowdel_mark)
(
// input:
size_t kdel, // row/column index to delete
cholmod_sparse *R, // NULL, or the nonzero pattern of kth row of L
double yk [2], // kth entry in the solution to A*y=b
Int *colmark, // Int array of size 1. See cholmod_updown.c
// input/output:
cholmod_factor *L, // factor to modify
cholmod_dense *X, // solution to Lx=b (size n-by-1)
cholmod_dense *DeltaB, // change in b, zero on output
cholmod_common *Common
)
{
//--------------------------------------------------------------------------
// check inputs
//--------------------------------------------------------------------------
RETURN_IF_NULL_COMMON (FALSE) ;
RETURN_IF_NULL (L, FALSE) ;
RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_REAL, FALSE) ;
Int k = kdel ;
if (kdel >= L->n || k < 0)
{
ERROR (CHOLMOD_INVALID, "k invalid") ;
return (FALSE) ;
}
if (R != NULL)
{
RETURN_IF_XTYPE_INVALID (R, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
if (R->ncol != 1 || R->nrow != L->n)
{
ERROR (CHOLMOD_INVALID, "R invalid") ;
return (FALSE) ;
}
}
if ((X != NULL) && (DeltaB != NULL))
{
RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
RETURN_IF_XTYPE_INVALID (DeltaB, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
if (X->nrow != L->n || X->ncol != 1 ||
DeltaB->nrow != L->n || DeltaB->ncol != 1 ||
X->dtype != L->dtype || DeltaB->dtype != L->dtype)
{
ERROR (CHOLMOD_INVALID, "X and/or DeltaB invalid") ;
return (FALSE) ;
}
}
Common->status = CHOLMOD_OK ;
//--------------------------------------------------------------------------
// allocate workspace
//--------------------------------------------------------------------------
// s = 2*n
int ok = TRUE ;
size_t s = CHOLMOD(mult_size_t) (L->n, 2, &ok) ;
if (!ok)
{
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
return (FALSE) ;
}
CHOLMOD(alloc_work) (L->n, s, s, L->dtype, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (FALSE) ;
}
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, L->dtype, Common)) ;
//--------------------------------------------------------------------------
// convert to simplicial numeric LDL' factor, if not already
//--------------------------------------------------------------------------
if (L->xtype == CHOLMOD_PATTERN || L->is_super || L->is_ll)
{
// can only update/downdate a simplicial LDL' factorization
CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, FALSE, FALSE, L,
Common) ;
if (Common->status < CHOLMOD_OK)
{
// out of memory, L is returned unchanged
return (FALSE) ;
}
}
//--------------------------------------------------------------------------
// update L and X
//--------------------------------------------------------------------------
float s_yk [2] ;
s_yk [0] = (float) yk [0] ;
s_yk [1] = (float) yk [1] ;
switch (L->dtype & 4)
{
case CHOLMOD_SINGLE:
ok = rs_cholmod_rowdel_worker (k, R, s_yk, colmark, L, X, DeltaB,
Common) ;
break ;
case CHOLMOD_DOUBLE:
ok = rd_cholmod_rowdel_worker (k, R, yk, colmark, L, X, DeltaB,
Common) ;
break ;
}
//--------------------------------------------------------------------------
// return result
//--------------------------------------------------------------------------
DEBUG (CHOLMOD(dump_factor) (L, "LDL factorization, L:", Common)) ;
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, s, L->dtype, Common)) ;
return (ok) ;
}
#endif
#endif
|
