Actual source code: lmvmimpl.c

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

  3: /*------------------------------------------------------------*/

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

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

 30: /*------------------------------------------------------------*/

 32: PetscErrorCode MatAllocate_LMVM(Mat B, Vec X, Vec F)
 33: {
 34:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
 35:   PetscBool         same, allocate = PETSC_FALSE;
 36:   PetscInt          m, n, M, N;
 37:   VecType           type;

 39:   if (lmvm->allocated) {
 40:     VecCheckMatCompatible(B, X, 2, F, 3);
 41:     VecGetType(X, &type);
 42:     PetscObjectTypeCompare((PetscObject)lmvm->Xprev, type, &same);
 43:     if (!same) {
 44:       /* Given X vector has a different type than allocated X-type data structures.
 45:          We need to destroy all of this and duplicate again out of the given vector. */
 46:       allocate = PETSC_TRUE;
 47:       MatLMVMReset(B, PETSC_TRUE);
 48:     }
 49:   } else {
 50:     allocate = PETSC_TRUE;
 51:   }
 52:   if (allocate) {
 53:     VecGetLocalSize(X, &n);
 54:     VecGetSize(X, &N);
 55:     VecGetLocalSize(F, &m);
 56:     VecGetSize(F, &M);
 57:     B->rmap->n = m;
 58:     B->cmap->n = n;
 59:     B->rmap->N = M > -1 ? M : B->rmap->N;
 60:     B->cmap->N = N > -1 ? N : B->cmap->N;
 61:     VecDuplicate(X, &lmvm->Xprev);
 62:     VecDuplicate(F, &lmvm->Fprev);
 63:     if (lmvm->m > 0) {
 64:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
 65:       VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
 66:     }
 67:     lmvm->m_old = lmvm->m;
 68:     lmvm->allocated = PETSC_TRUE;
 69:     B->preallocated = PETSC_TRUE;
 70:     B->assembled = PETSC_TRUE;
 71:   }
 72:   return 0;
 73: }

 75: /*------------------------------------------------------------*/

 77: PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
 78: {
 79:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
 80:   PetscInt          i;
 81:   Vec               Stmp, Ytmp;

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

104: /*------------------------------------------------------------*/

106: PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
107: {
108:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;

110:   if (!lmvm->m) return 0;
111:   if (lmvm->prev_set) {
112:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
113:     VecAXPBY(lmvm->Xprev, 1.0, -1.0, X);
114:     VecAXPBY(lmvm->Fprev, 1.0, -1.0, F);
115:     /* Update S and Y */
116:     MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
117:   }

119:   /* Save the solution and function to be used in the next update */
120:   VecCopy(X, lmvm->Xprev);
121:   VecCopy(F, lmvm->Fprev);
122:   lmvm->prev_set = PETSC_TRUE;
123:   return 0;
124: }

126: /*------------------------------------------------------------*/

128: static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
129: {
130:   MatMult(B, X, Z);
131:   VecAXPY(Z, 1.0, Y);
132:   return 0;
133: }

135: /*------------------------------------------------------------*/

137: static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
138: {
139:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;

141:   VecCheckSameSize(X, 2, Y, 3);
142:   VecCheckMatCompatible(B, X, 2, Y, 3);
144:   (*lmvm->ops->mult)(B, X, Y);
145:   if (lmvm->shift != 0.0) {
146:     VecAXPY(Y, lmvm->shift, X);
147:   }
148:   return 0;
149: }

151: /*------------------------------------------------------------*/

153: static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
154: {
155:   Mat_LMVM          *bctx = (Mat_LMVM*)B->data;
156:   Mat_LMVM          *mctx;
157:   PetscInt          i;
158:   PetscBool         allocatedM;

160:   if (str == DIFFERENT_NONZERO_PATTERN) {
161:     MatLMVMReset(M, PETSC_TRUE);
162:     MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev);
163:   } else {
164:     MatLMVMIsAllocated(M, &allocatedM);
166:     MatCheckSameSize(B, 1, M, 2);
167:   }

169:   mctx = (Mat_LMVM*)M->data;
170:   if (bctx->user_pc) {
171:     MatLMVMSetJ0PC(M, bctx->J0pc);
172:   } else if (bctx->user_ksp) {
173:     MatLMVMSetJ0KSP(M, bctx->J0ksp);
174:   } else if (bctx->J0) {
175:     MatLMVMSetJ0(M, bctx->J0);
176:   } else if (bctx->user_scale) {
177:     if (bctx->J0diag) {
178:       MatLMVMSetJ0Diag(M, bctx->J0diag);
179:     } else {
180:       MatLMVMSetJ0Scale(M, bctx->J0scalar);
181:     }
182:   }
183:   mctx->nupdates = bctx->nupdates;
184:   mctx->nrejects = bctx->nrejects;
185:   mctx->k = bctx->k;
186:   for (i=0; i<=bctx->k; ++i) {
187:     VecCopy(bctx->S[i], mctx->S[i]);
188:     VecCopy(bctx->Y[i], mctx->Y[i]);
189:     VecCopy(bctx->Xprev, mctx->Xprev);
190:     VecCopy(bctx->Fprev, mctx->Fprev);
191:   }
192:   if (bctx->ops->copy) {
193:     (*bctx->ops->copy)(B, M, str);
194:   }
195:   return 0;
196: }

198: /*------------------------------------------------------------*/

200: static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
201: {
202:   Mat_LMVM          *bctx = (Mat_LMVM*)B->data;
203:   Mat_LMVM          *mctx;
204:   MatType           lmvmType;
205:   Mat               A;

207:   MatGetType(B, &lmvmType);
208:   MatCreate(PetscObjectComm((PetscObject)B), mat);
209:   MatSetType(*mat, lmvmType);

211:   A = *mat;
212:   mctx = (Mat_LMVM*)A->data;
213:   mctx->m = bctx->m;
214:   mctx->ksp_max_it = bctx->ksp_max_it;
215:   mctx->ksp_rtol = bctx->ksp_rtol;
216:   mctx->ksp_atol = bctx->ksp_atol;
217:   mctx->shift = bctx->shift;
218:   KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it);

220:   MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev);
221:   if (op == MAT_COPY_VALUES) {
222:     MatCopy(B, *mat, SAME_NONZERO_PATTERN);
223:   }
224:   return 0;
225: }

227: /*------------------------------------------------------------*/

229: static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
230: {
231:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;

234:   lmvm->shift += PetscRealPart(a);
235:   return 0;
236: }

238: /*------------------------------------------------------------*/

240: static PetscErrorCode MatGetVecs_LMVM(Mat B, Vec *L, Vec *R)
241: {
242:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;

245:   VecDuplicate(lmvm->Xprev, L);
246:   VecDuplicate(lmvm->Fprev, R);
247:   return 0;
248: }

250: /*------------------------------------------------------------*/

252: PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
253: {
254:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
255:   PetscBool         isascii;
256:   MatType           type;

258:   PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
259:   if (isascii) {
260:     MatGetType(B, &type);
261:     PetscViewerASCIIPrintf(pv,"Max. storage: %D\n",lmvm->m);
262:     PetscViewerASCIIPrintf(pv,"Used storage: %D\n",lmvm->k+1);
263:     PetscViewerASCIIPrintf(pv,"Number of updates: %D\n",lmvm->nupdates);
264:     PetscViewerASCIIPrintf(pv,"Number of rejects: %D\n",lmvm->nrejects);
265:     PetscViewerASCIIPrintf(pv,"Number of resets: %D\n",lmvm->nresets);
266:     if (lmvm->J0) {
267:       PetscViewerASCIIPrintf(pv,"J0 Matrix:\n");
268:       PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO);
269:       MatView(lmvm->J0, pv);
270:       PetscViewerPopFormat(pv);
271:     }
272:   }
273:   return 0;
274: }

276: /*------------------------------------------------------------*/

278: PetscErrorCode MatSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject, Mat B)
279: {
280:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;

282:   PetscOptionsHead(PetscOptionsObject,"Limited-memory Variable Metric matrix for approximating Jacobians");
283:   PetscOptionsInt("-mat_lmvm_hist_size","number of past updates kept in memory for the approximation","",lmvm->m,&lmvm->m,NULL);
284:   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);
285:   PetscOptionsReal("-mat_lmvm_eps","(developer) machine zero definition","",lmvm->eps,&lmvm->eps,NULL);
286:   PetscOptionsTail();
287:   KSPSetFromOptions(lmvm->J0ksp);
288:   return 0;
289: }

291: /*------------------------------------------------------------*/

293: PetscErrorCode MatSetUp_LMVM(Mat B)
294: {
295:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
296:   PetscInt          m, n, M, N;
297:   PetscMPIInt       size;
298:   MPI_Comm          comm = PetscObjectComm((PetscObject)B);

300:   MatGetSize(B, &M, &N);
302:   if (!lmvm->allocated) {
303:     MPI_Comm_size(comm, &size);
304:     if (size == 1) {
305:       VecCreateSeq(comm, N, &lmvm->Xprev);
306:       VecCreateSeq(comm, M, &lmvm->Fprev);
307:     } else {
308:       MatGetLocalSize(B, &m, &n);
309:       VecCreateMPI(comm, n, N, &lmvm->Xprev);
310:       VecCreateMPI(comm, m, M, &lmvm->Fprev);
311:     }
312:     if (lmvm->m > 0) {
313:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
314:       VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
315:     }
316:     lmvm->m_old = lmvm->m;
317:     lmvm->allocated = PETSC_TRUE;
318:     B->preallocated = PETSC_TRUE;
319:     B->assembled = PETSC_TRUE;
320:   }
321:   return 0;
322: }

324: /*------------------------------------------------------------*/

326: PetscErrorCode MatDestroy_LMVM(Mat B)
327: {
328:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;

330:   if (lmvm->allocated) {
331:     VecDestroyVecs(lmvm->m, &lmvm->S);
332:     VecDestroyVecs(lmvm->m, &lmvm->Y);
333:     VecDestroy(&lmvm->Xprev);
334:     VecDestroy(&lmvm->Fprev);
335:   }
336:   KSPDestroy(&lmvm->J0ksp);
337:   MatLMVMClearJ0(B);
338:   PetscFree(B->data);
339:   return 0;
340: }

342: /*------------------------------------------------------------*/

344: PetscErrorCode MatCreate_LMVM(Mat B)
345: {
346:   Mat_LMVM          *lmvm;

348:   PetscNewLog(B, &lmvm);
349:   B->data = (void*)lmvm;

351:   lmvm->m_old = 0;
352:   lmvm->m = 5;
353:   lmvm->k = -1;
354:   lmvm->nupdates = 0;
355:   lmvm->nrejects = 0;
356:   lmvm->nresets = 0;

358:   lmvm->ksp_max_it = 20;
359:   lmvm->ksp_rtol = 0.0;
360:   lmvm->ksp_atol = 0.0;

362:   lmvm->shift = 0.0;

364:   lmvm->eps = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0);
365:   lmvm->allocated = PETSC_FALSE;
366:   lmvm->prev_set = PETSC_FALSE;
367:   lmvm->user_scale = PETSC_FALSE;
368:   lmvm->user_pc = PETSC_FALSE;
369:   lmvm->user_ksp = PETSC_FALSE;
370:   lmvm->square = PETSC_FALSE;

372:   B->ops->destroy = MatDestroy_LMVM;
373:   B->ops->setfromoptions = MatSetFromOptions_LMVM;
374:   B->ops->view = MatView_LMVM;
375:   B->ops->setup = MatSetUp_LMVM;
376:   B->ops->getvecs = MatGetVecs_LMVM;
377:   B->ops->shift = MatShift_LMVM;
378:   B->ops->duplicate = MatDuplicate_LMVM;
379:   B->ops->mult = MatMult_LMVM;
380:   B->ops->multadd = MatMultAdd_LMVM;
381:   B->ops->copy = MatCopy_LMVM;

383:   lmvm->ops->update = MatUpdate_LMVM;
384:   lmvm->ops->allocate = MatAllocate_LMVM;
385:   lmvm->ops->reset = MatReset_LMVM;

387:   KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp);
388:   PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1);
389:   KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_");
390:   KSPSetType(lmvm->J0ksp, KSPGMRES);
391:   KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it);
392:   return 0;
393: }