Actual source code: lcd.c
1: #include <../src/ksp/ksp/impls/lcd/lcdimpl.h>
3: static PetscErrorCode KSPSetUp_LCD(KSP ksp)
4: {
5: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
6: PetscInt restart = lcd->restart;
8: PetscFunctionBegin;
9: /* get work vectors needed by LCD */
10: PetscCall(KSPSetWorkVecs(ksp, 2));
12: PetscCall(VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->P));
13: PetscCall(VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->Q));
14: PetscFunctionReturn(PETSC_SUCCESS);
15: }
17: /* KSPSolve_LCD - This routine actually applies the left conjugate
18: direction method
20: Input Parameter:
21: . ksp - the Krylov space object that was set to use LCD, by, for
22: example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPLCD);
24: Output Parameter:
25: . its - number of iterations used
27: */
28: static PetscErrorCode KSPSolve_LCD(KSP ksp)
29: {
30: PetscInt it, j, max_k;
31: PetscScalar alfa, beta, num, den, mone;
32: PetscReal rnorm = 0.0;
33: Vec X, B, R, Z;
34: KSP_LCD *lcd;
35: Mat Amat, Pmat;
36: PetscBool diagonalscale;
38: PetscFunctionBegin;
39: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
40: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
42: lcd = (KSP_LCD *)ksp->data;
43: X = ksp->vec_sol;
44: B = ksp->vec_rhs;
45: R = ksp->work[0];
46: Z = ksp->work[1];
47: max_k = lcd->restart;
48: mone = -1;
50: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
52: ksp->its = 0;
53: if (!ksp->guess_zero) {
54: PetscCall(KSP_MatMult(ksp, Amat, X, Z)); /* z <- b - Ax */
55: PetscCall(VecAYPX(Z, mone, B));
56: } else {
57: PetscCall(VecCopy(B, Z)); /* z <- b (x is 0) */
58: }
60: PetscCall(KSP_PCApply(ksp, Z, R)); /* r <- M^-1z */
61: if (ksp->normtype != KSP_NORM_NONE) {
62: PetscCall(VecNorm(R, NORM_2, &rnorm));
63: KSPCheckNorm(ksp, rnorm);
64: }
65: PetscCall(KSPLogResidualHistory(ksp, rnorm));
66: PetscCall(KSPMonitor(ksp, 0, rnorm));
67: ksp->rnorm = rnorm;
69: /* test for convergence */
70: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
71: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
73: PetscCall(VecCopy(R, lcd->P[0]));
75: while (!ksp->reason && ksp->its < ksp->max_it) {
76: it = 0;
77: PetscCall(KSP_MatMult(ksp, Amat, lcd->P[it], Z));
78: PetscCall(KSP_PCApply(ksp, Z, lcd->Q[it]));
80: while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
81: ksp->its++;
82: PetscCall(VecDot(lcd->P[it], R, &num));
83: PetscCall(VecDot(lcd->P[it], lcd->Q[it], &den));
84: KSPCheckDot(ksp, den);
85: alfa = num / den;
86: PetscCall(VecAXPY(X, alfa, lcd->P[it]));
87: PetscCall(VecAXPY(R, -alfa, lcd->Q[it]));
88: if (ksp->normtype != KSP_NORM_NONE) {
89: PetscCall(VecNorm(R, NORM_2, &rnorm));
90: KSPCheckNorm(ksp, rnorm);
91: }
93: ksp->rnorm = rnorm;
94: PetscCall(KSPLogResidualHistory(ksp, rnorm));
95: PetscCall(KSPMonitor(ksp, ksp->its, rnorm));
96: PetscCall((*ksp->converged)(ksp, ksp->its, rnorm, &ksp->reason, ksp->cnvP));
98: if (ksp->reason) break;
100: PetscCall(VecCopy(R, lcd->P[it + 1]));
101: PetscCall(KSP_MatMult(ksp, Amat, lcd->P[it + 1], Z));
102: PetscCall(KSP_PCApply(ksp, Z, lcd->Q[it + 1]));
104: for (j = 0; j <= it; j++) {
105: PetscCall(VecDot(lcd->P[j], lcd->Q[it + 1], &num));
106: KSPCheckDot(ksp, num);
107: PetscCall(VecDot(lcd->P[j], lcd->Q[j], &den));
108: beta = -num / den;
109: PetscCall(VecAXPY(lcd->P[it + 1], beta, lcd->P[j]));
110: PetscCall(VecAXPY(lcd->Q[it + 1], beta, lcd->Q[j]));
111: }
112: it++;
113: }
114: PetscCall(VecCopy(lcd->P[it], lcd->P[0]));
115: }
116: if (ksp->its >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
117: PetscCall(VecCopy(X, ksp->vec_sol));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
120: /*
121: KSPDestroy_LCD - Frees all memory space used by the Krylov method
123: */
124: static PetscErrorCode KSPReset_LCD(KSP ksp)
125: {
126: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
128: PetscFunctionBegin;
129: if (lcd->P) PetscCall(VecDestroyVecs(lcd->restart + 1, &lcd->P));
130: if (lcd->Q) PetscCall(VecDestroyVecs(lcd->restart + 1, &lcd->Q));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: static PetscErrorCode KSPDestroy_LCD(KSP ksp)
135: {
136: PetscFunctionBegin;
137: PetscCall(KSPReset_LCD(ksp));
138: PetscCall(PetscFree(ksp->data));
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /*
143: KSPView_LCD - Prints information about the current Krylov method being used
145: Currently this only prints information to a file (or stdout) about the
146: symmetry of the problem. If your Krylov method has special options or
147: flags that information should be printed here.
149: */
150: static PetscErrorCode KSPView_LCD(KSP ksp, PetscViewer viewer)
151: {
152: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
153: PetscBool iascii;
155: PetscFunctionBegin;
156: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
157: if (iascii) {
158: PetscCall(PetscViewerASCIIPrintf(viewer, " restart=%" PetscInt_FMT "\n", lcd->restart));
159: PetscCall(PetscViewerASCIIPrintf(viewer, " happy breakdown tolerance %g\n", (double)lcd->haptol));
160: }
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*
165: KSPSetFromOptions_LCD - Checks the options database for options related to the
166: LCD method.
167: */
168: static PetscErrorCode KSPSetFromOptions_LCD(KSP ksp, PetscOptionItems *PetscOptionsObject)
169: {
170: PetscBool flg;
171: KSP_LCD *lcd = (KSP_LCD *)ksp->data;
173: PetscFunctionBegin;
174: PetscOptionsHeadBegin(PetscOptionsObject, "KSP LCD options");
175: PetscCall(PetscOptionsBoundedInt("-ksp_lcd_restart", "Number of vectors conjugate", "KSPLCDSetRestart", lcd->restart, &lcd->restart, &flg, 1));
176: PetscCall(PetscOptionsBoundedReal("-ksp_lcd_haptol", "Tolerance for exact convergence (happy ending)", "KSPLCDSetHapTol", lcd->haptol, &lcd->haptol, &flg, 0.0));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*MC
181: KSPLCD - Implements the LCD (left conjugate direction) method
183: Options Database Keys:
184: + -ksp_lcd_restart - number of vectors conjugate
185: - -ksp_lcd_haptol - tolerance for exact convergence (happy ending)
187: Level: beginner
189: Notes:
190: Support only for left preconditioning
192: See {cite}`yuan2004semi`, {cite}`dai2004study`, {cite}`catabriga2004evaluating`, and {cite}`catabriga2006performance`
194: Contributed by:
195: Lucia Catabriga <luciac@ices.utexas.edu>
197: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`,
198: `KSPCGSetType()`, `KSPLCDSetRestart()`, `KSPLCDSetHapTol()`
199: M*/
201: PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
202: {
203: KSP_LCD *lcd;
205: PetscFunctionBegin;
206: PetscCall(PetscNew(&lcd));
207: ksp->data = (void *)lcd;
208: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
209: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
210: lcd->restart = 30;
211: lcd->haptol = 1.0e-30;
213: /*
214: Sets the functions that are associated with this data structure
215: (in C++ this is the same as defining virtual functions)
216: */
217: ksp->ops->setup = KSPSetUp_LCD;
218: ksp->ops->solve = KSPSolve_LCD;
219: ksp->ops->reset = KSPReset_LCD;
220: ksp->ops->destroy = KSPDestroy_LCD;
221: ksp->ops->view = KSPView_LCD;
222: ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
223: ksp->ops->buildsolution = KSPBuildSolutionDefault;
224: ksp->ops->buildresidual = KSPBuildResidualDefault;
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }