Actual source code: cgne.c


  2: /*
  3:        cgimpl.h defines the simple data structured used to store information
  4:     related to the type of matrix (e.g. complex symmetric) being solved and
  5:     data used during the optional Lanczo process used to compute eigenvalues
  6: */
  7: #include <../src/ksp/ksp/impls/cg/cgimpl.h>
  8: extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP, PetscReal *, PetscReal *);
  9: extern PetscErrorCode KSPComputeEigenvalues_CG(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);

 11: static PetscErrorCode KSPCGSetType_CGNE(KSP ksp, KSPCGType type)
 12: {
 13:   KSP_CG *cg = (KSP_CG *)ksp->data;

 15:   cg->type = type;
 16:   return 0;
 17: }

 19: /*
 20:      KSPSetUp_CGNE - Sets up the workspace needed by the CGNE method.

 22:      IDENTICAL TO THE CG ONE EXCEPT for one extra work vector!
 23: */
 24: static PetscErrorCode KSPSetUp_CGNE(KSP ksp)
 25: {
 26:   KSP_CG  *cgP   = (KSP_CG *)ksp->data;
 27:   PetscInt maxit = ksp->max_it;

 29:   /* get work vectors needed by CGNE */
 30:   KSPSetWorkVecs(ksp, 4);

 32:   /*
 33:      If user requested computations of eigenvalues then allocate work
 34:      work space needed
 35:   */
 36:   if (ksp->calc_sings) {
 37:     /* get space to store tridiagonal matrix for Lanczos */
 38:     PetscMalloc4(maxit, &cgP->e, maxit, &cgP->d, maxit, &cgP->ee, maxit, &cgP->dd);

 40:     ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
 41:     ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_CG;
 42:   }
 43:   return 0;
 44: }

 46: /*
 47:        KSPSolve_CGNE - This routine actually applies the conjugate gradient
 48:     method

 50:    Input Parameter:
 51: .     ksp - the Krylov space object that was set to use conjugate gradient, by, for
 52:             example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG);

 54:     Virtually identical to the KSPSolve_CG, it should definitely reuse the same code.

 56: */
 57: static PetscErrorCode KSPSolve_CGNE(KSP ksp)
 58: {
 59:   PetscInt    i, stored_max_it, eigs;
 60:   PetscScalar dpi, a = 1.0, beta, betaold = 1.0, b = 0, *e = NULL, *d = NULL;
 61:   PetscReal   dp = 0.0;
 62:   Vec         X, B, Z, R, P, T;
 63:   KSP_CG     *cg;
 64:   Mat         Amat, Pmat;
 65:   PetscBool   diagonalscale, transpose_pc;

 67:   PCGetDiagonalScale(ksp->pc, &diagonalscale);
 69:   PCApplyTransposeExists(ksp->pc, &transpose_pc);

 71:   cg            = (KSP_CG *)ksp->data;
 72:   eigs          = ksp->calc_sings;
 73:   stored_max_it = ksp->max_it;
 74:   X             = ksp->vec_sol;
 75:   B             = ksp->vec_rhs;
 76:   R             = ksp->work[0];
 77:   Z             = ksp->work[1];
 78:   P             = ksp->work[2];
 79:   T             = ksp->work[3];

 81: #define VecXDot(x, y, a) (((cg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x, y, a) : VecTDot(x, y, a))

 83:   if (eigs) {
 84:     e    = cg->e;
 85:     d    = cg->d;
 86:     e[0] = 0.0;
 87:   }
 88:   PCGetOperators(ksp->pc, &Amat, &Pmat);

 90:   ksp->its = 0;
 91:   KSP_MatMultTranspose(ksp, Amat, B, T);
 92:   if (!ksp->guess_zero) {
 93:     KSP_MatMult(ksp, Amat, X, P);
 94:     KSP_MatMultTranspose(ksp, Amat, P, R);
 95:     VecAYPX(R, -1.0, T);
 96:   } else {
 97:     VecCopy(T, R); /*     r <- b (x is 0) */
 98:   }
 99:   if (transpose_pc) {
100:     KSP_PCApplyTranspose(ksp, R, T);
101:   } else {
102:     KSP_PCApply(ksp, R, T);
103:   }
104:   KSP_PCApply(ksp, T, Z);

106:   if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
107:     VecNorm(Z, NORM_2, &dp); /*    dp <- z'*z       */
108:   } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
109:     VecNorm(R, NORM_2, &dp); /*    dp <- r'*r       */
110:   } else if (ksp->normtype == KSP_NORM_NATURAL) {
111:     VecXDot(Z, R, &beta);
112:     KSPCheckDot(ksp, beta);
113:     dp = PetscSqrtReal(PetscAbsScalar(beta));
114:   } else dp = 0.0;
115:   KSPLogResidualHistory(ksp, dp);
116:   KSPMonitor(ksp, 0, dp);
117:   ksp->rnorm = dp;
118:   (*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP); /* test for convergence */
119:   if (ksp->reason) return 0;

121:   i = 0;
122:   do {
123:     ksp->its = i + 1;
124:     VecXDot(Z, R, &beta); /*     beta <- r'z     */
125:     KSPCheckDot(ksp, beta);
126:     if (beta == 0.0) {
127:       ksp->reason = KSP_CONVERGED_ATOL;
128:       PetscInfo(ksp, "converged due to beta = 0\n");
129:       break;
130: #if !defined(PETSC_USE_COMPLEX)
131:     } else if (beta < 0.0) {
132:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
133:       PetscInfo(ksp, "diverging due to indefinite preconditioner\n");
134:       break;
135: #endif
136:     }
137:     if (!i) {
138:       VecCopy(Z, P); /*     p <- z          */
139:       b = 0.0;
140:     } else {
141:       b = beta / betaold;
142:       if (eigs) {
144:         e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
145:       }
146:       VecAYPX(P, b, Z); /*     p <- z + b* p   */
147:     }
148:     betaold = beta;
149:     KSP_MatMult(ksp, Amat, P, T);
150:     KSP_MatMultTranspose(ksp, Amat, T, Z);
151:     VecXDot(P, Z, &dpi); /*     dpi <- z'p      */
152:     KSPCheckDot(ksp, dpi);
153:     a = beta / dpi; /*     a = beta/p'z    */
154:     if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
155:     VecAXPY(X, a, P);  /*     x <- x + ap     */
156:     VecAXPY(R, -a, Z); /*     r <- r - az     */
157:     if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
158:       if (transpose_pc) {
159:         KSP_PCApplyTranspose(ksp, R, T);
160:       } else {
161:         KSP_PCApply(ksp, R, T);
162:       }
163:       KSP_PCApply(ksp, T, Z);
164:       VecNorm(Z, NORM_2, &dp); /*    dp <- z'*z       */
165:     } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
166:       VecNorm(R, NORM_2, &dp);
167:     } else if (ksp->normtype == KSP_NORM_NATURAL) {
168:       dp = PetscSqrtReal(PetscAbsScalar(beta));
169:     } else dp = 0.0;
170:     ksp->rnorm = dp;
171:     KSPLogResidualHistory(ksp, dp);
172:     KSPMonitor(ksp, i + 1, dp);
173:     (*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP);
174:     if (ksp->reason) break;
175:     if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
176:       if (transpose_pc) {
177:         KSP_PCApplyTranspose(ksp, R, T);
178:       } else {
179:         KSP_PCApply(ksp, R, T);
180:       }
181:       KSP_PCApply(ksp, T, Z);
182:     }
183:     i++;
184:   } while (i < ksp->max_it);
185:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
186:   return 0;
187: }

189: /*
190:     KSPCreate_CGNE - Creates the data structure for the Krylov method CGNE and sets the
191:        function pointers for all the routines it needs to call (KSPSolve_CGNE() etc)

193:     It must be labeled as PETSC_EXTERN to be dynamically linkable in C++
194: */

196: /*MC
197:      KSPCGNE - Applies the preconditioned conjugate gradient method to the normal equations
198:           without explicitly forming A^t*A

200:    Options Database Keys:
201: .   -ksp_cg_type <Hermitian or symmetric - (for complex matrices only) indicates the matrix is Hermitian or symmetric

203:    Level: beginner

205:    Notes:
206:    Eigenvalue computation routines including `KSPSetComputeEigenvalues()` and `KSPComputeEigenvalues()` will return information about the
207:     spectrum of A^t*A, rather than A.

209:    `KSPCGNE` is a general-purpose non-symmetric method. It works well when the singular values are much better behaved than
210:    eigenvalues. A unitary matrix is a classic example where `KSPCGNE` converges in one iteration, but `KSPGMRES` and `KSPCGS` need N
211:    iterations, see [1]. If you intend to solve least squares problems, use `KSPLSQR`.

213:    This is NOT a different algorithm than used with `KSPCG`, it merely uses that algorithm with the
214:    matrix defined by A^t*A and preconditioner defined by B^t*B where B is the preconditioner for A.

216:    This method requires that one be able to apply the transpose of the preconditioner and operator
217:    as well as the operator and preconditioner. If the transpose of the preconditioner is not available then
218:    the preconditioner is used in its place so one ends up preconditioning A'A with B B. Seems odd?

220:    This only supports left preconditioning.

222:    Reference:
223: .   [1] -  Nachtigal, Reddy, and Trefethen, "How fast are nonsymmetric matrix iterations", 1992

225:    Developer Note:
226:    This object is subclassed off of `KSPCG`

228: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, 'KSPCG', `KSPLSQR', 'KSPCGLS`,
229:           `KSPCGSetType()`, `KSPBICG`, `KSPSetComputeEigenvalues()`, `KSPComputeEigenvalues()`
230: M*/

232: PETSC_EXTERN PetscErrorCode KSPCreate_CGNE(KSP ksp)
233: {
234:   KSP_CG *cg;

236:   PetscNew(&cg);
237: #if !defined(PETSC_USE_COMPLEX)
238:   cg->type = KSP_CG_SYMMETRIC;
239: #else
240:   cg->type = KSP_CG_HERMITIAN;
241: #endif
242:   ksp->data = (void *)cg;
243:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
244:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);
245:   KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
246:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);

248:   /*
249:        Sets the functions that are associated with this data structure
250:        (in C++ this is the same as defining virtual functions)
251:   */
252:   ksp->ops->setup          = KSPSetUp_CGNE;
253:   ksp->ops->solve          = KSPSolve_CGNE;
254:   ksp->ops->destroy        = KSPDestroy_CG;
255:   ksp->ops->view           = KSPView_CG;
256:   ksp->ops->setfromoptions = KSPSetFromOptions_CG;
257:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
258:   ksp->ops->buildresidual  = KSPBuildResidualDefault;

260:   /*
261:       Attach the function KSPCGSetType_CGNE() to this object. The routine
262:       KSPCGSetType() checks for this attached function and calls it if it finds
263:       it. (Sort of like a dynamic member function that can be added at run time
264:   */
265:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", KSPCGSetType_CGNE);
266:   return 0;
267: }