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:   lmvm->k        = -1;
  8:   lmvm->prev_set = PETSC_FALSE;
  9:   lmvm->shift    = 0.0;
 10:   if (destructive && lmvm->allocated) {
 11:     MatLMVMClearJ0(B);
 12:     B->rmap->n = B->rmap->N = B->cmap->n = B->cmap->N = 0;
 13:     VecDestroyVecs(lmvm->m, &lmvm->S);
 14:     VecDestroyVecs(lmvm->m, &lmvm->Y);
 15:     VecDestroy(&lmvm->Xprev);
 16:     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:   return 0;
 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:   PetscInt  m, n, M, N;
 33:   VecType   type;

 35:   if (lmvm->allocated) {
 36:     VecCheckMatCompatible(B, X, 2, F, 3);
 37:     VecGetType(X, &type);
 38:     PetscObjectTypeCompare((PetscObject)lmvm->Xprev, type, &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:       MatLMVMReset(B, PETSC_TRUE);
 44:     }
 45:   } else {
 46:     allocate = PETSC_TRUE;
 47:   }
 48:   if (allocate) {
 49:     VecGetLocalSize(X, &n);
 50:     VecGetSize(X, &N);
 51:     VecGetLocalSize(F, &m);
 52:     VecGetSize(F, &M);
 53:     B->rmap->n = m;
 54:     B->cmap->n = n;
 55:     B->rmap->N = M > -1 ? M : B->rmap->N;
 56:     B->cmap->N = N > -1 ? N : B->cmap->N;
 57:     VecDuplicate(X, &lmvm->Xprev);
 58:     VecDuplicate(F, &lmvm->Fprev);
 59:     if (lmvm->m > 0) {
 60:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
 61:       VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
 62:     }
 63:     lmvm->m_old     = lmvm->m;
 64:     lmvm->allocated = PETSC_TRUE;
 65:     B->preallocated = PETSC_TRUE;
 66:     B->assembled    = PETSC_TRUE;
 67:   }
 68:   return 0;
 69: }

 71: PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
 72: {
 73:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
 74:   PetscInt  i;
 75:   Vec       Stmp, Ytmp;

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

 98: PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
 99: {
100:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

102:   if (!lmvm->m) return 0;
103:   if (lmvm->prev_set) {
104:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
105:     VecAXPBY(lmvm->Xprev, 1.0, -1.0, X);
106:     VecAXPBY(lmvm->Fprev, 1.0, -1.0, F);
107:     /* Update S and Y */
108:     MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
109:   }

111:   /* Save the solution and function to be used in the next update */
112:   VecCopy(X, lmvm->Xprev);
113:   VecCopy(F, lmvm->Fprev);
114:   lmvm->prev_set = PETSC_TRUE;
115:   return 0;
116: }

118: static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
119: {
120:   MatMult(B, X, Z);
121:   VecAXPY(Z, 1.0, Y);
122:   return 0;
123: }

125: static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
126: {
127:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

129:   VecCheckSameSize(X, 2, Y, 3);
130:   VecCheckMatCompatible(B, X, 2, Y, 3);
132:   (*lmvm->ops->mult)(B, X, Y);
133:   if (lmvm->shift != 0.0) VecAXPY(Y, lmvm->shift, X);
134:   return 0;
135: }

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

144:   if (str == DIFFERENT_NONZERO_PATTERN) {
145:     MatLMVMReset(M, PETSC_TRUE);
146:     MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev);
147:   } else {
148:     MatLMVMIsAllocated(M, &allocatedM);
150:     MatCheckSameSize(B, 1, M, 2);
151:   }

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

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

187:   MatGetType(B, &lmvmType);
188:   MatCreate(PetscObjectComm((PetscObject)B), mat);
189:   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:   KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it);

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

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

210:   lmvm->shift += PetscRealPart(a);
211:   return 0;
212: }

214: static PetscErrorCode MatCreateVecs_LMVM(Mat B, Vec *L, Vec *R)
215: {
216:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

219:   if (L) VecDuplicate(lmvm->Xprev, L);
220:   if (R) VecDuplicate(lmvm->Fprev, R);
221:   return 0;
222: }

224: PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
225: {
226:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
227:   PetscBool isascii;
228:   MatType   type;

230:   PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii);
231:   if (isascii) {
232:     MatGetType(B, &type);
233:     PetscViewerASCIIPrintf(pv, "Max. storage: %" PetscInt_FMT "\n", lmvm->m);
234:     PetscViewerASCIIPrintf(pv, "Used storage: %" PetscInt_FMT "\n", lmvm->k + 1);
235:     PetscViewerASCIIPrintf(pv, "Number of updates: %" PetscInt_FMT "\n", lmvm->nupdates);
236:     PetscViewerASCIIPrintf(pv, "Number of rejects: %" PetscInt_FMT "\n", lmvm->nrejects);
237:     PetscViewerASCIIPrintf(pv, "Number of resets: %" PetscInt_FMT "\n", lmvm->nresets);
238:     if (lmvm->J0) {
239:       PetscViewerASCIIPrintf(pv, "J0 Matrix:\n");
240:       PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO);
241:       MatView(lmvm->J0, pv);
242:       PetscViewerPopFormat(pv);
243:     }
244:   }
245:   return 0;
246: }

248: PetscErrorCode MatSetFromOptions_LMVM(Mat B, PetscOptionItems *PetscOptionsObject)
249: {
250:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

252:   PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory Variable Metric matrix for approximating Jacobians");
253:   PetscOptionsInt("-mat_lmvm_hist_size", "number of past updates kept in memory for the approximation", "", lmvm->m, &lmvm->m, NULL);
254:   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);
255:   PetscOptionsReal("-mat_lmvm_eps", "(developer) machine zero definition", "", lmvm->eps, &lmvm->eps, NULL);
256:   PetscOptionsHeadEnd();
257:   KSPSetFromOptions(lmvm->J0ksp);
258:   return 0;
259: }

261: PetscErrorCode MatSetUp_LMVM(Mat B)
262: {
263:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
264:   PetscInt    m, n, M, N;
265:   PetscMPIInt size;
266:   MPI_Comm    comm = PetscObjectComm((PetscObject)B);

268:   MatGetSize(B, &M, &N);
270:   if (!lmvm->allocated) {
271:     MPI_Comm_size(comm, &size);
272:     if (size == 1) {
273:       VecCreateSeq(comm, N, &lmvm->Xprev);
274:       VecCreateSeq(comm, M, &lmvm->Fprev);
275:     } else {
276:       MatGetLocalSize(B, &m, &n);
277:       VecCreateMPI(comm, n, N, &lmvm->Xprev);
278:       VecCreateMPI(comm, m, M, &lmvm->Fprev);
279:     }
280:     if (lmvm->m > 0) {
281:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
282:       VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
283:     }
284:     lmvm->m_old     = lmvm->m;
285:     lmvm->allocated = PETSC_TRUE;
286:     B->preallocated = PETSC_TRUE;
287:     B->assembled    = PETSC_TRUE;
288:   }
289:   return 0;
290: }

292: PetscErrorCode MatDestroy_LMVM(Mat B)
293: {
294:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

296:   if (lmvm->allocated) {
297:     VecDestroyVecs(lmvm->m, &lmvm->S);
298:     VecDestroyVecs(lmvm->m, &lmvm->Y);
299:     VecDestroy(&lmvm->Xprev);
300:     VecDestroy(&lmvm->Fprev);
301:   }
302:   KSPDestroy(&lmvm->J0ksp);
303:   MatLMVMClearJ0(B);
304:   PetscFree(B->data);
305:   return 0;
306: }

308: PetscErrorCode MatCreate_LMVM(Mat B)
309: {
310:   Mat_LMVM *lmvm;

312:   PetscNew(&lmvm);
313:   B->data = (void *)lmvm;

315:   lmvm->m_old    = 0;
316:   lmvm->m        = 5;
317:   lmvm->k        = -1;
318:   lmvm->nupdates = 0;
319:   lmvm->nrejects = 0;
320:   lmvm->nresets  = 0;

322:   lmvm->ksp_max_it = 20;
323:   lmvm->ksp_rtol   = 0.0;
324:   lmvm->ksp_atol   = 0.0;

326:   lmvm->shift = 0.0;

328:   lmvm->eps        = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0);
329:   lmvm->allocated  = PETSC_FALSE;
330:   lmvm->prev_set   = PETSC_FALSE;
331:   lmvm->user_scale = PETSC_FALSE;
332:   lmvm->user_pc    = PETSC_FALSE;
333:   lmvm->user_ksp   = PETSC_FALSE;
334:   lmvm->square     = PETSC_FALSE;

336:   B->ops->destroy        = MatDestroy_LMVM;
337:   B->ops->setfromoptions = MatSetFromOptions_LMVM;
338:   B->ops->view           = MatView_LMVM;
339:   B->ops->setup          = MatSetUp_LMVM;
340:   B->ops->getvecs        = MatCreateVecs_LMVM;
341:   B->ops->shift          = MatShift_LMVM;
342:   B->ops->duplicate      = MatDuplicate_LMVM;
343:   B->ops->mult           = MatMult_LMVM;
344:   B->ops->multadd        = MatMultAdd_LMVM;
345:   B->ops->copy           = MatCopy_LMVM;

347:   lmvm->ops->update   = MatUpdate_LMVM;
348:   lmvm->ops->allocate = MatAllocate_LMVM;
349:   lmvm->ops->reset    = MatReset_LMVM;

351:   KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp);
352:   PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1);
353:   KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_");
354:   KSPSetType(lmvm->J0ksp, KSPGMRES);
355:   KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it);
356:   return 0;
357: }