2: #define PETSCKSP_DLL 4: #include <../src/ksp/ksp/impls/gmres/agmres/agmresimpl.h> 6: #define AGMRES_DEFAULT_MAXK 30 7: #define AGMRES_DELTA_DIRECTIONS 10 8: static PetscErrorCode KSPAGMRESBuildSoln(KSP,PetscInt);
9: static PetscErrorCode KSPAGMRESBuildBasis(KSP);
10: static PetscErrorCode KSPAGMRESBuildHessenberg(KSP);
12: PetscLogEvent KSP_AGMRESComputeDeflationData, KSP_AGMRESBuildBasis, KSP_AGMRESComputeShifts, KSP_AGMRESRoddec;
14: extern PetscErrorCode KSPSetUp_DGMRES(KSP);
15: extern PetscErrorCode KSPBuildSolution_DGMRES(KSP,Vec,Vec*);
16: extern PetscErrorCode KSPSolve_DGMRES(KSP);
17: extern PetscErrorCode KSPDGMRESComputeDeflationData_DGMRES(KSP,PetscInt*);
18: extern PetscErrorCode KSPDGMRESComputeSchurForm_DGMRES(KSP,PetscInt*);
19: extern PetscErrorCode KSPDGMRESApplyDeflation_DGMRES(KSP,Vec,Vec);
20: extern PetscErrorCode KSPDestroy_DGMRES(KSP);
21: extern PetscErrorCode KSPSetFromOptions_DGMRES(PetscOptionItems *,KSP);
22: extern PetscErrorCode KSPDGMRESSetEigen_DGMRES(KSP,PetscInt);
23: /*
24: * This function allocates data for the Newton basis GMRES implementation.
25: * Note that most data are allocated in KSPSetUp_DGMRES and KSPSetUp_GMRES, including the space for the basis vectors, the various Hessenberg matrices and the Givens rotations coefficients
26: *
27: */
28: static PetscErrorCode KSPSetUp_AGMRES(KSP ksp) 29: {
30: PetscErrorCode ierr;
31: PetscInt hes;
32: PetscInt nloc;
33: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
34: PetscInt neig = agmres->neig;
35: const PetscInt max_k = agmres->max_k;
36: PetscInt N = MAXKSPSIZE;
37: PetscInt lwork = PetscMax(8 * N + 16, 4 * neig * (N - neig));
40: if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"no symmetric preconditioning for KSPAGMRES");
41: N = MAXKSPSIZE;
42: /* Preallocate space during the call to KSPSetup_GMRES for the Krylov basis */
43: agmres->q_preallocate = PETSC_TRUE; /* No allocation on the fly */
44: /* Preallocate space to compute later the eigenvalues in GMRES */
45: ksp->calc_sings = PETSC_TRUE;
46: agmres->max_k = N; /* Set the augmented size to be allocated in KSPSetup_GMRES */
47: KSPSetUp_DGMRES(ksp);
48: agmres->max_k = max_k;
49: hes = (N + 1) * (N + 1);
51: /* Data for the Newton basis GMRES */
52: PetscCalloc4(max_k,&agmres->Rshift,max_k,&agmres->Ishift,hes,&agmres->Rloc,((N+1)*4),&agmres->wbufptr);
53: PetscMalloc7((N+1),&agmres->Scale,(N+1),&agmres->sgn,(N+1),&agmres->tloc,(N+1),&agmres->temp,(N+1),&agmres->tau,lwork,&agmres->work,(N+1),&agmres->nrs);
54: PetscMemzero(agmres->Scale, (N+1)*sizeof(PetscScalar));
55: PetscMemzero(agmres->sgn, (N+1)*sizeof(PetscScalar));
56: PetscMemzero(agmres->tloc, (N+1)*sizeof(PetscScalar));
57: PetscMemzero(agmres->temp, (N+1)*sizeof(PetscScalar));
59: /* Allocate space for the vectors in the orthogonalized basis*/
60: VecGetLocalSize(agmres->vecs[0], &nloc);
61: PetscMalloc1(nloc*(N+1), &agmres->Qloc);
63: /* Init the ring of processors for the roddec orthogonalization */
64: KSPAGMRESRoddecInitNeighboor(ksp);
66: if (agmres->neig < 1) return(0);
68: /* Allocate space for the deflation */
69: PetscMalloc1(N, &agmres->select);
70: VecDuplicateVecs(VEC_V(0), N, &agmres->TmpU);
71: PetscMalloc2(N*N, &agmres->MatEigL, N*N, &agmres->MatEigR);
72: /* PetscMalloc6(N*N, &agmres->Q, N*N, &agmres->Z, N, &agmres->wr, N, &agmres->wi, N, &agmres->beta, N, &agmres->modul); */
73: PetscMalloc3(N*N, &agmres->Q, N*N, &agmres->Z, N, &agmres->beta);
74: PetscMalloc2((N+1),&agmres->perm,(2*neig*N),&agmres->iwork);
75: return(0);
76: }
78: /* This function returns the current solution from the private data structure of AGMRES back to ptr.
79: * This function is provided to be compliant with the KSP GMRES scheme.
80: *
81: */
82: static PetscErrorCode KSPBuildSolution_AGMRES(KSP ksp,Vec ptr, Vec *result) 83: {
84: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
88: if (!ptr) {
89: if (!agmres->sol_temp) {
90: VecDuplicate(ksp->vec_sol,&agmres->sol_temp);
91: VecCopy(ksp->vec_sol,agmres->sol_temp);
92: PetscLogObjectParent((PetscObject)ksp,(PetscObject)agmres->sol_temp);
93: }
94: ptr = agmres->sol_temp;
95: } else {
96: VecCopy(ksp->vec_sol, ptr);
97: }
98: if (result) *result = ptr;
99: return(0);
100: }
102: /* This function computes the shifts needed to generate stable basis vectors (through the Newton polynomials)
103: * At input, the operators (matrix and preconditioners) are used to create a new GMRES KSP.
104: * One cycle of GMRES with the Arnoldi process is performed and the eigenvalues of the induced Hessenberg matrix (the Ritz values) are computed.
105: * NOTE: This function is not currently used; the next function is rather used when the eigenvectors are needed next to augment the basis
106: */
107: PetscErrorCode KSPComputeShifts_GMRES(KSP ksp)108: {
109: PetscErrorCode ierr;
110: KSP_AGMRES *agmres = (KSP_AGMRES*)(ksp->data);
111: KSP kspgmres;
112: Mat Amat, Pmat;
113: const PetscInt max_k = agmres->max_k;
114: PC pc;
115: PetscInt m;
116: PetscScalar *Rshift, *Ishift;
117: PetscBool flg;
120: /* Perform one cycle of classical GMRES (with the Arnoldi process) to get the Hessenberg matrix
121: We assume here that the ksp is AGMRES and that the operators for the
122: linear system have been set in this ksp */
123: KSPCreate(PetscObjectComm((PetscObject)ksp), &kspgmres);
124: if (!ksp->pc) { KSPGetPC(ksp,&ksp->pc); }
125: PCGetOperators(ksp->pc, &Amat, &Pmat);
126: KSPSetOperators(kspgmres, Amat, Pmat);
127: KSPSetFromOptions(kspgmres);
128: PetscOptionsHasName(NULL,((PetscObject)ksp)->prefix, "-ksp_view", &flg);
129: if (flg) { PetscOptionsClearValue(NULL,"-ksp_view"); }
130: KSPSetType(kspgmres, KSPGMRES);
131: KSPGMRESSetRestart(kspgmres, max_k);
132: KSPGetPC(ksp, &pc);
133: KSPSetPC(kspgmres, pc);
134: /* Copy common options */
135: kspgmres->pc_side = ksp->pc_side;
136: /* Setup KSP context */
137: KSPSetComputeEigenvalues(kspgmres, PETSC_TRUE);
138: KSPSetUp(kspgmres);
140: kspgmres->max_it = max_k; /* Restrict the maximum number of iterations to one cycle of GMRES */
141: kspgmres->rtol = ksp->rtol;
143: KSPSolve(kspgmres, ksp->vec_rhs, ksp->vec_sol);
145: ksp->guess_zero = PETSC_FALSE;
146: ksp->rnorm = kspgmres->rnorm;
147: ksp->its = kspgmres->its;
148: if (kspgmres->reason == KSP_CONVERGED_RTOL) {
149: ksp->reason = KSP_CONVERGED_RTOL;
150: return(0);
151: } else ksp->reason = KSP_CONVERGED_ITERATING;
152: /* Now, compute the Shifts values */
153: PetscMalloc2(max_k,&Rshift,max_k,&Ishift);
154: KSPComputeEigenvalues(kspgmres, max_k, Rshift, Ishift, &m);
155: if (m < max_k) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB, "Unable to compute the Shifts for the Newton basis");
156: else {
157: KSPAGMRESLejaOrdering(Rshift, Ishift, agmres->Rshift, agmres->Ishift, max_k);
159: agmres->HasShifts = PETSC_TRUE;
160: }
161: /* Restore KSP view options */
162: if (flg) { PetscOptionsSetValue(NULL,"-ksp_view", ""); }
163: return(0);
164: }
166: /* This function computes the shift values (Ritz values) needed to generate stable basis vectors
167: * One cycle of DGMRES is performed to find the eigenvalues. The same data structures are used since AGMRES extends DGMRES
168: * Note that when the basis is to be augmented, then this function computes the harmonic Ritz vectors from this first cycle.
169: * Input :
170: * - The operators (matrix, preconditioners and right hand side) are normally required.
171: * - max_k : the size of the (non augmented) basis.
172: * - neig: The number of eigenvectors to augment, if deflation is needed
173: * Output :
174: * - The shifts as complex pair of arrays in wr and wi (size max_k).
175: * - The harmonic Ritz vectors (agmres->U) if deflation is needed.
176: */
177: static PetscErrorCode KSPComputeShifts_DGMRES(KSP ksp)178: {
180: KSP_AGMRES *agmres = (KSP_AGMRES*)(ksp->data);
181: PetscInt max_k = agmres->max_k; /* size of the (non augmented) Krylov subspace */
182: PetscInt Neig = 0;
183: const PetscInt max_it = ksp->max_it;
184: PetscBool flg;
186: /* Perform one cycle of dgmres to find the eigenvalues and compute the first approximations of the eigenvectors */
189: PetscLogEventBegin(KSP_AGMRESComputeShifts, ksp, 0,0,0);
190: /* Send the size of the augmented basis to DGMRES */
191: ksp->max_it = max_k; /* set this to have DGMRES performing only one cycle */
192: ksp->ops->buildsolution = KSPBuildSolution_DGMRES;
193: KSPSolve_DGMRES(ksp);
194: ksp->guess_zero = PETSC_FALSE;
195: if (ksp->reason == KSP_CONVERGED_RTOL) {
196: PetscLogEventEnd(KSP_AGMRESComputeShifts, ksp, 0,0,0);
197: return(0);
198: } else ksp->reason = KSP_CONVERGED_ITERATING;
200: if ((agmres->r == 0) && (agmres->neig > 0)) { /* Compute the eigenvalues for the shifts and the eigenvectors (to augment the Newton basis) */
201: agmres->HasSchur = PETSC_FALSE;
202: KSPDGMRESComputeDeflationData_DGMRES (ksp, &Neig);CHKERRQ (ierr);
203: Neig = max_k;
204: } else { /* From DGMRES, compute only the eigenvalues needed as Shifts for the Newton Basis */
205: KSPDGMRESComputeSchurForm_DGMRES(ksp, &Neig);
206: }
208: /* It may happen that the Ritz values from one cycle of GMRES are not accurate enough to provide a good stability. In this case, another cycle of GMRES is performed. The two sets of values thus generated are sorted and the most accurate are kept as shifts */
209: PetscOptionsHasName(NULL,NULL, "-ksp_agmres_ImproveShifts", &flg);
210: if (!flg) {
211: KSPAGMRESLejaOrdering(agmres->wr, agmres->wi, agmres->Rshift, agmres->Ishift, max_k);
212: } else { /* Perform another cycle of DGMRES to find another set of eigenvalues */
213: PetscInt i;
214: PetscScalar *wr, *wi,*Rshift, *Ishift;
215: PetscMalloc4(2*max_k, &wr, 2*max_k, &wi, 2*max_k, &Rshift, 2*max_k, &Ishift);
216: for (i = 0; i < max_k; i++) {
217: wr[i] = agmres->wr[i];
218: wi[i] = agmres->wi[i];
219: }
221: KSPSolve_DGMRES(ksp);
223: ksp->guess_zero = PETSC_FALSE;
224: if (ksp->reason == KSP_CONVERGED_RTOL) return(0);
225: else ksp->reason = KSP_CONVERGED_ITERATING;
226: if (agmres->neig > 0) { /* Compute the eigenvalues for the shifts) and the eigenvectors (to augment the Newton basis */
227: agmres->HasSchur = PETSC_FALSE;
229: KSPDGMRESComputeDeflationData_DGMRES(ksp, &Neig);
230: Neig = max_k;
231: } else { /* From DGMRES, compute only the eigenvalues needed as Shifts for the Newton Basis */
232: KSPDGMRESComputeSchurForm_DGMRES(ksp, &Neig);
233: }
234: for (i = 0; i < max_k; i++) {
235: wr[max_k+i] = agmres->wr[i];
236: wi[max_k+i] = agmres->wi[i];
237: }
238: KSPAGMRESLejaOrdering(wr, wi, Rshift, Ishift, 2*max_k);
239: for (i = 0; i< max_k; i++) {
240: agmres->Rshift[i] = Rshift[i];
241: agmres->Ishift[i] = Ishift[i];
242: }
243: PetscFree(Rshift);
244: PetscFree(wr);
245: PetscFree(Ishift);
246: PetscFree(wi);
247: }
249: agmres->HasShifts = PETSC_TRUE;
250: ksp->max_it = max_it;
251: PetscLogEventEnd(KSP_AGMRESComputeShifts, ksp, 0,0,0);
252: return(0);
253: }
255: /*
256: * Generate the basis vectors from the Newton polynomials with shifts and scaling factors
257: * The scaling factors are computed to obtain unit vectors. Note that this step can be avoided with the preprocessing option KSP_AGMRES_NONORM.
258: * Inputs :
259: * - Operators (Matrix and preconditioners and the first basis vector in VEC_V(0)
260: * - Shifts values in agmres->Rshift and agmres->Ishift.
261: * Output :
262: * - agmres->vecs or VEC_V : basis vectors
263: * - agmres->Scale : Scaling factors (equal to 1 if no scaling is done)
264: */
265: static PetscErrorCode KSPAGMRESBuildBasis(KSP ksp)266: {
268: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
269: PetscReal *Rshift = agmres->Rshift;
270: PetscReal *Ishift = agmres->Ishift;
271: PetscReal *Scale = agmres->Scale;
272: const PetscInt max_k = agmres->max_k;
273: PetscInt KspSize = KSPSIZE; /* if max_k == KspSizen then the basis should not be augmented */
274: PetscInt j = 1;
277: PetscLogEventBegin(KSP_AGMRESBuildBasis, ksp, 0,0,0);
278: Scale[0] = 1.0;
279: while (j <= max_k) {
280: if (Ishift[j-1] == 0) {
281: if ((ksp->pc_side == PC_LEFT) && agmres->r && agmres->DeflPrecond) {
282: /* Apply the precond-matrix operators */
283: KSP_PCApplyBAorAB(ksp, VEC_V(j-1), VEC_TMP, VEC_TMP_MATOP);
284: /* Then apply deflation as a preconditioner */
285: KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_TMP, VEC_V(j));
286: } else if ((ksp->pc_side == PC_RIGHT) && agmres->r && agmres->DeflPrecond) {
287: KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_V(j-1), VEC_TMP);
288: KSP_PCApplyBAorAB(ksp, VEC_TMP, VEC_V(j), VEC_TMP_MATOP);
289: } else {
290: KSP_PCApplyBAorAB(ksp, VEC_V(j-1), VEC_V(j), VEC_TMP_MATOP);
291: }
292: VecAXPY(VEC_V(j), -Rshift[j-1], VEC_V(j-1));
293: #if defined(KSP_AGMRES_NONORM)
294: Scale[j] = 1.0;
295: #else
296: VecScale(VEC_V(j), Scale[j-1]); /* This step can be postponed until all vectors are built */
297: VecNorm(VEC_V(j), NORM_2, &(Scale[j]));
298: Scale[j] = 1.0/Scale[j];
299: #endif
301: agmres->matvecs += 1;
302: j++;
303: } else {
304: if ((ksp->pc_side == PC_LEFT) && agmres->r && agmres->DeflPrecond) {
305: /* Apply the precond-matrix operators */
306: KSP_PCApplyBAorAB(ksp, VEC_V(j-1), VEC_TMP, VEC_TMP_MATOP);
307: /* Then apply deflation as a preconditioner */
308: KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_TMP, VEC_V(j));
309: } else if ((ksp->pc_side == PC_RIGHT) && agmres->r && agmres->DeflPrecond) {
310: KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_V(j-1), VEC_TMP);
311: KSP_PCApplyBAorAB(ksp, VEC_TMP, VEC_V(j), VEC_TMP_MATOP);
312: } else {
313: KSP_PCApplyBAorAB(ksp, VEC_V(j-1), VEC_V(j), VEC_TMP_MATOP);
314: }
315: VecAXPY(VEC_V(j), -Rshift[j-1], VEC_V(j-1));
316: #if defined(KSP_AGMRES_NONORM)
317: Scale[j] = 1.0;
318: #else
319: VecScale(VEC_V(j), Scale[j-1]);
320: VecNorm(VEC_V(j), NORM_2, &(Scale[j]));
321: Scale[j] = 1.0/Scale[j];
322: #endif
323: agmres->matvecs += 1;
324: j++;
325: if ((ksp->pc_side == PC_LEFT) && agmres->r && agmres->DeflPrecond) {
326: /* Apply the precond-matrix operators */
327: KSP_PCApplyBAorAB(ksp, VEC_V(j-1), VEC_TMP, VEC_TMP_MATOP);
328: /* Then apply deflation as a preconditioner */
329: KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_TMP, VEC_V(j));
330: } else if ((ksp->pc_side == PC_RIGHT) && agmres->r && agmres->DeflPrecond) {
331: KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_V(j-1), VEC_TMP);
332: KSP_PCApplyBAorAB(ksp, VEC_TMP, VEC_V(j), VEC_TMP_MATOP);
333: } else {
334: KSP_PCApplyBAorAB(ksp, VEC_V(j-1), VEC_V(j), VEC_TMP_MATOP);
335: }
336: VecAXPY(VEC_V(j), -Rshift[j-2], VEC_V(j-1));
337: VecAXPY(VEC_V(j), Scale[j-2]*Ishift[j-2]*Ishift[j-2], VEC_V(j-2));
338: #if defined(KSP_AGMRES_NONORM)
339: Scale[j] = 1.0;
340: #else
341: VecNorm(VEC_V(j), NORM_2, &(Scale[j]));
342: Scale[j] = 1.0/Scale[j];
343: #endif
344: agmres->matvecs += 1;
345: j++;
346: }
347: }
348: /* Augment the subspace with the eigenvectors*/
349: while (j <= KspSize) {
350: KSP_PCApplyBAorAB(ksp, agmres->U[j - max_k - 1], VEC_V(j), VEC_TMP_MATOP);
351: #if defined(KSP_AGMRES_NONORM)
352: Scale[j] = 1.0;
353: #else
354: VecScale(VEC_V(j), Scale[j-1]);
355: VecNorm(VEC_V(j), NORM_2, &(Scale[j]));
356: Scale[j] = 1.0/Scale[j];
357: #endif
358: agmres->matvecs += 1;
359: j++;
360: }
361: PetscLogEventEnd(KSP_AGMRESBuildBasis, ksp, 0,0,0);
362: return(0);
363: }
365: /* Form the Hessenberg matrix for the Arnoldi-like relation.
366: * Inputs :
367: * - Shifts values in agmres->Rshift and agmres->Ishift
368: * - RLoc : Triangular matrix from the RODDEC orthogonalization
369: * Outputs :
370: * - H = agmres->hh_origin : The Hessenberg matrix.
371: *
372: * NOTE: Note that the computed Hessenberg matrix is not mathematically equivalent to that in the real Arnoldi process (in KSP GMRES). If it is needed, it can be explicitly formed as H <-- H * RLoc^-1.
373: *
374: */
375: static PetscErrorCode KSPAGMRESBuildHessenberg(KSP ksp)376: {
377: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
378: PetscScalar *Rshift = agmres->Rshift;
379: PetscScalar *Ishift = agmres->Ishift;
380: PetscScalar *Scale = agmres->Scale;
382: PetscInt i = 0, j = 0;
383: const PetscInt max_k = agmres->max_k;
384: PetscInt KspSize = KSPSIZE;
385: PetscInt N = MAXKSPSIZE+1;
388: PetscMemzero(agmres->hh_origin, (N+1)*N*sizeof(PetscScalar));
389: while (j < max_k) {
390: /* Real shifts */
391: if (Ishift[j] == 0) {
392: for (i = 0; i <= j; i++) {
393: *H(i,j) = *RLOC(i,j+1)/Scale[j] + (Rshift[j] * *RLOC(i,j));
394: }
395: *H(j+1,j) = *RLOC(j+1,j+1)/Scale[j];
396: j++;
397: } else if (Ishift[j] > 0) {
398: for (i = 0; i <= j; i++) {
399: *H(i,j) = *RLOC(i,j+1)/Scale[j] + Rshift[j] * *RLOC(i, j);
400: }
401: *H(j+1,j) = *RLOC(j+1, j+1)/Scale[j];
402: j++;
403: for (i = 0; i <= j; i++) {
404: *H(i,j) = (*RLOC(i,j+1) + Rshift[j-1] * *RLOC(i,j) - Scale[j-1] * Ishift[j-1]*Ishift[j-1]* *RLOC(i,j-1));
405: }
406: *H(j,j) = (*RLOC(j,j+1) + Rshift[j-1] * *RLOC(j,j));
407: *H(j+1,j) = *RLOC(j+1,j+1);
408: j++;
409: } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "BAD ORDERING OF THE SHIFTS VALUES IN THE NEWTON BASIS");
410: }
411: for (j = max_k; j< KspSize; j++) { /* take into account the norm of the augmented vectors */
412: for (i = 0; i <= j+1; i++) *H(i,j) = *RLOC(i, j+1)/Scale[j];
413: }
414: return(0);
415: }
417: /*
418: * Form the new approximate solution from the least-square problem
419: *
420: */
422: static PetscErrorCode KSPAGMRESBuildSoln(KSP ksp,PetscInt it)423: {
424: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
426: const PetscInt max_k = agmres->max_k; /* Size of the non-augmented Krylov basis */
427: PetscInt i, j;
428: PetscInt r = agmres->r; /* current number of augmented eigenvectors */
429: PetscBLASInt KspSize;
430: PetscBLASInt lC;
431: PetscBLASInt N;
432: PetscBLASInt ldH = it + 1;
433: PetscBLASInt lwork;
434: PetscBLASInt info, nrhs = 1;
437: PetscBLASIntCast(KSPSIZE,&KspSize);
438: PetscBLASIntCast(4 * (KspSize+1),&lwork);
439: PetscBLASIntCast(KspSize+1,&lC);
440: PetscBLASIntCast(MAXKSPSIZE + 1,&N);
441: PetscBLASIntCast(N + 1,&ldH);
442: /* Save a copy of the Hessenberg matrix */
443: for (j = 0; j < N-1; j++) {
444: for (i = 0; i < N; i++) {
445: *HS(i,j) = *H(i,j);
446: }
447: }
448: /* QR factorize the Hessenberg matrix */
449: #if defined(PETSC_MISSING_LAPACK_GEQRF)
450: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEQRF - Lapack routine is unavailable.");
451: #else
452: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&lC, &KspSize, agmres->hh_origin, &ldH, agmres->tau, agmres->work, &lwork, &info));
453: if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XGEQRF INFO=%d", info);
454: #endif
455: /* Update the right hand side of the least square problem */
456: PetscMemzero(agmres->nrs, N*sizeof(PetscScalar));
458: agmres->nrs[0] = ksp->rnorm;
459: #if defined(PETSC_MISSING_LAPACK_ORMQR)
460: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEQRF - Lapack routine is unavailable.");
461: #else
462: PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", "T", &lC, &nrhs, &KspSize, agmres->hh_origin, &ldH, agmres->tau, agmres->nrs, &N, agmres->work, &lwork, &info));
463: if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XORMQR INFO=%d",info);
464: #endif
465: ksp->rnorm = PetscAbsScalar(agmres->nrs[KspSize]);
466: /* solve the least-square problem */
467: #if defined(PETSC_MISSING_LAPACK_TRTRS)
468: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"TRTRS - Lapack routine is unavailable.");
469: #else
470: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "N", "N", &KspSize, &nrhs, agmres->hh_origin, &ldH, agmres->nrs, &N, &info));
471: if (info) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB,"Error in LAPACK routine XTRTRS INFO=%d",info);
472: #endif
473: /* Accumulate the correction to the solution of the preconditioned problem in VEC_TMP */
474: VecZeroEntries(VEC_TMP);
475: VecMAXPY(VEC_TMP, max_k, agmres->nrs, &VEC_V(0));
476: if (!agmres->DeflPrecond) { VecMAXPY(VEC_TMP, r, &agmres->nrs[max_k], agmres->U); }
478: if ((ksp->pc_side == PC_RIGHT) && agmres->r && agmres->DeflPrecond) {
479: KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_TMP, VEC_TMP_MATOP);
480: VecCopy(VEC_TMP_MATOP, VEC_TMP);
481: }
482: KSPUnwindPreconditioner(ksp, VEC_TMP, VEC_TMP_MATOP);
483: /* add the solution to the previous one */
484: VecAXPY(ksp->vec_sol, 1.0, VEC_TMP);
485: return(0);
486: }
488: /*
489: * Run one cycle of the Newton-basis gmres, possibly augmented with eigenvectors.
490: *
491: * Return residual history if requested.
492: * Input :
493: * - The vector VEC_V(0) is the initia residual
494: * Output :
495: * - the solution vector is in agmres->vec_sol
496: * - itcount : number of inner iterations
497: * - res : the new residual norm
498: .
499: NOTE: Unlike GMRES where the residual norm is available at each (inner) iteration, here it is available at the end of the cycle.
500: */
501: static PetscErrorCode KSPAGMRESCycle(PetscInt *itcount,KSP ksp)502: {
503: KSP_AGMRES *agmres = (KSP_AGMRES*)(ksp->data);
504: PetscReal res;
506: PetscInt KspSize = KSPSIZE;
509: /* check for the convergence */
510: res = ksp->rnorm; /* Norm of the initial residual vector */
511: if (!res) {
512: if (itcount) *itcount = 0;
513: ksp->reason = KSP_CONVERGED_ATOL;
514: PetscInfo(ksp,"Converged due to zero residual norm on entry\n");
515: return(0);
516: }
517: (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
518: /* Build the Krylov basis with Newton polynomials */
519: KSPAGMRESBuildBasis(ksp);
520: /* QR Factorize the basis with RODDEC */
521: KSPAGMRESRoddec(ksp, KspSize+1);
523: /* Recover a (partial) Hessenberg matrix for the Arnoldi-like relation */
524: KSPAGMRESBuildHessenberg(ksp);
525: /* Solve the least square problem and unwind the preconditioner */
526: KSPAGMRESBuildSoln(ksp, KspSize);
528: res = ksp->rnorm;
529: ksp->its += KspSize;
530: agmres->it = KspSize-1;
531: /* Test for the convergence */
532: (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
533: KSPLogResidualHistory(ksp,res);
534: KSPMonitor(ksp,ksp->its,res);
537: *itcount = KspSize;
538: return(0);
539: }
541: static PetscErrorCode KSPSolve_AGMRES(KSP ksp)542: {
544: PetscInt its;
545: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
546: PetscBool guess_zero = ksp->guess_zero;
547: PetscReal res_old, res;
548: PetscInt test;
551: PetscObjectSAWsTakeAccess((PetscObject)ksp);
552: ksp->its = 0;
553: PetscObjectSAWsGrantAccess((PetscObject)ksp);
555: ksp->reason = KSP_CONVERGED_ITERATING;
556: if (!agmres->HasShifts) { /* Compute Shifts for the Newton basis */
557: KSPComputeShifts_DGMRES(ksp);
558: }
559: /* NOTE: At this step, the initial guess is not equal to zero since one cycle of the classical GMRES is performed to compute the shifts */
560: (*ksp->converged)(ksp,0,ksp->rnorm,&ksp->reason,ksp->cnvP);
561: while (!ksp->reason) {
562: KSPInitialResidual(ksp,ksp->vec_sol,VEC_TMP,VEC_TMP_MATOP,VEC_V(0),ksp->vec_rhs);
563: if ((ksp->pc_side == PC_LEFT) && agmres->r && agmres->DeflPrecond) {
564: KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_V(0), VEC_TMP);
565: VecCopy(VEC_TMP, VEC_V(0));
567: agmres->matvecs += 1;
568: }
569: VecNormalize(VEC_V(0),&(ksp->rnorm));
570: KSPCheckNorm(ksp,ksp->rnorm);
571: res_old = ksp->rnorm; /* Record the residual norm to test if deflation is needed */
573: ksp->ops->buildsolution = KSPBuildSolution_AGMRES;
575: KSPAGMRESCycle(&its,ksp);
576: if (ksp->its >= ksp->max_it) {
577: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
578: break;
579: }
580: /* compute the eigenvectors to augment the subspace : use an adaptive strategy */
581: res = ksp->rnorm;
582: if (!ksp->reason && agmres->neig > 0) {
583: test = agmres->max_k * PetscLogReal(ksp->rtol/res) / PetscLogReal(res/res_old); /* estimate the remaining number of steps */
584: if ((test > agmres->smv*(ksp->max_it-ksp->its)) || agmres->force) {
585: if (!agmres->force && ((test > agmres->bgv*(ksp->max_it-ksp->its)) && ((agmres->r + 1) < agmres->max_neig))) {
586: agmres->neig += 1; /* Augment the number of eigenvalues to deflate if the convergence is too slow */
587: }
588: KSPDGMRESComputeDeflationData_DGMRES(ksp,&agmres->neig);
589: }
590: }
591: ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
592: }
593: ksp->guess_zero = guess_zero; /* restore if user has provided nonzero initial guess */
594: return(0);
595: }
597: static PetscErrorCode KSPDestroy_AGMRES(KSP ksp)598: {
600: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
603: PetscFree(agmres->hh_origin);
604: PetscFree(agmres->nrs);
605: PetscFree(agmres->Qloc);
606: PetscFree(agmres->Rloc);
607: PetscFree(agmres->sgn);
608: PetscFree(agmres->tloc);
609: PetscFree(agmres->Rshift);
610: PetscFree(agmres->Ishift);
611: PetscFree(agmres->Scale);
612: PetscFree(agmres->wbufptr);
613: PetscFree(agmres->tau);
614: PetscFree(agmres->work);
615: PetscFree(agmres->temp);
616: PetscFree(agmres->select);
617: PetscFree(agmres->wr);
618: PetscFree(agmres->wi);
619: if (agmres->neig) {
620: VecDestroyVecs(MAXKSPSIZE,&agmres->TmpU);
621: PetscFree(agmres->perm);
622: PetscFree(agmres->MatEigL);
623: PetscFree(agmres->MatEigR);
624: PetscFree(agmres->Q);
625: PetscFree(agmres->Z);
626: PetscFree(agmres->beta);
627: }
628: KSPDestroy_DGMRES(ksp);
629: return(0);
630: }
633: static PetscErrorCode KSPView_AGMRES(KSP ksp,PetscViewer viewer)634: {
635: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
636: const char *cstr = "RODDEC ORTHOGONOLIZATION";
637: char ritzvec[25];
639: PetscBool iascii,isstring;
640: #if defined(KSP_AGMRES_NONORM)
641: const char *Nstr = "SCALING FACTORS : NO";
642: #else
643: const char *Nstr = "SCALING FACTORS : YES";
644: #endif
647: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
648: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
650: if (iascii) {
651: PetscViewerASCIIPrintf(viewer, " restart=%d using %s\n", agmres->max_k, cstr);
652: PetscViewerASCIIPrintf(viewer, " %s\n", Nstr);
653: PetscViewerASCIIPrintf(viewer, " Number of matvecs : %D\n", agmres->matvecs);
654: if (agmres->force) {PetscViewerASCIIPrintf (viewer, " Adaptive strategy is used: FALSE\n");}
655: else PetscViewerASCIIPrintf(viewer, " Adaptive strategy is used: TRUE\n");
656: if (agmres->DeflPrecond) {
657: ierr=PetscViewerASCIIPrintf(viewer, " STRATEGY OF DEFLATION: PRECONDITIONER \n");
658: ierr=PetscViewerASCIIPrintf(viewer, " Frequency of extracted eigenvalues = %D\n", agmres->neig);
659: ierr=PetscViewerASCIIPrintf(viewer, " Total number of extracted eigenvalues = %D\n", agmres->r);
660: ierr=PetscViewerASCIIPrintf(viewer, " Maximum number of eigenvalues set to be extracted = %D\n", agmres->max_neig);
661: } else {
662: if (agmres->ritz) sprintf(ritzvec, "Ritz vectors");
663: else sprintf(ritzvec, "Harmonic Ritz vectors");
664: PetscViewerASCIIPrintf(viewer, " STRATEGY OF DEFLATION: AUGMENT\n");
665: PetscViewerASCIIPrintf(viewer," augmented vectors %d at frequency %d with %s\n", agmres->r, agmres->neig, ritzvec);
666: }
667: ierr=PetscViewerASCIIPrintf(viewer, " Minimum relaxation parameter for the adaptive strategy(smv) = %g\n", agmres->smv);
668: ierr=PetscViewerASCIIPrintf(viewer, " Maximum relaxation parameter for the adaptive strategy(bgv) = %g\n", agmres->bgv);
669: } else if (isstring) {
670: PetscViewerStringSPrintf(viewer,"%s restart %D",cstr,agmres->max_k);
671: }
672: return(0);
673: }
675: static PetscErrorCode KSPSetFromOptions_AGMRES(PetscOptionItems *PetscOptionsObject,KSP ksp)676: {
678: PetscInt neig;
679: KSP_AGMRES *agmres = (KSP_AGMRES*)ksp->data;
680: PetscBool flg;
683: KSPSetFromOptions_DGMRES(PetscOptionsObject,ksp); /* Set common options from DGMRES and GMRES */
684: PetscOptionsHead(PetscOptionsObject,"KSP AGMRES Options");
685: PetscOptionsInt("-ksp_agmres_eigen", "Number of eigenvalues to deflate", "KSPDGMRESSetEigen", agmres->neig, &neig, &flg);
686: if (flg) {
687: KSPDGMRESSetEigen_DGMRES(ksp, neig);
688: agmres->r = 0;
689: } else agmres->neig = 0;
690: PetscOptionsInt("-ksp_agmres_maxeigen", "Maximum number of eigenvalues to deflate", "KSPDGMRESSetMaxEigen", agmres->max_neig, &neig, &flg);
691: if (flg) agmres->max_neig = neig+EIG_OFFSET;
692: else agmres->max_neig = agmres->neig+EIG_OFFSET;
693: PetscOptionsBool("-ksp_agmres_DeflPrecond", "Determine if the deflation should be applied as a preconditioner -- similar to KSP DGMRES", "KSPGMRESDeflPrecond",agmres->DeflPrecond,&agmres->DeflPrecond,NULL);
694: PetscOptionsBool("-ksp_agmres_ritz", "Compute the Ritz vectors instead of the Harmonic Ritz vectors ", "KSPGMRESHarmonic",agmres->ritz,&agmres->ritz ,&flg);
695: PetscOptionsReal("-ksp_agmres_MinRatio", "Relaxation parameter in the adaptive strategy; smallest multiple of the remaining number of steps allowed", "KSPGMRESSetMinRatio", agmres->smv, &agmres->smv, NULL);
696: PetscOptionsReal("-ksp_agmres_MaxRatio", "Relaxation parameter in the adaptive strategy; Largest multiple of the remaining number of steps allowed", "KSPGMRESSetMaxRatio",agmres->bgv,&agmres->bgv, &flg);
697: PetscOptionsTail();
698: return(0);
699: }
702: /*MC
703: KSPAGMRES - Newton basis GMRES implementation with adaptive augmented eigenvectors
705: The techniques used are best described in [1]. The contribution of this work is that it combines many of the previous work to reduce the amount of MPI messages and improve the robustness of the global approach by using deflation techniques. It has been successfully applied to a class of real and industrial problems. Please see [1] for numerical experiments and [2] for a description of these problems.
706: There are many ongoing work that aim at avoiding (or minimizing) the communication in Krylov subspace methods. This code can be used as an experimental framework to combine several techniques in the particular case of GMRES. For instance, the computation of the shifts can be improved with techniques described in [3]. The orthogonalization technique can be replaced by TSQR [4]. The generation of the basis can be done using s-steps approaches[5].
709: Options Database Keys:
710: + -ksp_gmres_restart <restart> - the number of Krylov directions
711: . -ksp_gmres_krylov_monitor - plot the Krylov space generated
712: . -ksp_agmres_eigen <neig> - Number of eigenvalues to deflate (Number of vectors to augment)
713: . -ksp_agmres_maxeigen <max_neig> - Maximum number of eigenvalues to deflate
714: . -ksp_agmres_MinRatio <1> - Relaxation parameter in the adaptive strategy; smallest multiple of the remaining number of steps allowed
715: . -ksp_agmres_MaxRatio <1> - Relaxation parameter in the adaptive strategy; Largest multiple of the remaining number of steps allowed
716: . --ksp_agmres_DeflPrecond Apply deflation as a preconditioner, this is similar to DGMRES but it rather builds a Newton basis. This is an experimental option.
717: . -ksp_dgmres_force <0, 1> - Force the deflation at each restart.
718: . - There are many experimental parameters. Run with -help option to see the whole list
720: Level: beginner
722: Notes:
723: Left and right preconditioning are supported, but not symmetric preconditioning. Complex arithmetic is not supported
725: Developer Notes:
726: This object is subclassed off of KSPDGMRES728: Contributed by Desire NUENTSA WAKAM, INRIA <desire.nuentsa_wakam@inria.fr>
729: Inputs from Guy Atenekeng <atenekeng@yahoo.com> and R.B. Sidje <roger.b.sidje@ua.edu>
731: References :
732: + [1] D. Nuentsa Wakam and J. Erhel, Parallelism and robustness in GMRES with the Newton basis and the deflated restarting. Research report INRIA RR-7787, November 2011,http://hal.inria.fr/inria-00638247/en, in revision for ETNA.
733: . [2] D. NUENTSA WAKAM and F. PACULL, Memory Efficient Hybrid Algebraic Solvers for Linear Systems Arising from Compressible Flows, Computers and Fluids, In Press, http://dx.doi.org/10.1016/j.compfluid.2012.03.023
734: . [3] B. Philippe and L. Reichel, On the generation of Krylov subspace bases, Applied Numerical
735: Mathematics, 62(9), pp. 1171-1186, 2012
736: . [4] J. Demmel, L. Grigori, M. F. Hoemmen, and J. Langou, Communication-optimal parallel and sequential QR and LU factorizations, SIAM journal on Scientific Computing, 34(1), A206-A239, 2012
737: . [5] M. Mohiyuddin, M. Hoemmen, J. Demmel, and K. Yelick, Minimizing communication in sparse matrix solvers, in SC '09: Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, New York, NY, USA, 2009, ACM, pp. 1154-1171.
738: . Sidje, Roger B. Alternatives for parallel Krylov subspace basis computation. Numer. Linear Algebra Appl. 4 (1997), no. 4, 305-331
740: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPDGMRES, KSPPGMRES,
741: KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
742: KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
743: KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPSetPCSide()
744: M*/
746: PETSC_EXTERN PetscErrorCode KSPCreate_AGMRES(KSP ksp)747: {
748: KSP_AGMRES *agmres;
752: PetscNewLog(ksp,&agmres);
753: ksp->data = (void*)agmres;
755: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
756: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
757: ksp->ops->buildsolution = KSPBuildSolution_AGMRES;
758: ksp->ops->setup = KSPSetUp_AGMRES;
759: ksp->ops->solve = KSPSolve_AGMRES;
760: ksp->ops->destroy = KSPDestroy_AGMRES;
761: ksp->ops->view = KSPView_AGMRES;
762: ksp->ops->setfromoptions = KSPSetFromOptions_AGMRES;
763: ksp->guess_zero = PETSC_TRUE;
764: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
765: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
768: PetscObjectComposeFunction((PetscObject) ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);
769: PetscObjectComposeFunction((PetscObject) ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);
770: PetscObjectComposeFunction((PetscObject) ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);
771: PetscObjectComposeFunction((PetscObject) ksp,"KSPGMRESSetHapTol_C",KSPGMRESSetHapTol_GMRES);
772: PetscObjectComposeFunction((PetscObject) ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);
773: /* -- New functions defined in DGMRES -- */
774: ierr=PetscObjectComposeFunction((PetscObject) ksp, "KSPDGMRESSetEigen_C",KSPDGMRESSetEigen_DGMRES);
775: ierr=PetscObjectComposeFunction((PetscObject) ksp, "KSPDGMRESComputeSchurForm_C",KSPDGMRESComputeSchurForm_DGMRES);
776: ierr=PetscObjectComposeFunction((PetscObject) ksp, "KSPDGMRESComputeDeflationData_C",KSPDGMRESComputeDeflationData_DGMRES);
777: ierr=PetscObjectComposeFunction((PetscObject) ksp, "KSPDGMRESApplyDeflation_C",KSPDGMRESApplyDeflation_DGMRES);
779: PetscLogEventRegister("AGMRESCompDefl", KSP_CLASSID, &KSP_AGMRESComputeDeflationData);
780: PetscLogEventRegister("AGMRESBuildBasis", KSP_CLASSID, &KSP_AGMRESBuildBasis);
781: PetscLogEventRegister("AGMRESCompShifts", KSP_CLASSID, &KSP_AGMRESComputeShifts);
782: PetscLogEventRegister("AGMRESOrthog", KSP_CLASSID, &KSP_AGMRESRoddec);
784: agmres->haptol = 1.0e-30;
785: agmres->q_preallocate = 0;
786: agmres->delta_allocate = AGMRES_DELTA_DIRECTIONS;
787: agmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
788: agmres->nrs = 0;
789: agmres->sol_temp = 0;
790: agmres->max_k = AGMRES_DEFAULT_MAXK;
791: agmres->Rsvd = 0;
792: agmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
793: agmres->orthogwork = 0;
795: /* Default values for the deflation */
796: agmres->r = 0;
797: agmres->neig = 0;
798: agmres->max_neig = 0;
799: agmres->lambdaN = 0.0;
800: agmres->smv = SMV;
801: agmres->bgv = 1;
802: agmres->force = PETSC_FALSE;
803: agmres->matvecs = 0;
804: agmres->improve = PETSC_FALSE;
805: agmres->HasShifts = PETSC_FALSE;
806: agmres->r = 0;
807: agmres->HasSchur = PETSC_FALSE;
808: agmres->DeflPrecond = PETSC_FALSE;
809: PetscObjectGetNewTag((PetscObject)ksp,&agmres->tag);
810: return(0);
811: }
813: /* LocalWords: iascii
814: */