Lines Matching refs:ksp
33 must call KSPCGSetType(ksp,KSP_CG_SYMMETRIC) or use the option
43 #include <../src/ksp/ksp/impls/cg/cgimpl.h> /*I "petscksp.h" I*/
47 static PetscErrorCode KSPCGSetObjectiveTarget_CG(KSP ksp, PetscReal obj_min)
49 KSP_CG *cg = (KSP_CG *)ksp->data;
56 static PetscErrorCode KSPCGSetRadius_CG(KSP ksp, PetscReal radius)
58 KSP_CG *cg = (KSP_CG *)ksp->data;
65 static PetscErrorCode KSPCGGetObjFcn_CG(KSP ksp, PetscReal *obj)
67 KSP_CG *cg = (KSP_CG *)ksp->data;
80 static PetscErrorCode KSPSetUp_CG(KSP ksp)
82 KSP_CG *cgP = (KSP_CG *)ksp->data;
83 PetscInt maxit = ksp->max_it, nwork = 3;
88 PetscCall(KSPSetWorkVecs(ksp, nwork));
94 if (ksp->calc_sings) {
98 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
99 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG;
116 . ksp - the Krylov space object that was set to use conjugate gradient, by, for
117 example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG);
119 static PetscErrorCode KSPSolve_CG(KSP ksp)
131 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
132 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
134 cg = (KSP_CG *)ksp->data;
135 eigs = ksp->calc_sings;
136 stored_max_it = ksp->max_it;
137 X = ksp->vec_sol;
138 B = ksp->vec_rhs;
139 R = ksp->work[0];
140 Z = ksp->work[1];
141 P = ksp->work[2];
150 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
152 ksp->its = 0;
153 if (!ksp->guess_zero) {
154 PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - Ax */
166 PetscCall(VecFlag(R, ksp->reason == KSP_DIVERGED_PC_FAILED));
168 switch (ksp->normtype) {
170 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
172 KSPCheckNorm(ksp, dp);
176 KSPCheckNorm(ksp, dp);
179 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
181 KSPCheckDot(ksp, beta);
188 SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
199 if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " obj %g\n", ksp->its, (double)cg->obj));
200 PetscCall(KSPLogResidualHistory(ksp, dp));
201 PetscCall(KSPMonitor(ksp, ksp->its, dp));
202 ksp->rnorm = dp;
204 PetscCall((*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
206 if (!ksp->reason && testobj && cg->obj <= cg->obj_min) {
207 PetscCall(PetscInfo(ksp, "converged to objective target minimum\n"));
208 ksp->reason = KSP_CONVERGED_ATOL;
211 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
213 if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
214 if (ksp->normtype != KSP_NORM_NATURAL) {
216 KSPCheckDot(ksp, beta);
221 ksp->its = i + 1;
223 ksp->reason = KSP_CONVERGED_ATOL;
224 PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
228 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite preconditioner, beta %g, betaold %g", (double)PetscRealPart(beta), (double)PetscRealPart(betaold));
229 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
230 PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
240 if (!ksp->guess_zero) PetscCall(VecDotRealPart(X, P, &dMp));
246 PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
257 PetscCall(KSP_MatMult(ksp, Amat, P, W)); /* w <- Ap */
259 KSPCheckDot(ksp, dpi);
278 if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " N obj %g\n", i + 1, (double)cg->obj));
279 if (ksp->converged_neg_curve) {
280 PetscCall(PetscInfo(ksp, "converged due to negative curvature: %g\n", (double)(PetscRealPart(dpi))));
281 ksp->reason = KSP_CONVERGED_NEG_CURVE;
283 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite matrix, dpi %g, dpiold %g", (double)PetscRealPart(dpi), (double)PetscRealPart(dpiold));
284 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
285 PetscCall(PetscInfo(ksp, "diverging due to indefinite matrix\n"));
294 ksp->reason = KSP_CONVERGED_STEP_LENGTH;
295 PetscCall(PetscInfo(ksp, "converged to the trust region radius %g\n", (double)cg->radius));
301 if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " R obj %g\n", i + 1, (double)cg->obj));
307 if (ksp->normtype == KSP_NORM_PRECONDITIONED && ksp->chknorm < i + 2) {
308 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
310 KSPCheckNorm(ksp, dp);
311 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED && ksp->chknorm < i + 2) {
313 KSPCheckNorm(ksp, dp);
314 } else if (ksp->normtype == KSP_NORM_NATURAL) {
315 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
317 KSPCheckDot(ksp, beta);
323 if (testobj) PetscCall(PetscInfo(ksp, "it %" PetscInt_FMT " obj %g\n", i + 1, (double)cg->obj));
325 ksp->rnorm = dp;
326 PetscCall(KSPLogResidualHistory(ksp, dp));
327 PetscCall(KSPMonitor(ksp, i + 1, dp));
328 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
330 if (!ksp->reason && testobj && cg->obj <= cg->obj_min) {
331 PetscCall(PetscInfo(ksp, "converged to objective target minimum\n"));
332 ksp->reason = KSP_CONVERGED_ATOL;
335 if (ksp->reason) break;
342 if ((ksp->normtype != KSP_NORM_PRECONDITIONED && ksp->normtype != KSP_NORM_NATURAL) || ksp->chknorm >= i + 2) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
343 if (ksp->normtype != KSP_NORM_NATURAL || ksp->chknorm >= i + 2) {
345 KSPCheckDot(ksp, beta);
349 } while (i < ksp->max_it);
350 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
364 static PetscErrorCode KSPSolve_CG_SingleReduction(KSP ksp)
375 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
376 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
378 cg = (KSP_CG *)ksp->data;
379 eigs = ksp->calc_sings;
380 stored_max_it = ksp->max_it;
381 X = ksp->vec_sol;
382 B = ksp->vec_rhs;
383 R = ksp->work[0];
384 Z = ksp->work[1];
385 P = ksp->work[2];
386 S = ksp->work[3];
387 W = ksp->work[4];
394 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
396 ksp->its = 0;
397 if (!ksp->guess_zero) {
398 PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - Ax */
404 switch (ksp->normtype) {
406 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
408 KSPCheckNorm(ksp, dp);
412 KSPCheckNorm(ksp, dp);
415 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
416 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
419 KSPCheckDot(ksp, beta);
426 SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
428 PetscCall(KSPLogResidualHistory(ksp, dp));
429 PetscCall(KSPMonitor(ksp, 0, dp));
430 ksp->rnorm = dp;
432 PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
433 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
435 if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
436 if (ksp->normtype != KSP_NORM_NATURAL) {
437 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
440 KSPCheckDot(ksp, beta);
445 ksp->its = i + 1;
447 ksp->reason = KSP_CONVERGED_ATOL;
448 PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
452 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite preconditioner");
453 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
454 PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
464 PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
471 PetscCall(KSP_MatMult(ksp, Amat, P, W)); /* w <- Ap */
478 KSPCheckDot(ksp, beta);
481 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite matrix");
482 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
483 PetscCall(PetscInfo(ksp, "diverging due to indefinite or negative definite matrix\n"));
490 if (ksp->normtype == KSP_NORM_PRECONDITIONED && ksp->chknorm < i + 2) {
491 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
492 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
494 KSPCheckNorm(ksp, dp);
495 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED && ksp->chknorm < i + 2) {
497 KSPCheckNorm(ksp, dp);
498 } else if (ksp->normtype == KSP_NORM_NATURAL) {
499 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
502 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
506 KSPCheckDot(ksp, beta);
511 ksp->rnorm = dp;
512 PetscCall(KSPLogResidualHistory(ksp, dp));
513 PetscCall(KSPMonitor(ksp, i + 1, dp));
514 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
515 if (ksp->reason) break;
517 if ((ksp->normtype != KSP_NORM_PRECONDITIONED && ksp->normtype != KSP_NORM_NATURAL) || ksp->chknorm >= i + 2) {
518 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
519 PetscCall(KSP_MatMult(ksp, Amat, Z, S));
521 if (ksp->normtype != KSP_NORM_NATURAL || ksp->chknorm >= i + 2) {
527 KSPCheckDot(ksp, beta); /* beta <- z'*r */
531 } while (i < ksp->max_it);
532 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
542 PetscErrorCode KSPDestroy_CG(KSP ksp)
544 KSP_CG *cg = (KSP_CG *)ksp->data;
548 PetscCall(KSPDestroyDefault(ksp));
549 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetObjectiveTarget_C", NULL));
550 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
551 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", NULL));
552 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGUseSingleReduction_C", NULL));
553 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
562 PetscErrorCode KSPView_CG(KSP ksp, PetscViewer viewer)
564 KSP_CG *cg = (KSP_CG *)ksp->data;
582 PetscErrorCode KSPSetFromOptions_CG(KSP ksp, PetscOptionItems PetscOptionsObject)
584 KSP_CG *cg = (KSP_CG *)ksp->data;
593 if (flg) PetscCall(KSPCGUseSingleReduction(ksp, cg->singlereduction));
603 PetscErrorCode KSPCGSetType_CG(KSP ksp, KSPCGType type)
605 KSP_CG *cg = (KSP_CG *)ksp->data;
619 static PetscErrorCode KSPCGUseSingleReduction_CG(KSP ksp, PetscBool flg)
621 KSP_CG *cg = (KSP_CG *)ksp->data;
626 ksp->ops->solve = KSPSolve_CG_SingleReduction;
627 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
629 ksp->ops->solve = KSPSolve_CG;
630 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_CG));
635 PETSC_INTERN PetscErrorCode KSPBuildResidual_CG(KSP ksp, Vec t, Vec v, Vec *V)
638 PetscCall(VecCopy(ksp->work[0], v));
691 PETSC_EXTERN PetscErrorCode KSPCreate_CG(KSP ksp)
699 ksp->data = (void *)cg;
701 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
702 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
703 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
704 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
710 ksp->ops->setup = KSPSetUp_CG;
711 ksp->ops->solve = KSPSolve_CG;
712 ksp->ops->destroy = KSPDestroy_CG;
713 ksp->ops->view = KSPView_CG;
714 ksp->ops->setfromoptions = KSPSetFromOptions_CG;
715 ksp->ops->buildsolution = KSPBuildSolutionDefault;
716 ksp->ops->buildresidual = KSPBuildResidual_CG;
723 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", KSPCGSetType_CG));
724 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGUseSingleReduction_C", KSPCGUseSingleReduction_CG));
725 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_CG));
726 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetObjectiveTarget_C", KSPCGSetObjectiveTarget_CG));
727 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_CG));