Actual source code: pod.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
  1: #include <petsc/private/kspimpl.h>
  2: #include <petsc/private/matimpl.h>
  3: #include <petscblaslapack.h>
  4: static PetscBool  cited = PETSC_FALSE;
  5: static const char citation[] =
  6: "@phdthesis{zampini2010non,\n"
  7: "  title={Non-overlapping Domain Decomposition Methods for Cardiac Reaction-Diffusion Models and Applications},\n"
  8: "  author={Zampini, S},\n"
  9: "  year={2010},\n"
 10: "  school={PhD thesis, Universita degli Studi di Milano}\n"
 11: "}\n";

 13: typedef struct {
 14:   PetscInt     maxn;             /* maximum number of snapshots */
 15:   PetscInt     n;                /* number of active snapshots */
 16:   PetscInt     curr;             /* current tip of snapshots set */
 17:   Vec          *xsnap;           /* snapshots */
 18:   Vec          *bsnap;           /* rhs snapshots */
 19:   Vec          *work;            /* parallel work vectors */
 20:   PetscScalar  *dots_iallreduce;
 21:   MPI_Request  req_iallreduce;
 22:   PetscInt     ndots_iallreduce; /* if we have iallreduce we can hide the VecMDot communications */
 23:   PetscReal    tol;              /* relative tolerance to retain eigenvalues */
 24:   PetscBool    Aspd;             /* if true, uses the SPD operator as inner product */
 25:   PetscScalar  *corr;            /* correlation matrix */
 26:   PetscReal    *eigs;            /* eigenvalues */
 27:   PetscScalar  *eigv;            /* eigenvectors */
 28:   PetscBLASInt nen;              /* dimension of lower dimensional system */
 29:   PetscInt     st;               /* first eigenvector of correlation matrix to be retained */
 30:   PetscBLASInt *iwork;           /* integer work vector */
 31:   PetscScalar  *yhay;            /* Y^H * A * Y */
 32:   PetscScalar  *low;             /* lower dimensional linear system */
 33: #if defined(PETSC_USE_COMPLEX)
 34:   PetscReal    *rwork;
 35: #endif
 36:   PetscBLASInt lwork;
 37:   PetscScalar  *swork;
 38:   PetscBool    monitor;
 39: } KSPGuessPOD;

 41: static PetscErrorCode KSPGuessReset_POD(KSPGuess guess)
 42: {
 43:   KSPGuessPOD    *pod = (KSPGuessPOD*)guess->data;
 45:   PetscLayout    Alay = NULL,vlay = NULL;
 46:   PetscBool      cong;

 49:   pod->nen  = 0;
 50:   pod->n    = 0;
 51:   pod->curr = 0;
 52:   /* need to wait for completion of outstanding requests */
 53:   if (pod->ndots_iallreduce) {
 54:     MPI_Wait(&pod->req_iallreduce,MPI_STATUS_IGNORE);
 55:   }
 56:   pod->ndots_iallreduce = 0;
 57:   /* destroy vectors if the size of the linear system has changed */
 58:   if (guess->A) {
 59:     MatGetLayouts(guess->A,&Alay,NULL);
 60:   }
 61:   if (pod->xsnap) {
 62:     VecGetLayout(pod->xsnap[0],&vlay);
 63:   }
 64:   cong = PETSC_FALSE;
 65:   if (vlay && Alay) {
 66:     PetscLayoutCompare(Alay,vlay,&cong);
 67:   }
 68:   if (!cong) {
 69:     VecDestroyVecs(pod->maxn,&pod->xsnap);
 70:     VecDestroyVecs(pod->maxn,&pod->bsnap);
 71:     VecDestroyVecs(1,&pod->work);
 72:   }
 73:   return(0);
 74: }

 76: static PetscErrorCode KSPGuessSetUp_POD(KSPGuess guess)
 77: {
 78:   KSPGuessPOD    *pod = (KSPGuessPOD*)guess->data;

 82:   if (!pod->corr) {
 83:     PetscScalar  sdummy;
 84:     PetscReal    rdummy = 0;
 85:     PetscBLASInt bN,lierr,idummy;

 87:     PetscCalloc6(pod->maxn*pod->maxn,&pod->corr,pod->maxn,&pod->eigs,pod->maxn*pod->maxn,&pod->eigv,
 88:                         6*pod->maxn,&pod->iwork,pod->maxn*pod->maxn,&pod->yhay,pod->maxn*pod->maxn,&pod->low);
 89: #if defined(PETSC_USE_COMPLEX)
 90:     PetscMalloc1(7*pod->maxn,&pod->rwork);
 91: #endif
 92: #if defined(PETSC_HAVE_MPI_IALLREDUCE)
 93:     PetscMalloc1(3*pod->maxn,&pod->dots_iallreduce);
 94: #endif
 95:     pod->lwork = -1;
 96:     PetscBLASIntCast(pod->maxn,&bN);
 97: #if !defined(PETSC_USE_COMPLEX)
 98:     PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","L",&bN,pod->corr,&bN,&rdummy,&rdummy,&idummy,&idummy,
 99:                                                   &rdummy,&idummy,pod->eigs,pod->eigv,&bN,&sdummy,&pod->lwork,pod->iwork,pod->iwork+5*bN,&lierr));
100: #else
101:     PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","L",&bN,pod->corr,&bN,&rdummy,&rdummy,&idummy,&idummy,
102:                                                   &rdummy,&idummy,pod->eigs,pod->eigv,&bN,&sdummy,&pod->lwork,pod->rwork,pod->iwork,pod->iwork+5*bN,&lierr));
103: #endif
104:     if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SYEV Lapack routine %d",(int)lierr);
105:     PetscBLASIntCast((PetscInt)PetscRealPart(sdummy),&pod->lwork);
106:     PetscMalloc1(pod->lwork+PetscMax(bN*bN,6*bN),&pod->swork);
107:   }
108:   /* work vectors are sequential, we explicitly use MPI_Allreduce */
109:   if (!pod->xsnap) {
110:     VecType   type;
111:     Vec       *v,vseq;
112:     PetscInt  n;

114:     KSPCreateVecs(guess->ksp,1,&v,0,NULL);
115:     VecCreate(PETSC_COMM_SELF,&vseq);
116:     VecGetLocalSize(v[0],&n);
117:     VecSetSizes(vseq,n,n);
118:     VecGetType(v[0],&type);
119:     VecSetType(vseq,type);
120:     VecDestroyVecs(1,&v);
121:     VecDuplicateVecs(vseq,pod->maxn,&pod->xsnap);
122:     VecDestroy(&vseq);
123:     PetscLogObjectParents(guess,pod->maxn,pod->xsnap);
124:   }
125:   if (!pod->bsnap) {
126:     VecDuplicateVecs(pod->xsnap[0],pod->maxn,&pod->bsnap);
127:     PetscLogObjectParents(guess,pod->maxn,pod->bsnap);
128:   }
129:   if (!pod->work) {
130:     KSPCreateVecs(guess->ksp,1,&pod->work,0,NULL);
131:   }
132:   return(0);
133: }

135: static PetscErrorCode KSPGuessDestroy_POD(KSPGuess guess)
136: {
137:   KSPGuessPOD *pod = (KSPGuessPOD*)guess->data;
138:   PetscErrorCode  ierr;

141:   PetscFree6(pod->corr,pod->eigs,pod->eigv,pod->iwork,
142:                     pod->yhay,pod->low);
143: #if defined(PETSC_USE_COMPLEX)
144:   PetscFree(pod->rwork);
145: #endif
146:   /* need to wait for completion before destroying dots_iallreduce */
147:   if (pod->ndots_iallreduce) {
148:     MPI_Wait(&pod->req_iallreduce,MPI_STATUS_IGNORE);
149:   }
150:   PetscFree(pod->dots_iallreduce);
151:   PetscFree(pod->swork);
152:   VecDestroyVecs(pod->maxn,&pod->bsnap);
153:   VecDestroyVecs(pod->maxn,&pod->xsnap);
154:   VecDestroyVecs(1,&pod->work);
155:   PetscFree(pod);
156:   return(0);
157: }

159: static PetscErrorCode KSPGuessUpdate_POD(KSPGuess,Vec,Vec);

161: static PetscErrorCode KSPGuessFormGuess_POD(KSPGuess guess,Vec b,Vec x)
162: {
163:   KSPGuessPOD    *pod = (KSPGuessPOD*)guess->data;
165:   PetscScalar    one = 1, zero = 0, *array;
166:   PetscBLASInt   bN,ione = 1,bNen,lierr;
167:   PetscInt       i;

170:   PetscCitationsRegister(citation,&cited);
171:   if (pod->ndots_iallreduce) { /* complete communication and project the linear system */
172:     KSPGuessUpdate_POD(guess,NULL,NULL);
173:   }
174:   if (!pod->nen) return(0);
175:   /* b_low = S * V^T * X^T * b */
176:   VecGetArrayRead(b,(const PetscScalar**)&array);
177:   VecPlaceArray(pod->bsnap[pod->curr],array);
178:   VecRestoreArrayRead(b,(const PetscScalar**)&array);
179:   VecMDot(pod->bsnap[pod->curr],pod->n,pod->xsnap,pod->swork);
180:   VecResetArray(pod->bsnap[pod->curr]);
181:   MPIU_Allreduce(pod->swork,pod->swork + pod->n,pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
182:   PetscBLASIntCast(pod->n,&bN);
183:   PetscBLASIntCast(pod->nen,&bNen);
184:   PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&bN,&bNen,&one,pod->eigv+pod->st*pod->n,&bN,pod->swork+pod->n,&ione,&zero,pod->swork,&ione));
185:   if (pod->monitor) {
186:     PetscPrintf(PetscObjectComm((PetscObject)guess),"  KSPGuessPOD alphas = ");
187:     for (i=0; i<pod->nen; i++) {
188: #if defined(PETSC_USE_COMPLEX)
189:       PetscPrintf(PetscObjectComm((PetscObject)guess),"%g + %g i",(double)PetscRealPart(pod->swork[i]),(double)PetscImaginaryPart(pod->swork[i]));
190: #else
191:       PetscPrintf(PetscObjectComm((PetscObject)guess),"%g ",(double)pod->swork[i]);
192: #endif
193:     }
194:     PetscPrintf(PetscObjectComm((PetscObject)guess),"\n");
195:   }
196:   /* A_low x_low = b_low */
197:   if (!pod->Aspd) { /* A is spd -> LOW = Identity */
198:     KSP       pksp = guess->ksp;
199:     PetscBool tsolve,symm;

201:     if (pod->monitor) {
202:       PetscMPIInt rank;
203:       Mat         L;

205:       MPI_Comm_rank(PetscObjectComm((PetscObject)guess),&rank);
206:       MatCreateSeqDense(PETSC_COMM_SELF,pod->nen,pod->nen,pod->low,&L);
207:       if (!rank) {
208:         MatView(L,NULL);
209:       }
210:       MatDestroy(&L);
211:     }
212:     MatGetOption(guess->A,MAT_SYMMETRIC,&symm);
213:     tsolve = symm ? PETSC_FALSE : pksp->transpose_solve;
214:     PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&bNen,&bNen,pod->low,&bNen,pod->iwork,&lierr));
215:     if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)lierr);
216:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_(tsolve ? "T" : "N",&bNen,&ione,pod->low,&bNen,pod->iwork,pod->swork,&bNen,&lierr));
217:     if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRS Lapack routine %d",(int)lierr);
218:   }
219:   /* x = X * V * S * x_low */
220:   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&bN,&bNen,&one,pod->eigv+pod->st*pod->n,&bN,pod->swork,&ione,&zero,pod->swork+pod->n,&ione));
221:   if (pod->monitor) {
222:     PetscPrintf(PetscObjectComm((PetscObject)guess),"  KSPGuessPOD sol = ");
223:     for (i=0; i<pod->nen; i++) {
224: #if defined(PETSC_USE_COMPLEX)
225:       PetscPrintf(PetscObjectComm((PetscObject)guess),"%g + %g i",(double)PetscRealPart(pod->swork[i+pod->n]),(double)PetscImaginaryPart(pod->swork[i+pod->n]));
226: #else
227:       PetscPrintf(PetscObjectComm((PetscObject)guess),"%g ",(double)pod->swork[i+pod->n]);
228: #endif
229:     }
230:     PetscPrintf(PetscObjectComm((PetscObject)guess),"\n");
231:   }
232:   VecGetArray(x,&array);
233:   VecPlaceArray(pod->bsnap[pod->curr],array);
234:   VecSet(pod->bsnap[pod->curr],0);
235:   VecMAXPY(pod->bsnap[pod->curr],pod->n,pod->swork+pod->n,pod->xsnap);
236:   VecResetArray(pod->bsnap[pod->curr]);
237:   VecRestoreArray(x,&array);
238:   return(0);
239: }

241: static PetscErrorCode KSPGuessUpdate_POD(KSPGuess guess, Vec b, Vec x)
242: {
243:   KSPGuessPOD    *pod = (KSPGuessPOD*)guess->data;
244:   PetscScalar    one = 1, zero = 0,*array;
245:   PetscReal      toten, parten, reps = 0; /* dlamch? */
246:   PetscBLASInt   bN,lierr,idummy;
247:   PetscInt       i;

251:   if (pod->ndots_iallreduce) goto complete_request;
252:   pod->n = pod->n < pod->maxn ? pod->n+1 : pod->maxn;
253:   VecCopy(x,pod->xsnap[pod->curr]);
254:   VecGetArray(pod->bsnap[pod->curr],&array);
255:   VecPlaceArray(pod->work[0],array);
256:   KSP_MatMult(guess->ksp,guess->A,x,pod->work[0]);
257:   VecResetArray(pod->work[0]);
258:   VecRestoreArray(pod->bsnap[pod->curr],&array);
259:   if (pod->Aspd) {
260:     VecMDot(pod->xsnap[pod->curr],pod->n,pod->bsnap,pod->swork);
261: #if !defined(PETSC_HAVE_MPI_IALLREDUCE)
262:     MPIU_Allreduce(pod->swork,pod->swork + 3*pod->n,pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
263: #else
264:     MPI_Iallreduce(pod->swork,pod->dots_iallreduce,pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess),&pod->req_iallreduce);
265:     pod->ndots_iallreduce = 1;
266: #endif
267:   } else {
268:     PetscInt  off;
269:     PetscBool herm;

271: #if defined(PETSC_USE_COMPLEX)
272:     MatGetOption(guess->A,MAT_HERMITIAN,&herm);
273: #else
274:     MatGetOption(guess->A,MAT_SYMMETRIC,&herm);
275: #endif
276:     off = (guess->ksp->transpose_solve && !herm) ? 2*pod->n : pod->n;

278:     /* TODO: we may want to use a user-defined dot for the correlation matrix */
279:     VecMDot(pod->xsnap[pod->curr],pod->n,pod->xsnap,pod->swork);
280:     VecMDot(pod->bsnap[pod->curr],pod->n,pod->xsnap,pod->swork + off);
281:     if (!herm) {
282:       off  = (off == pod->n) ? 2*pod->n : pod->n;
283:       VecMDot(pod->xsnap[pod->curr],pod->n,pod->bsnap,pod->swork + off);
284: #if !defined(PETSC_HAVE_MPI_IALLREDUCE)
285:       MPIU_Allreduce(pod->swork,pod->swork + 3*pod->n,3*pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
286: #else
287:       MPI_Iallreduce(pod->swork,pod->dots_iallreduce,3*pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess),&pod->req_iallreduce);
288:       pod->ndots_iallreduce = 3;
289: #endif
290:     } else {
291: #if !defined(PETSC_HAVE_MPI_IALLREDUCE)
292:       MPIU_Allreduce(pod->swork,pod->swork + 3*pod->n,2*pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
293:       for (i=0;i<pod->n;i++) pod->swork[5*pod->n + i] = pod->swork[4*pod->n + i];
294: #else
295:       MPI_Iallreduce(pod->swork,pod->dots_iallreduce,2*pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess),&pod->req_iallreduce);
296:       pod->ndots_iallreduce = 2;
297: #endif
298:     }
299:   }
300:   if (pod->ndots_iallreduce) return(0);

302: complete_request:
303:   if (pod->ndots_iallreduce) {
304:     MPI_Wait(&pod->req_iallreduce,MPI_STATUS_IGNORE);
305:     switch (pod->ndots_iallreduce) {
306:     case 3:
307:       for (i=0;i<pod->n;i++) pod->swork[3*pod->n + i] = pod->dots_iallreduce[         i];
308:       for (i=0;i<pod->n;i++) pod->swork[4*pod->n + i] = pod->dots_iallreduce[  pod->n+i];
309:       for (i=0;i<pod->n;i++) pod->swork[5*pod->n + i] = pod->dots_iallreduce[2*pod->n+i];
310:       break;
311:     case 2:
312:       for (i=0;i<pod->n;i++) pod->swork[3*pod->n + i] = pod->dots_iallreduce[       i];
313:       for (i=0;i<pod->n;i++) pod->swork[4*pod->n + i] = pod->dots_iallreduce[pod->n+i];
314:       for (i=0;i<pod->n;i++) pod->swork[5*pod->n + i] = pod->dots_iallreduce[pod->n+i];
315:       break;
316:     case 1:
317:       for (i=0;i<pod->n;i++) pod->swork[3*pod->n + i] = pod->dots_iallreduce[i];
318:       break;
319:     default:
320:       SETERRQ1(PetscObjectComm((PetscObject)guess),PETSC_ERR_PLIB,"Invalid number of outstanding dots operations: %D",pod->ndots_iallreduce);
321:     }
322:   }
323:   pod->ndots_iallreduce = 0;

325:   /* correlation matrix and Y^H A Y (Galerkin) */
326:   for (i=0;i<pod->n;i++) {
327:     pod->corr[pod->curr*pod->maxn+i] = pod->swork[3*pod->n + i];
328:     pod->corr[i*pod->maxn+pod->curr] = PetscConj(pod->swork[3*pod->n + i]);
329:     if (!pod->Aspd) {
330:       pod->yhay[pod->curr*pod->maxn+i] = pod->swork[4*pod->n + i];
331:       pod->yhay[i*pod->maxn+pod->curr] = PetscConj(pod->swork[5*pod->n + i]);
332:     }
333:   }
334:   /* syevx change the input matrix */
335:   for (i=0;i<pod->n;i++) {
336:     PetscInt j;
337:     for (j=i;j<pod->n;j++) pod->swork[i*pod->n+j] = pod->corr[i*pod->maxn+j];
338:   }
339:   PetscBLASIntCast(pod->n,&bN);
340: #if !defined(PETSC_USE_COMPLEX)
341:   PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","L",&bN,pod->swork,&bN,
342:                                                 &reps,&reps,&idummy,&idummy,
343:                                                 &reps,&idummy,pod->eigs,pod->eigv,&bN,
344:                                                 pod->swork+bN*bN,&pod->lwork,pod->iwork,pod->iwork+5*bN,&lierr));
345: #else
346:   PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","L",&bN,pod->swork,&bN,
347:                                                 &reps,&reps,&idummy,&idummy,
348:                                                 &reps,&idummy,pod->eigs,pod->eigv,&bN,
349:                                                 pod->swork+bN*bN,&pod->lwork,pod->rwork,pod->iwork,pod->iwork+5*bN,&lierr));
350: #endif
351:   if (lierr<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYEV Lapack routine: illegal argument %d",-(int)lierr);
352:   else if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYEV Lapack routine: %d eigenvectors failed to converge",(int)lierr);

354:   /* dimension of lower dimensional system */
355:   pod->st = -1;
356:   for (i=0,toten=0;i<pod->n;i++) {
357:     pod->eigs[i] = PetscMax(pod->eigs[i],0.0);
358:     toten += pod->eigs[i];
359:     if (!pod->eigs[i]) pod->st = i;
360:   }
361:   pod->nen = 0;
362:   for (i=pod->n-1,parten=0;i>pod->st && toten > 0;i--) {
363:     pod->nen++;
364:     parten += pod->eigs[i];
365:     if (parten + toten*pod->tol >= toten) break;
366:   }
367:   pod->st = pod->n - pod->nen;

369:   /* Compute eigv = V * S */
370:   for (i=pod->st;i<pod->n;i++) {
371:     const PetscReal v = 1.0/PetscSqrtReal(pod->eigs[i]);
372:     const PetscInt  st = pod->n*i;
373:     PetscInt        j;

375:     for (j=0;j<pod->n;j++) pod->eigv[st+j] *= v;
376:   }

378:   /* compute S * V^T * X^T * A * X * V * S if needed */
379:   if (pod->nen && !pod->Aspd) {
380:     PetscBLASInt bNen,bMaxN;
381:     PetscInt     st = pod->st*pod->n;
382:     PetscBLASIntCast(pod->nen,&bNen);
383:     PetscBLASIntCast(pod->maxn,&bMaxN);
384:     PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&bNen,&bN,&bN,&one,pod->eigv+st,&bN,pod->yhay,&bMaxN,&zero,pod->swork,&bNen));
385:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&bNen,&bNen,&bN,&one,pod->swork,&bNen,pod->eigv+st,&bN,&zero,pod->low,&bNen));
386:   }

388:   if (pod->monitor) {
389:     PetscPrintf(PetscObjectComm((PetscObject)guess),"  KSPGuessPOD: basis %D, energy fractions = ",pod->nen);
390:     for (i=pod->n-1;i>=0;i--) {
391:       PetscPrintf(PetscObjectComm((PetscObject)guess),"%1.6e (%d) ",pod->eigs[i]/toten,i >= pod->st ? 1 : 0);
392:     }
393:     PetscPrintf(PetscObjectComm((PetscObject)guess),"\n");
394:     if (PetscDefined(USE_DEBUG)) {
395:       for (i=0;i<pod->n;i++) {
396:         Vec v;
397:         PetscInt j;
398:         PetscBLASInt bNen,ione = 1;

400:         VecDuplicate(pod->xsnap[i],&v);
401:         VecCopy(pod->xsnap[i],v);
402:         PetscBLASIntCast(pod->nen,&bNen);
403:         PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&bN,&bNen,&one,pod->eigv+pod->st*pod->n,&bN,pod->corr+pod->maxn*i,&ione,&zero,pod->swork,&ione));
404:         PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&bN,&bNen,&one,pod->eigv+pod->st*pod->n,&bN,pod->swork,&ione,&zero,pod->swork+pod->n,&ione));
405:         for (j=0;j<pod->n;j++) pod->swork[j] = -pod->swork[pod->n+j];
406:         VecMAXPY(v,pod->n,pod->swork,pod->xsnap);
407:         VecDot(v,v,pod->swork);
408:         MPIU_Allreduce(pod->swork,pod->swork + 1,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
409:         PetscPrintf(PetscObjectComm((PetscObject)guess),"  Error projection %D: %g (expected lower than %g)\n",i,(double)PetscRealPart(pod->swork[1]),(double)(toten-parten));
410:         VecDestroy(&v);
411:       }
412:     }
413:   }
414:   /* new tip */
415:   pod->curr = (pod->curr+1)%pod->maxn;
416:   return(0);
417: }

419: static PetscErrorCode KSPGuessSetFromOptions_POD(KSPGuess guess)
420: {
421:   KSPGuessPOD    *pod = (KSPGuessPOD *)guess->data;

425:   PetscOptionsBegin(PetscObjectComm((PetscObject)guess),((PetscObject)guess)->prefix,"POD initial guess options","KSPGuess");
426:   PetscOptionsInt("-ksp_guess_pod_size","Number of snapshots",NULL,pod->maxn,&pod->maxn,NULL);
427:   PetscOptionsBool("-ksp_guess_pod_monitor","Monitor initial guess generator",NULL,pod->monitor,&pod->monitor,NULL);
428:   PetscOptionsReal("-ksp_guess_pod_tol","Tolerance to retain eigenvectors",NULL,pod->tol,&pod->tol,NULL);
429:   PetscOptionsBool("-ksp_guess_pod_Ainner","Use the operator as inner product (must be SPD)",NULL,pod->Aspd,&pod->Aspd,NULL);
430:   PetscOptionsEnd();
431:   return(0);
432: }

434: static PetscErrorCode KSPGuessView_POD(KSPGuess guess,PetscViewer viewer)
435: {
436:   KSPGuessPOD    *pod = (KSPGuessPOD*)guess->data;
437:   PetscBool      isascii;

441:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
442:   if (isascii) {
443:     PetscViewerASCIIPrintf(viewer,"Max size %D, tolerance %g, Ainner %d\n",pod->maxn,pod->tol,pod->Aspd);
444:   }
445:   return(0);
446: }

448: /*
449:     KSPGUESSPOD - Implements a proper orthogonal decomposition based Galerkin scheme for repeated linear system solves.

451:   The initial guess is obtained by solving a small and dense linear system, obtained by Galerkin projection on a lower dimensional space generated by the previous solutions.
452:   The number of solutions to be retained and the energy tolerance to construct the lower dimensional basis can be specified at command line by -ksp_guess_pod_tol <real> and -ksp_guess_pod_size <int>.

454:   References:
455: .   1. - http://www.math.uni-konstanz.de/numerik/personen/volkwein/teaching/POD-Book.pdf

457:     Level: intermediate

459: .seealso: KSPGuess, KSPGuessType, KSPGuessCreate(), KSPSetGuess(), KSPGetGuess()
460: @*/
461: PetscErrorCode KSPGuessCreate_POD(KSPGuess guess)
462: {
463:   KSPGuessPOD    *pod;

467:   PetscNewLog(guess,&pod);
468:   pod->maxn   = 10;
469:   pod->tol    = PETSC_MACHINE_EPSILON;
470:   guess->data = pod;

472:   guess->ops->setfromoptions = KSPGuessSetFromOptions_POD;
473:   guess->ops->destroy        = KSPGuessDestroy_POD;
474:   guess->ops->setup          = KSPGuessSetUp_POD;
475:   guess->ops->view           = KSPGuessView_POD;
476:   guess->ops->reset          = KSPGuessReset_POD;
477:   guess->ops->update         = KSPGuessUpdate_POD;
478:   guess->ops->formguess      = KSPGuessFormGuess_POD;
479:   return(0);
480: }