1 #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h> /*I "petscksp.h" I*/
2
3 /*------------------------------------------------------------*/
4
5 /*
6 The solution method is the matrix-free implementation of the inverse Hessian
7 representation in page 312 of Griewank "Broyden Updating, The Good and The Bad!"
8 (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).
9
10 Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
11 the matrix is updated with a new (S[i], Y[i]) pair. This allows
12 repeated calls of MatSolve without incurring redundant computation.
13
14 dX <- J0^{-1} * F
15
16 for i=0,1,2,...,k
17 # Q[i] = (B_i)^{-1} * Y[i]
18 tau = (S[i]^T dX) / (S[i]^T Q[i])
19 dX <- dX + (tau * (S[i] - Q[i]))
20 end
21 */
22
MatSolve_LMVMBrdn(Mat B,Vec F,Vec dX)23 static PetscErrorCode MatSolve_LMVMBrdn(Mat B, Vec F, Vec dX)
24 {
25 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
26 Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
27 PetscErrorCode ierr;
28 PetscInt i, j;
29 PetscScalar sjtqi, stx, stq;
30
31 PetscFunctionBegin;
32 VecCheckSameSize(F, 2, dX, 3);
33 VecCheckMatCompatible(B, dX, 3, F, 2);
34
35 if (lbrdn->needQ) {
36 /* Pre-compute (Q[i] = (B_i)^{-1} * Y[i]) */
37 for (i = 0; i <= lmvm->k; ++i) {
38 ierr = MatLMVMApplyJ0Inv(B, lmvm->Y[i], lbrdn->Q[i]);CHKERRQ(ierr);
39 for (j = 0; j <= i-1; ++j) {
40 ierr = VecDot(lmvm->S[j], lbrdn->Q[i], &sjtqi);CHKERRQ(ierr);
41 ierr = VecAXPBYPCZ(lbrdn->Q[i], PetscRealPart(sjtqi)/lbrdn->stq[j], -PetscRealPart(sjtqi)/lbrdn->stq[j], 1.0, lmvm->S[j], lbrdn->Q[j]);CHKERRQ(ierr);
42 }
43 ierr = VecDot(lmvm->S[i], lbrdn->Q[i], &stq);CHKERRQ(ierr);
44 lbrdn->stq[i] = PetscRealPart(stq);
45 }
46 lbrdn->needQ = PETSC_FALSE;
47 }
48
49 ierr = MatLMVMApplyJ0Inv(B, F, dX);CHKERRQ(ierr);
50 for (i = 0; i <= lmvm->k; ++i) {
51 ierr = VecDot(lmvm->S[i], dX, &stx);CHKERRQ(ierr);
52 ierr = VecAXPBYPCZ(dX, PetscRealPart(stx)/lbrdn->stq[i], -PetscRealPart(stx)/lbrdn->stq[i], 1.0, lmvm->S[i], lbrdn->Q[i]);CHKERRQ(ierr);
53 }
54 PetscFunctionReturn(0);
55 }
56
57 /*------------------------------------------------------------*/
58
59 /*
60 The forward product is the matrix-free implementation of Equation 2 in
61 page 302 of Griewank "Broyden Updating, The Good and The Bad!"
62 (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).
63
64 P[i] = (B_i)*S[i] terms are computed ahead of time whenever
65 the matrix is updated with a new (S[i], Y[i]) pair. This allows
66 repeated calls of MatMult inside KSP solvers without unnecessarily
67 recomputing P[i] terms in expensive nested-loops.
68
69 Z <- J0 * X
70
71 for i=0,1,2,...,k
72 # P[i] = B_i * S[i]
73 tau = (S[i]^T X) / (S[i]^T S[i])
74 dX <- dX + (tau * (Y[i] - P[i]))
75 end
76 */
77
MatMult_LMVMBrdn(Mat B,Vec X,Vec Z)78 static PetscErrorCode MatMult_LMVMBrdn(Mat B, Vec X, Vec Z)
79 {
80 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
81 Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
82 PetscErrorCode ierr;
83 PetscInt i, j;
84 PetscScalar sjtsi, stx;
85
86 PetscFunctionBegin;
87 VecCheckSameSize(X, 2, Z, 3);
88 VecCheckMatCompatible(B, X, 2, Z, 3);
89
90 if (lbrdn->needP) {
91 /* Pre-compute (P[i] = (B_i) * S[i]) */
92 for (i = 0; i <= lmvm->k; ++i) {
93 ierr = MatLMVMApplyJ0Fwd(B, lmvm->S[i], lbrdn->P[i]);CHKERRQ(ierr);
94 for (j = 0; j <= i-1; ++j) {
95 ierr = VecDot(lmvm->S[j], lmvm->S[i], &sjtsi);CHKERRQ(ierr);
96 ierr = VecAXPBYPCZ(lbrdn->P[i], PetscRealPart(sjtsi)/lbrdn->sts[j], -PetscRealPart(sjtsi)/lbrdn->sts[j], 1.0, lmvm->Y[j], lbrdn->P[j]);CHKERRQ(ierr);
97 }
98 }
99 lbrdn->needP = PETSC_FALSE;
100 }
101
102 ierr = MatLMVMApplyJ0Fwd(B, X, Z);CHKERRQ(ierr);
103 for (i = 0; i <= lmvm->k; ++i) {
104 ierr = VecDot(lmvm->S[i], X, &stx);CHKERRQ(ierr);
105 ierr = VecAXPBYPCZ(Z, PetscRealPart(stx)/lbrdn->sts[i], -PetscRealPart(stx)/lbrdn->sts[i], 1.0, lmvm->Y[i], lbrdn->P[i]);CHKERRQ(ierr);
106 }
107 PetscFunctionReturn(0);
108 }
109
110 /*------------------------------------------------------------*/
111
MatUpdate_LMVMBrdn(Mat B,Vec X,Vec F)112 static PetscErrorCode MatUpdate_LMVMBrdn(Mat B, Vec X, Vec F)
113 {
114 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
115 Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
116 PetscErrorCode ierr;
117 PetscInt old_k, i;
118 PetscScalar sts;
119
120 PetscFunctionBegin;
121 if (!lmvm->m) PetscFunctionReturn(0);
122 if (lmvm->prev_set) {
123 /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
124 ierr = VecAYPX(lmvm->Xprev, -1.0, X);CHKERRQ(ierr);
125 ierr = VecAYPX(lmvm->Fprev, -1.0, F);CHKERRQ(ierr);
126 /* Accept the update */
127 lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
128 old_k = lmvm->k;
129 ierr = MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);CHKERRQ(ierr);
130 /* If we hit the memory limit, shift the sts array */
131 if (old_k == lmvm->k) {
132 for (i = 0; i <= lmvm->k-1; ++i) {
133 lbrdn->sts[i] = lbrdn->sts[i+1];
134 }
135 }
136 ierr = VecDot(lmvm->S[lmvm->k], lmvm->S[lmvm->k], &sts);CHKERRQ(ierr);
137 lbrdn->sts[lmvm->k] = PetscRealPart(sts);
138 }
139 /* Save the solution and function to be used in the next update */
140 ierr = VecCopy(X, lmvm->Xprev);CHKERRQ(ierr);
141 ierr = VecCopy(F, lmvm->Fprev);CHKERRQ(ierr);
142 lmvm->prev_set = PETSC_TRUE;
143 PetscFunctionReturn(0);
144 }
145
146 /*------------------------------------------------------------*/
147
MatCopy_LMVMBrdn(Mat B,Mat M,MatStructure str)148 static PetscErrorCode MatCopy_LMVMBrdn(Mat B, Mat M, MatStructure str)
149 {
150 Mat_LMVM *bdata = (Mat_LMVM*)B->data;
151 Mat_Brdn *bctx = (Mat_Brdn*)bdata->ctx;
152 Mat_LMVM *mdata = (Mat_LMVM*)M->data;
153 Mat_Brdn *mctx = (Mat_Brdn*)mdata->ctx;
154 PetscErrorCode ierr;
155 PetscInt i;
156
157 PetscFunctionBegin;
158 mctx->needP = bctx->needP;
159 mctx->needQ = bctx->needQ;
160 for (i=0; i<=bdata->k; ++i) {
161 mctx->sts[i] = bctx->sts[i];
162 mctx->stq[i] = bctx->stq[i];
163 ierr = VecCopy(bctx->P[i], mctx->P[i]);CHKERRQ(ierr);
164 ierr = VecCopy(bctx->Q[i], mctx->Q[i]);CHKERRQ(ierr);
165 }
166 PetscFunctionReturn(0);
167 }
168
169 /*------------------------------------------------------------*/
170
MatReset_LMVMBrdn(Mat B,PetscBool destructive)171 static PetscErrorCode MatReset_LMVMBrdn(Mat B, PetscBool destructive)
172 {
173 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
174 Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
175 PetscErrorCode ierr;
176
177 PetscFunctionBegin;
178 lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
179 if (destructive && lbrdn->allocated) {
180 ierr = PetscFree2(lbrdn->sts, lbrdn->stq);CHKERRQ(ierr);
181 ierr = VecDestroyVecs(lmvm->m, &lbrdn->P);CHKERRQ(ierr);
182 ierr = VecDestroyVecs(lmvm->m, &lbrdn->Q);CHKERRQ(ierr);
183 lbrdn->allocated = PETSC_FALSE;
184 }
185 ierr = MatReset_LMVM(B, destructive);CHKERRQ(ierr);
186 PetscFunctionReturn(0);
187 }
188
189 /*------------------------------------------------------------*/
190
MatAllocate_LMVMBrdn(Mat B,Vec X,Vec F)191 static PetscErrorCode MatAllocate_LMVMBrdn(Mat B, Vec X, Vec F)
192 {
193 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
194 Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
195 PetscErrorCode ierr;
196
197 PetscFunctionBegin;
198 ierr = MatAllocate_LMVM(B, X, F);CHKERRQ(ierr);
199 if (!lbrdn->allocated) {
200 ierr = PetscMalloc2(lmvm->m, &lbrdn->sts, lmvm->m, &lbrdn->stq);CHKERRQ(ierr);
201 if (lmvm->m > 0) {
202 ierr = VecDuplicateVecs(X, lmvm->m, &lbrdn->P);CHKERRQ(ierr);
203 ierr = VecDuplicateVecs(X, lmvm->m, &lbrdn->Q);CHKERRQ(ierr);
204 }
205 lbrdn->allocated = PETSC_TRUE;
206 }
207 PetscFunctionReturn(0);
208 }
209
210 /*------------------------------------------------------------*/
211
MatDestroy_LMVMBrdn(Mat B)212 static PetscErrorCode MatDestroy_LMVMBrdn(Mat B)
213 {
214 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
215 Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
216 PetscErrorCode ierr;
217
218 PetscFunctionBegin;
219 if (lbrdn->allocated) {
220 ierr = PetscFree2(lbrdn->sts, lbrdn->stq);CHKERRQ(ierr);
221 ierr = VecDestroyVecs(lmvm->m, &lbrdn->P);CHKERRQ(ierr);
222 ierr = VecDestroyVecs(lmvm->m, &lbrdn->Q);CHKERRQ(ierr);
223 lbrdn->allocated = PETSC_FALSE;
224 }
225 ierr = PetscFree(lmvm->ctx);CHKERRQ(ierr);
226 ierr = MatDestroy_LMVM(B);CHKERRQ(ierr);
227 PetscFunctionReturn(0);
228 }
229
230 /*------------------------------------------------------------*/
231
MatSetUp_LMVMBrdn(Mat B)232 static PetscErrorCode MatSetUp_LMVMBrdn(Mat B)
233 {
234 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
235 Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
236 PetscErrorCode ierr;
237
238 PetscFunctionBegin;
239 ierr = MatSetUp_LMVM(B);CHKERRQ(ierr);
240 if (!lbrdn->allocated) {
241 ierr = PetscMalloc2(lmvm->m, &lbrdn->sts, lmvm->m, &lbrdn->stq);CHKERRQ(ierr);
242 if (lmvm->m > 0) {
243 ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbrdn->P);CHKERRQ(ierr);
244 ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbrdn->Q);CHKERRQ(ierr);
245 }
246 lbrdn->allocated = PETSC_TRUE;
247 }
248 PetscFunctionReturn(0);
249 }
250
251 /*------------------------------------------------------------*/
252
MatCreate_LMVMBrdn(Mat B)253 PetscErrorCode MatCreate_LMVMBrdn(Mat B)
254 {
255 Mat_LMVM *lmvm;
256 Mat_Brdn *lbrdn;
257 PetscErrorCode ierr;
258
259 PetscFunctionBegin;
260 ierr = MatCreate_LMVM(B);CHKERRQ(ierr);
261 ierr = PetscObjectChangeTypeName((PetscObject)B, MATLMVMBROYDEN);CHKERRQ(ierr);
262 B->ops->setup = MatSetUp_LMVMBrdn;
263 B->ops->destroy = MatDestroy_LMVMBrdn;
264 B->ops->solve = MatSolve_LMVMBrdn;
265
266 lmvm = (Mat_LMVM*)B->data;
267 lmvm->square = PETSC_TRUE;
268 lmvm->ops->allocate = MatAllocate_LMVMBrdn;
269 lmvm->ops->reset = MatReset_LMVMBrdn;
270 lmvm->ops->mult = MatMult_LMVMBrdn;
271 lmvm->ops->update = MatUpdate_LMVMBrdn;
272 lmvm->ops->copy = MatCopy_LMVMBrdn;
273
274 ierr = PetscNewLog(B, &lbrdn);CHKERRQ(ierr);
275 lmvm->ctx = (void*)lbrdn;
276 lbrdn->allocated = PETSC_FALSE;
277 lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
278 PetscFunctionReturn(0);
279 }
280
281 /*------------------------------------------------------------*/
282
283 /*@
284 MatCreateLMVMBroyden - Creates a limited-memory "good" Broyden-type approximation
285 matrix used for a Jacobian. L-Brdn is not guaranteed to be symmetric or
286 positive-definite.
287
288 The provided local and global sizes must match the solution and function vectors
289 used with MatLMVMUpdate() and MatSolve(). The resulting L-Brdn matrix will have
290 storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
291 parallel. To use the L-Brdn matrix with other vector types, the matrix must be
292 created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
293 This ensures that the internal storage and work vectors are duplicated from the
294 correct type of vector.
295
296 Collective
297
298 Input Parameters:
299 + comm - MPI communicator, set to PETSC_COMM_SELF
300 . n - number of local rows for storage vectors
301 - N - global size of the storage vectors
302
303 Output Parameter:
304 . B - the matrix
305
306 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
307 paradigm instead of this routine directly.
308
309 Options Database Keys:
310 . -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
311
312 Level: intermediate
313
314 .seealso: MatCreate(), MATLMVM, MATLMVMBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
315 MatCreateLMVMBFGS(), MatCreateLMVMBadBrdn(), MatCreateLMVMSymBrdn()
316 @*/
MatCreateLMVMBroyden(MPI_Comm comm,PetscInt n,PetscInt N,Mat * B)317 PetscErrorCode MatCreateLMVMBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
318 {
319 PetscErrorCode ierr;
320
321 PetscFunctionBegin;
322 ierr = MatCreate(comm, B);CHKERRQ(ierr);
323 ierr = MatSetSizes(*B, n, n, N, N);CHKERRQ(ierr);
324 ierr = MatSetType(*B, MATLMVMBROYDEN);CHKERRQ(ierr);
325 ierr = MatSetUp(*B);CHKERRQ(ierr);
326 PetscFunctionReturn(0);
327 }
328