Actual source code: lmvmimpl.c

  1: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>

  3: PetscErrorCode MatReset_LMVM(Mat B, PetscBool destructive)
  4: {
  5:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

  7:   PetscFunctionBegin;
  8:   lmvm->k        = -1;
  9:   lmvm->prev_set = PETSC_FALSE;
 10:   lmvm->shift    = 0.0;
 11:   if (destructive && lmvm->allocated) {
 12:     PetscCall(MatLMVMClearJ0(B));
 13:     PetscCall(VecDestroyVecs(lmvm->m, &lmvm->S));
 14:     PetscCall(VecDestroyVecs(lmvm->m, &lmvm->Y));
 15:     PetscCall(VecDestroy(&lmvm->Xprev));
 16:     PetscCall(VecDestroy(&lmvm->Fprev));
 17:     lmvm->nupdates  = 0;
 18:     lmvm->nrejects  = 0;
 19:     lmvm->m_old     = 0;
 20:     lmvm->allocated = PETSC_FALSE;
 21:     B->preallocated = PETSC_FALSE;
 22:     B->assembled    = PETSC_FALSE;
 23:   }
 24:   ++lmvm->nresets;
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: PetscErrorCode MatAllocate_LMVM(Mat B, Vec X, Vec F)
 29: {
 30:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
 31:   PetscBool same, allocate = PETSC_FALSE;
 32:   VecType   vtype;

 34:   PetscFunctionBegin;
 35:   if (lmvm->allocated) {
 36:     VecCheckMatCompatible(B, X, 2, F, 3);
 37:     PetscCall(VecGetType(X, &vtype));
 38:     PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->Xprev, vtype, &same));
 39:     if (!same) {
 40:       /* Given X vector has a different type than allocated X-type data structures.
 41:          We need to destroy all of this and duplicate again out of the given vector. */
 42:       allocate = PETSC_TRUE;
 43:       PetscCall(MatLMVMReset(B, PETSC_TRUE));
 44:     }
 45:   } else allocate = PETSC_TRUE;
 46:   if (allocate) {
 47:     PetscCall(VecGetType(X, &vtype));
 48:     PetscCall(MatSetVecType(B, vtype));
 49:     PetscCall(PetscLayoutReference(F->map, &B->rmap));
 50:     PetscCall(PetscLayoutReference(X->map, &B->cmap));
 51:     PetscCall(VecDuplicate(X, &lmvm->Xprev));
 52:     PetscCall(VecDuplicate(F, &lmvm->Fprev));
 53:     if (lmvm->m > 0) {
 54:       PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S));
 55:       PetscCall(VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y));
 56:     }
 57:     lmvm->m_old     = lmvm->m;
 58:     lmvm->allocated = PETSC_TRUE;
 59:     B->preallocated = PETSC_TRUE;
 60:     B->assembled    = PETSC_TRUE;
 61:   }
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

 65: PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
 66: {
 67:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
 68:   PetscInt  i;
 69:   Vec       Stmp, Ytmp;

 71:   PetscFunctionBegin;
 72:   if (lmvm->k == lmvm->m - 1) {
 73:     /* We hit the memory limit, so shift all the vectors back one spot
 74:        and shift the oldest to the front to receive the latest update. */
 75:     Stmp = lmvm->S[0];
 76:     Ytmp = lmvm->Y[0];
 77:     for (i = 0; i < lmvm->k; ++i) {
 78:       lmvm->S[i] = lmvm->S[i + 1];
 79:       lmvm->Y[i] = lmvm->Y[i + 1];
 80:     }
 81:     lmvm->S[lmvm->k] = Stmp;
 82:     lmvm->Y[lmvm->k] = Ytmp;
 83:   } else {
 84:     ++lmvm->k;
 85:   }
 86:   /* Put the precomputed update into the last vector */
 87:   PetscCall(VecCopy(S, lmvm->S[lmvm->k]));
 88:   PetscCall(VecCopy(Y, lmvm->Y[lmvm->k]));
 89:   ++lmvm->nupdates;
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

 93: PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
 94: {
 95:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

 97:   PetscFunctionBegin;
 98:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
 99:   if (lmvm->prev_set) {
100:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
101:     PetscCall(VecAXPBY(lmvm->Xprev, 1.0, -1.0, X));
102:     PetscCall(VecAXPBY(lmvm->Fprev, 1.0, -1.0, F));
103:     /* Update S and Y */
104:     PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
105:   }

107:   /* Save the solution and function to be used in the next update */
108:   PetscCall(VecCopy(X, lmvm->Xprev));
109:   PetscCall(VecCopy(F, lmvm->Fprev));
110:   lmvm->prev_set = PETSC_TRUE;
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
115: {
116:   PetscFunctionBegin;
117:   PetscCall(MatMult(B, X, Z));
118:   PetscCall(VecAXPY(Z, 1.0, Y));
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
123: {
124:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

126:   PetscFunctionBegin;
127:   VecCheckSameSize(X, 2, Y, 3);
128:   VecCheckMatCompatible(B, X, 2, Y, 3);
129:   PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
130:   PetscCall((*lmvm->ops->mult)(B, X, Y));
131:   PetscCall(VecAXPY(Y, lmvm->shift, X));
132:   PetscFunctionReturn(PETSC_SUCCESS);
133: }

135: static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
136: {
137:   Mat_LMVM *bctx = (Mat_LMVM *)B->data;
138:   Mat_LMVM *mctx;
139:   PetscInt  i;
140:   PetscBool allocatedM;

142:   PetscFunctionBegin;
143:   if (str == DIFFERENT_NONZERO_PATTERN) {
144:     PetscCall(MatLMVMReset(M, PETSC_TRUE));
145:     PetscCall(MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev));
146:   } else {
147:     PetscCall(MatLMVMIsAllocated(M, &allocatedM));
148:     PetscCheck(allocatedM, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Target matrix must be allocated first");
149:     MatCheckSameSize(B, 1, M, 2);
150:   }

152:   mctx = (Mat_LMVM *)M->data;
153:   if (bctx->user_pc) {
154:     PetscCall(MatLMVMSetJ0PC(M, bctx->J0pc));
155:   } else if (bctx->user_ksp) {
156:     PetscCall(MatLMVMSetJ0KSP(M, bctx->J0ksp));
157:   } else if (bctx->J0) {
158:     PetscCall(MatLMVMSetJ0(M, bctx->J0));
159:   } else if (bctx->user_scale) {
160:     if (bctx->J0diag) {
161:       PetscCall(MatLMVMSetJ0Diag(M, bctx->J0diag));
162:     } else {
163:       PetscCall(MatLMVMSetJ0Scale(M, bctx->J0scalar));
164:     }
165:   }
166:   mctx->nupdates = bctx->nupdates;
167:   mctx->nrejects = bctx->nrejects;
168:   mctx->k        = bctx->k;
169:   for (i = 0; i <= bctx->k; ++i) {
170:     PetscCall(VecCopy(bctx->S[i], mctx->S[i]));
171:     PetscCall(VecCopy(bctx->Y[i], mctx->Y[i]));
172:     PetscCall(VecCopy(bctx->Xprev, mctx->Xprev));
173:     PetscCall(VecCopy(bctx->Fprev, mctx->Fprev));
174:   }
175:   if (bctx->ops->copy) PetscCall((*bctx->ops->copy)(B, M, str));
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

179: static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
180: {
181:   Mat_LMVM *bctx = (Mat_LMVM *)B->data;
182:   Mat_LMVM *mctx;
183:   MatType   lmvmType;
184:   Mat       A;

186:   PetscFunctionBegin;
187:   PetscCall(MatGetType(B, &lmvmType));
188:   PetscCall(MatCreate(PetscObjectComm((PetscObject)B), mat));
189:   PetscCall(MatSetType(*mat, lmvmType));

191:   A                = *mat;
192:   mctx             = (Mat_LMVM *)A->data;
193:   mctx->m          = bctx->m;
194:   mctx->ksp_max_it = bctx->ksp_max_it;
195:   mctx->ksp_rtol   = bctx->ksp_rtol;
196:   mctx->ksp_atol   = bctx->ksp_atol;
197:   mctx->shift      = bctx->shift;
198:   PetscCall(KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it));

200:   PetscCall(MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev));
201:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(B, *mat, SAME_NONZERO_PATTERN));
202:   PetscFunctionReturn(PETSC_SUCCESS);
203: }

205: static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
206: {
207:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

209:   PetscFunctionBegin;
210:   PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
211:   lmvm->shift += PetscRealPart(a);
212:   PetscFunctionReturn(PETSC_SUCCESS);
213: }

215: PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
216: {
217:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
218:   PetscBool isascii;
219:   MatType   type;

221:   PetscFunctionBegin;
222:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
223:   if (isascii) {
224:     PetscCall(MatGetType(B, &type));
225:     PetscCall(PetscViewerASCIIPrintf(pv, "Max. storage: %" PetscInt_FMT "\n", lmvm->m));
226:     PetscCall(PetscViewerASCIIPrintf(pv, "Used storage: %" PetscInt_FMT "\n", lmvm->k + 1));
227:     PetscCall(PetscViewerASCIIPrintf(pv, "Number of updates: %" PetscInt_FMT "\n", lmvm->nupdates));
228:     PetscCall(PetscViewerASCIIPrintf(pv, "Number of rejects: %" PetscInt_FMT "\n", lmvm->nrejects));
229:     PetscCall(PetscViewerASCIIPrintf(pv, "Number of resets: %" PetscInt_FMT "\n", lmvm->nresets));
230:     if (lmvm->J0) {
231:       PetscCall(PetscViewerASCIIPrintf(pv, "J0 Matrix:\n"));
232:       PetscCall(PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO));
233:       PetscCall(MatView(lmvm->J0, pv));
234:       PetscCall(PetscViewerPopFormat(pv));
235:     }
236:   }
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: PetscErrorCode MatSetFromOptions_LMVM(Mat B, PetscOptionItems *PetscOptionsObject)
241: {
242:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

244:   PetscFunctionBegin;
245:   PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory Variable Metric matrix for approximating Jacobians");
246:   PetscCall(PetscOptionsInt("-mat_lmvm_hist_size", "number of past updates kept in memory for the approximation", "", lmvm->m, &lmvm->m, NULL));
247:   PetscCall(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));
248:   PetscCall(PetscOptionsReal("-mat_lmvm_eps", "(developer) machine zero definition", "", lmvm->eps, &lmvm->eps, NULL));
249:   PetscOptionsHeadEnd();
250:   PetscCall(KSPSetFromOptions(lmvm->J0ksp));
251:   PetscFunctionReturn(PETSC_SUCCESS);
252: }

254: PetscErrorCode MatSetUp_LMVM(Mat B)
255: {
256:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

258:   PetscFunctionBegin;
259:   if (!lmvm->allocated) {
260:     PetscCall(PetscLayoutSetUp(B->rmap));
261:     PetscCall(PetscLayoutSetUp(B->cmap));
262:     PetscCall(MatCreateVecs(B, &lmvm->Xprev, &lmvm->Fprev));
263:     if (lmvm->m > 0) {
264:       PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S));
265:       PetscCall(VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y));
266:     }
267:     lmvm->m_old     = lmvm->m;
268:     lmvm->allocated = PETSC_TRUE;
269:     B->preallocated = PETSC_TRUE;
270:     B->assembled    = PETSC_TRUE;
271:   }
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: PetscErrorCode MatDestroy_LMVM(Mat B)
276: {
277:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

279:   PetscFunctionBegin;
280:   if (lmvm->allocated) {
281:     PetscCall(VecDestroyVecs(lmvm->m, &lmvm->S));
282:     PetscCall(VecDestroyVecs(lmvm->m, &lmvm->Y));
283:     PetscCall(VecDestroy(&lmvm->Xprev));
284:     PetscCall(VecDestroy(&lmvm->Fprev));
285:   }
286:   PetscCall(KSPDestroy(&lmvm->J0ksp));
287:   PetscCall(MatLMVMClearJ0(B));
288:   PetscCall(PetscFree(B->data));
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: PetscErrorCode MatCreate_LMVM(Mat B)
293: {
294:   Mat_LMVM *lmvm;

296:   PetscFunctionBegin;
297:   PetscCall(PetscNew(&lmvm));
298:   B->data = (void *)lmvm;

300:   lmvm->m_old    = 0;
301:   lmvm->m        = 5;
302:   lmvm->k        = -1;
303:   lmvm->nupdates = 0;
304:   lmvm->nrejects = 0;
305:   lmvm->nresets  = 0;

307:   lmvm->ksp_max_it = 20;
308:   lmvm->ksp_rtol   = 0.0;
309:   lmvm->ksp_atol   = 0.0;

311:   lmvm->shift = 0.0;

313:   lmvm->eps        = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0);
314:   lmvm->allocated  = PETSC_FALSE;
315:   lmvm->prev_set   = PETSC_FALSE;
316:   lmvm->user_scale = PETSC_FALSE;
317:   lmvm->user_pc    = PETSC_FALSE;
318:   lmvm->user_ksp   = PETSC_FALSE;
319:   lmvm->square     = PETSC_FALSE;

321:   B->ops->destroy        = MatDestroy_LMVM;
322:   B->ops->setfromoptions = MatSetFromOptions_LMVM;
323:   B->ops->view           = MatView_LMVM;
324:   B->ops->setup          = MatSetUp_LMVM;
325:   B->ops->shift          = MatShift_LMVM;
326:   B->ops->duplicate      = MatDuplicate_LMVM;
327:   B->ops->mult           = MatMult_LMVM;
328:   B->ops->multadd        = MatMultAdd_LMVM;
329:   B->ops->copy           = MatCopy_LMVM;

331:   lmvm->ops->update   = MatUpdate_LMVM;
332:   lmvm->ops->allocate = MatAllocate_LMVM;
333:   lmvm->ops->reset    = MatReset_LMVM;

335:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp));
336:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1));
337:   PetscCall(KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_"));
338:   PetscCall(KSPSetType(lmvm->J0ksp, KSPGMRES));
339:   PetscCall(KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it));
340:   PetscFunctionReturn(PETSC_SUCCESS);
341: }