1 #include <../src/ksp/ksp/utils/lmvm/lmvm.h> /*I "petscksp.h" I*/
2
3 /*------------------------------------------------------------*/
4
MatReset_LMVM(Mat B,PetscBool destructive)5 PetscErrorCode MatReset_LMVM(Mat B, PetscBool destructive)
6 {
7 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
8 PetscErrorCode ierr;
9
10 PetscFunctionBegin;
11 lmvm->k = -1;
12 lmvm->prev_set = PETSC_FALSE;
13 lmvm->shift = 0.0;
14 if (destructive && lmvm->allocated) {
15 ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
16 B->rmap->n = B->rmap->N = B->cmap->n = B->cmap->N = 0;
17 ierr = VecDestroyVecs(lmvm->m, &lmvm->S);CHKERRQ(ierr);
18 ierr = VecDestroyVecs(lmvm->m, &lmvm->Y);CHKERRQ(ierr);
19 ierr = VecDestroy(&lmvm->Xprev);CHKERRQ(ierr);
20 ierr = VecDestroy(&lmvm->Fprev);CHKERRQ(ierr);
21 lmvm->nupdates = 0;
22 lmvm->nrejects = 0;
23 lmvm->m_old = 0;
24 lmvm->allocated = PETSC_FALSE;
25 B->preallocated = PETSC_FALSE;
26 B->assembled = PETSC_FALSE;
27 }
28 ++lmvm->nresets;
29 PetscFunctionReturn(0);
30 }
31
32 /*------------------------------------------------------------*/
33
MatAllocate_LMVM(Mat B,Vec X,Vec F)34 PetscErrorCode MatAllocate_LMVM(Mat B, Vec X, Vec F)
35 {
36 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
37 PetscErrorCode ierr;
38 PetscBool same, allocate = PETSC_FALSE;
39 PetscInt m, n, M, N;
40 VecType type;
41
42 PetscFunctionBegin;
43 if (lmvm->allocated) {
44 VecCheckMatCompatible(B, X, 2, F, 3);
45 ierr = VecGetType(X, &type);CHKERRQ(ierr);
46 ierr = PetscObjectTypeCompare((PetscObject)lmvm->Xprev, type, &same);CHKERRQ(ierr);
47 if (!same) {
48 /* Given X vector has a different type than allocated X-type data structures.
49 We need to destroy all of this and duplicate again out of the given vector. */
50 allocate = PETSC_TRUE;
51 ierr = MatLMVMReset(B, PETSC_TRUE);CHKERRQ(ierr);
52 }
53 } else {
54 allocate = PETSC_TRUE;
55 }
56 if (allocate) {
57 ierr = VecGetLocalSize(X, &n);CHKERRQ(ierr);
58 ierr = VecGetSize(X, &N);CHKERRQ(ierr);
59 ierr = VecGetLocalSize(F, &m);CHKERRQ(ierr);
60 ierr = VecGetSize(F, &M);CHKERRQ(ierr);
61 B->rmap->n = m;
62 B->cmap->n = n;
63 B->rmap->N = M > -1 ? M : B->rmap->N;
64 B->cmap->N = N > -1 ? N : B->cmap->N;
65 ierr = VecDuplicate(X, &lmvm->Xprev);CHKERRQ(ierr);
66 ierr = VecDuplicate(F, &lmvm->Fprev);CHKERRQ(ierr);
67 if (lmvm->m > 0) {
68 ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);CHKERRQ(ierr);
69 ierr = VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);CHKERRQ(ierr);
70 }
71 lmvm->m_old = lmvm->m;
72 lmvm->allocated = PETSC_TRUE;
73 B->preallocated = PETSC_TRUE;
74 B->assembled = PETSC_TRUE;
75 }
76 PetscFunctionReturn(0);
77 }
78
79 /*------------------------------------------------------------*/
80
MatUpdateKernel_LMVM(Mat B,Vec S,Vec Y)81 PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
82 {
83 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
84 PetscErrorCode ierr;
85 PetscInt i;
86 Vec Stmp, Ytmp;
87
88 PetscFunctionBegin;
89 if (lmvm->k == lmvm->m-1) {
90 /* We hit the memory limit, so shift all the vectors back one spot
91 and shift the oldest to the front to receive the latest update. */
92 Stmp = lmvm->S[0];
93 Ytmp = lmvm->Y[0];
94 for (i = 0; i < lmvm->k; ++i) {
95 lmvm->S[i] = lmvm->S[i+1];
96 lmvm->Y[i] = lmvm->Y[i+1];
97 }
98 lmvm->S[lmvm->k] = Stmp;
99 lmvm->Y[lmvm->k] = Ytmp;
100 } else {
101 ++lmvm->k;
102 }
103 /* Put the precomputed update into the last vector */
104 ierr = VecCopy(S, lmvm->S[lmvm->k]);CHKERRQ(ierr);
105 ierr = VecCopy(Y, lmvm->Y[lmvm->k]);CHKERRQ(ierr);
106 ++lmvm->nupdates;
107 PetscFunctionReturn(0);
108 }
109
110 /*------------------------------------------------------------*/
111
MatUpdate_LMVM(Mat B,Vec X,Vec F)112 PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
113 {
114 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
115 PetscErrorCode ierr;
116
117 PetscFunctionBegin;
118 if (!lmvm->m) PetscFunctionReturn(0);
119 if (lmvm->prev_set) {
120 /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
121 ierr = VecAXPBY(lmvm->Xprev, 1.0, -1.0, X);CHKERRQ(ierr);
122 ierr = VecAXPBY(lmvm->Fprev, 1.0, -1.0, F);CHKERRQ(ierr);
123 /* Update S and Y */
124 ierr = MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);CHKERRQ(ierr);
125 }
126
127 /* Save the solution and function to be used in the next update */
128 ierr = VecCopy(X, lmvm->Xprev);CHKERRQ(ierr);
129 ierr = VecCopy(F, lmvm->Fprev);CHKERRQ(ierr);
130 lmvm->prev_set = PETSC_TRUE;
131 PetscFunctionReturn(0);
132 }
133
134 /*------------------------------------------------------------*/
135
MatMultAdd_LMVM(Mat B,Vec X,Vec Y,Vec Z)136 static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
137 {
138 PetscErrorCode ierr;
139
140 PetscFunctionBegin;
141 ierr = MatMult(B, X, Z);CHKERRQ(ierr);
142 ierr = VecAXPY(Z, 1.0, Y);CHKERRQ(ierr);
143 PetscFunctionReturn(0);
144 }
145
146 /*------------------------------------------------------------*/
147
MatMult_LMVM(Mat B,Vec X,Vec Y)148 static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
149 {
150 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
151 PetscErrorCode ierr;
152
153 PetscFunctionBegin;
154 VecCheckSameSize(X, 2, Y, 3);
155 VecCheckMatCompatible(B, X, 2, Y, 3);
156 if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
157 ierr = (*lmvm->ops->mult)(B, X, Y);CHKERRQ(ierr);
158 if (lmvm->shift != 0.0) {
159 ierr = VecAXPY(Y, lmvm->shift, X);CHKERRQ(ierr);
160 }
161 PetscFunctionReturn(0);
162 }
163
164 /*------------------------------------------------------------*/
165
MatCopy_LMVM(Mat B,Mat M,MatStructure str)166 static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
167 {
168 Mat_LMVM *bctx = (Mat_LMVM*)B->data;
169 Mat_LMVM *mctx;
170 PetscErrorCode ierr;
171 PetscInt i;
172 PetscBool allocatedM;
173
174 PetscFunctionBegin;
175 if (str == DIFFERENT_NONZERO_PATTERN) {
176 ierr = MatLMVMReset(M, PETSC_TRUE);CHKERRQ(ierr);
177 ierr = MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev);CHKERRQ(ierr);
178 } else {
179 ierr = MatLMVMIsAllocated(M, &allocatedM);CHKERRQ(ierr);
180 if (!allocatedM) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Target matrix must be allocated first");
181 MatCheckSameSize(B, 1, M, 2);
182 }
183
184 mctx = (Mat_LMVM*)M->data;
185 if (bctx->user_pc) {
186 ierr = MatLMVMSetJ0PC(M, bctx->J0pc);CHKERRQ(ierr);
187 } else if (bctx->user_ksp) {
188 ierr = MatLMVMSetJ0KSP(M, bctx->J0ksp);CHKERRQ(ierr);
189 } else if (bctx->J0) {
190 ierr = MatLMVMSetJ0(M, bctx->J0);CHKERRQ(ierr);
191 } else if (bctx->user_scale) {
192 if (bctx->J0diag) {
193 ierr = MatLMVMSetJ0Diag(M, bctx->J0diag);CHKERRQ(ierr);
194 } else {
195 ierr = MatLMVMSetJ0Scale(M, bctx->J0scalar);CHKERRQ(ierr);
196 }
197 }
198 mctx->nupdates = bctx->nupdates;
199 mctx->nrejects = bctx->nrejects;
200 mctx->k = bctx->k;
201 for (i=0; i<=bctx->k; ++i) {
202 ierr = VecCopy(bctx->S[i], mctx->S[i]);CHKERRQ(ierr);
203 ierr = VecCopy(bctx->Y[i], mctx->Y[i]);CHKERRQ(ierr);
204 ierr = VecCopy(bctx->Xprev, mctx->Xprev);CHKERRQ(ierr);
205 ierr = VecCopy(bctx->Fprev, mctx->Fprev);CHKERRQ(ierr);
206 }
207 if (bctx->ops->copy) {
208 ierr = (*bctx->ops->copy)(B, M, str);CHKERRQ(ierr);
209 }
210 PetscFunctionReturn(0);
211 }
212
213 /*------------------------------------------------------------*/
214
MatDuplicate_LMVM(Mat B,MatDuplicateOption op,Mat * mat)215 static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
216 {
217 Mat_LMVM *bctx = (Mat_LMVM*)B->data;
218 Mat_LMVM *mctx;
219 PetscErrorCode ierr;
220 MatType lmvmType;
221 Mat A;
222
223 PetscFunctionBegin;
224 ierr = MatGetType(B, &lmvmType);CHKERRQ(ierr);
225 ierr = MatCreate(PetscObjectComm((PetscObject)B), mat);CHKERRQ(ierr);
226 ierr = MatSetType(*mat, lmvmType);CHKERRQ(ierr);
227
228 A = *mat;
229 mctx = (Mat_LMVM*)A->data;
230 mctx->m = bctx->m;
231 mctx->ksp_max_it = bctx->ksp_max_it;
232 mctx->ksp_rtol = bctx->ksp_rtol;
233 mctx->ksp_atol = bctx->ksp_atol;
234 mctx->shift = bctx->shift;
235 ierr = KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it);CHKERRQ(ierr);
236
237 ierr = MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev);CHKERRQ(ierr);
238 if (op == MAT_COPY_VALUES) {
239 ierr = MatCopy(B, *mat, SAME_NONZERO_PATTERN);CHKERRQ(ierr);
240 }
241 PetscFunctionReturn(0);
242 }
243
244 /*------------------------------------------------------------*/
245
MatShift_LMVM(Mat B,PetscScalar a)246 static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
247 {
248 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
249
250 PetscFunctionBegin;
251 if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
252 lmvm->shift += PetscRealPart(a);
253 PetscFunctionReturn(0);
254 }
255
256 /*------------------------------------------------------------*/
257
MatGetVecs_LMVM(Mat B,Vec * L,Vec * R)258 static PetscErrorCode MatGetVecs_LMVM(Mat B, Vec *L, Vec *R)
259 {
260 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
261 PetscErrorCode ierr;
262
263 PetscFunctionBegin;
264 if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
265 ierr = VecDuplicate(lmvm->Xprev, L);CHKERRQ(ierr);
266 ierr = VecDuplicate(lmvm->Fprev, R);CHKERRQ(ierr);
267 PetscFunctionReturn(0);
268 }
269
270 /*------------------------------------------------------------*/
271
MatView_LMVM(Mat B,PetscViewer pv)272 PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
273 {
274 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
275 PetscErrorCode ierr;
276 PetscBool isascii;
277 MatType type;
278
279 PetscFunctionBegin;
280 ierr = PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
281 if (isascii) {
282 ierr = MatGetType(B, &type);CHKERRQ(ierr);
283 ierr = PetscViewerASCIIPrintf(pv,"Max. storage: %D\n",lmvm->m);CHKERRQ(ierr);
284 ierr = PetscViewerASCIIPrintf(pv,"Used storage: %D\n",lmvm->k+1);CHKERRQ(ierr);
285 ierr = PetscViewerASCIIPrintf(pv,"Number of updates: %D\n",lmvm->nupdates);CHKERRQ(ierr);
286 ierr = PetscViewerASCIIPrintf(pv,"Number of rejects: %D\n",lmvm->nrejects);CHKERRQ(ierr);
287 ierr = PetscViewerASCIIPrintf(pv,"Number of resets: %D\n",lmvm->nresets);CHKERRQ(ierr);
288 if (lmvm->J0) {
289 ierr = PetscViewerASCIIPrintf(pv,"J0 Matrix:\n");CHKERRQ(ierr);
290 ierr = PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
291 ierr = MatView(lmvm->J0, pv);CHKERRQ(ierr);
292 ierr = PetscViewerPopFormat(pv);CHKERRQ(ierr);
293 }
294 }
295 PetscFunctionReturn(0);
296 }
297
298 /*------------------------------------------------------------*/
299
MatSetFromOptions_LMVM(PetscOptionItems * PetscOptionsObject,Mat B)300 PetscErrorCode MatSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject, Mat B)
301 {
302 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
303 PetscErrorCode ierr;
304
305 PetscFunctionBegin;
306 ierr = PetscOptionsHead(PetscOptionsObject,"Limited-memory Variable Metric matrix for approximating Jacobians");CHKERRQ(ierr);
307 ierr = PetscOptionsInt("-mat_lmvm_hist_size","number of past updates kept in memory for the approximation","",lmvm->m,&lmvm->m,NULL);CHKERRQ(ierr);
308 ierr = PetscOptionsInt("-mat_lmvm_ksp_its","(developer) fixed number of KSP iterations to take when inverting J0","",lmvm->ksp_max_it,&lmvm->ksp_max_it,NULL);CHKERRQ(ierr);
309 ierr = PetscOptionsReal("-mat_lmvm_eps","(developer) machine zero definition","",lmvm->eps,&lmvm->eps,NULL);CHKERRQ(ierr);
310 ierr = PetscOptionsTail();CHKERRQ(ierr);
311 ierr = KSPSetFromOptions(lmvm->J0ksp);CHKERRQ(ierr);
312 PetscFunctionReturn(0);
313 }
314
315 /*------------------------------------------------------------*/
316
MatSetUp_LMVM(Mat B)317 PetscErrorCode MatSetUp_LMVM(Mat B)
318 {
319 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
320 PetscErrorCode ierr;
321 PetscInt m, n, M, N;
322 PetscMPIInt size;
323 MPI_Comm comm = PetscObjectComm((PetscObject)B);
324
325 PetscFunctionBegin;
326 ierr = MatGetSize(B, &M, &N);CHKERRQ(ierr);
327 if (M == 0 && N == 0) SETERRQ(comm, PETSC_ERR_ORDER, "MatSetSizes() must be called before MatSetUp()");
328 if (!lmvm->allocated) {
329 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
330 if (size == 1) {
331 ierr = VecCreateSeq(comm, N, &lmvm->Xprev);CHKERRQ(ierr);
332 ierr = VecCreateSeq(comm, M, &lmvm->Fprev);CHKERRQ(ierr);
333 } else {
334 ierr = MatGetLocalSize(B, &m, &n);CHKERRQ(ierr);
335 ierr = VecCreateMPI(comm, n, N, &lmvm->Xprev);CHKERRQ(ierr);
336 ierr = VecCreateMPI(comm, m, M, &lmvm->Fprev);CHKERRQ(ierr);
337 }
338 if (lmvm->m > 0) {
339 ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);CHKERRQ(ierr);
340 ierr = VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);CHKERRQ(ierr);
341 }
342 lmvm->m_old = lmvm->m;
343 lmvm->allocated = PETSC_TRUE;
344 B->preallocated = PETSC_TRUE;
345 B->assembled = PETSC_TRUE;
346 }
347 PetscFunctionReturn(0);
348 }
349
350 /*------------------------------------------------------------*/
351
MatDestroy_LMVM(Mat B)352 PetscErrorCode MatDestroy_LMVM(Mat B)
353 {
354 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
355 PetscErrorCode ierr;
356
357 PetscFunctionBegin;
358 if (lmvm->allocated) {
359 ierr = VecDestroyVecs(lmvm->m, &lmvm->S);CHKERRQ(ierr);
360 ierr = VecDestroyVecs(lmvm->m, &lmvm->Y);CHKERRQ(ierr);
361 ierr = VecDestroy(&lmvm->Xprev);CHKERRQ(ierr);
362 ierr = VecDestroy(&lmvm->Fprev);CHKERRQ(ierr);
363 }
364 ierr = KSPDestroy(&lmvm->J0ksp);CHKERRQ(ierr);
365 ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
366 ierr = PetscFree(B->data);CHKERRQ(ierr);
367 PetscFunctionReturn(0);
368 }
369
370 /*------------------------------------------------------------*/
371
MatCreate_LMVM(Mat B)372 PetscErrorCode MatCreate_LMVM(Mat B)
373 {
374 Mat_LMVM *lmvm;
375 PetscErrorCode ierr;
376
377 PetscFunctionBegin;
378 ierr = PetscNewLog(B, &lmvm);CHKERRQ(ierr);
379 B->data = (void*)lmvm;
380
381 lmvm->m_old = 0;
382 lmvm->m = 5;
383 lmvm->k = -1;
384 lmvm->nupdates = 0;
385 lmvm->nrejects = 0;
386 lmvm->nresets = 0;
387
388 lmvm->ksp_max_it = 20;
389 lmvm->ksp_rtol = 0.0;
390 lmvm->ksp_atol = 0.0;
391
392 lmvm->shift = 0.0;
393
394 lmvm->eps = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0);
395 lmvm->allocated = PETSC_FALSE;
396 lmvm->prev_set = PETSC_FALSE;
397 lmvm->user_scale = PETSC_FALSE;
398 lmvm->user_pc = PETSC_FALSE;
399 lmvm->user_ksp = PETSC_FALSE;
400 lmvm->square = PETSC_FALSE;
401
402 B->ops->destroy = MatDestroy_LMVM;
403 B->ops->setfromoptions = MatSetFromOptions_LMVM;
404 B->ops->view = MatView_LMVM;
405 B->ops->setup = MatSetUp_LMVM;
406 B->ops->getvecs = MatGetVecs_LMVM;
407 B->ops->shift = MatShift_LMVM;
408 B->ops->duplicate = MatDuplicate_LMVM;
409 B->ops->mult = MatMult_LMVM;
410 B->ops->multadd = MatMultAdd_LMVM;
411 B->ops->copy = MatCopy_LMVM;
412
413 lmvm->ops->update = MatUpdate_LMVM;
414 lmvm->ops->allocate = MatAllocate_LMVM;
415 lmvm->ops->reset = MatReset_LMVM;
416
417 ierr = KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp);CHKERRQ(ierr);
418 ierr = PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1);CHKERRQ(ierr);
419 ierr = KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_");CHKERRQ(ierr);
420 ierr = KSPSetType(lmvm->J0ksp, KSPGMRES);CHKERRQ(ierr);
421 ierr = KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it);CHKERRQ(ierr);
422 PetscFunctionReturn(0);
423 }
424