Actual source code: precond.c
slepc-main 2024-11-09
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: Implements the ST class for preconditioned eigenvalue methods
12: */
14: #include <slepc/private/stimpl.h>
16: typedef struct {
17: PetscBool ksphasmat; /* the KSP must have the same matrix as PC */
18: } ST_PRECOND;
20: static PetscErrorCode STSetDefaultKSP_Precond(ST st)
21: {
22: PC pc;
23: PCType pctype;
24: PetscBool t0,t1;
26: PetscFunctionBegin;
27: PetscCall(KSPGetPC(st->ksp,&pc));
28: PetscCall(PCGetType(pc,&pctype));
29: if (!pctype && st->A && st->A[0]) {
30: if (st->matmode == ST_MATMODE_SHELL) PetscCall(PCSetType(pc,PCJACOBI));
31: else {
32: PetscCall(MatHasOperation(st->A[0],MATOP_DUPLICATE,&t0));
33: if (st->nmat>1) PetscCall(MatHasOperation(st->A[0],MATOP_AXPY,&t1));
34: else t1 = PETSC_TRUE;
35: PetscCall(PCSetType(pc,(t0 && t1)?PCBJACOBI:PCNONE));
36: }
37: }
38: PetscCall(KSPSetErrorIfNotConverged(st->ksp,PETSC_FALSE));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static PetscErrorCode STPostSolve_Precond(ST st)
43: {
44: PetscFunctionBegin;
45: if (st->matmode == ST_MATMODE_INPLACE && !(st->Pmat || (PetscAbsScalar(st->sigma)>=PETSC_MAX_REAL && st->nmat>1))) {
46: if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
47: else PetscCall(MatShift(st->A[0],st->sigma));
48: st->Astate[0] = ((PetscObject)st->A[0])->state;
49: st->state = ST_STATE_INITIAL;
50: st->opready = PETSC_FALSE;
51: }
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: /*
56: Operator (precond):
57: Op P M
58: if nmat=1: --- A-sI NULL
59: if nmat=2: --- A-sB NULL
60: */
61: static PetscErrorCode STComputeOperator_Precond(ST st)
62: {
63: PetscFunctionBegin;
64: /* if the user did not set the shift, use the target value */
65: if (!st->sigma_set) st->sigma = st->defsigma;
66: st->M = NULL;
68: /* build custom preconditioner from the split matrices */
69: if (st->Psplit) {
70: if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
71: PetscCall(PetscObjectReference((PetscObject)st->Psplit[0]));
72: PetscCall(MatDestroy(&st->Pmat));
73: st->Pmat = st->Psplit[0];
74: } else if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL) PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
75: }
77: /* P = A-sigma*B */
78: if (st->Pmat) {
79: PetscCall(PetscObjectReference((PetscObject)st->Pmat));
80: PetscCall(MatDestroy(&st->P));
81: st->P = st->Pmat;
82: } else {
83: PetscCall(PetscObjectReference((PetscObject)st->A[1]));
84: PetscCall(MatDestroy(&st->T[0]));
85: st->T[0] = st->A[1];
86: if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
87: PetscCall(PetscObjectReference((PetscObject)st->T[0]));
88: PetscCall(MatDestroy(&st->P));
89: st->P = st->T[0];
90: } else if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL) {
91: PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[1]));
92: PetscCall(PetscObjectReference((PetscObject)st->T[1]));
93: PetscCall(MatDestroy(&st->P));
94: st->P = st->T[1];
95: }
96: }
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode STSetUp_Precond(ST st)
101: {
102: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
104: PetscFunctionBegin;
105: if (st->P) {
106: PetscCall(ST_KSPSetOperators(st,ctx->ksphasmat?st->P:NULL,st->P));
107: /* NOTE: we do not call KSPSetUp() here because some eigensolvers such as JD require a lazy setup */
108: }
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: static PetscErrorCode STSetShift_Precond(ST st,PetscScalar newshift)
113: {
114: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
116: PetscFunctionBegin;
117: if (st->Psplit) { /* update custom preconditioner from the split matrices */
118: if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL || st->nmat==1) PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PETSC_FALSE,PETSC_TRUE,&st->Pmat));
119: }
120: if (st->transform && !st->Pmat) {
121: PetscCall(STMatMAXPY_Private(st,-newshift,-st->sigma,0,NULL,PETSC_FALSE,PETSC_FALSE,&st->T[1]));
122: if (st->P!=st->T[1]) {
123: PetscCall(PetscObjectReference((PetscObject)st->T[1]));
124: PetscCall(MatDestroy(&st->P));
125: st->P = st->T[1];
126: }
127: }
128: if (st->P) PetscCall(ST_KSPSetOperators(st,ctx->ksphasmat?st->P:NULL,st->P));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: static PetscErrorCode STPrecondSetKSPHasMat_Precond(ST st,PetscBool ksphasmat)
133: {
134: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
136: PetscFunctionBegin;
137: if (ctx->ksphasmat != ksphasmat) {
138: ctx->ksphasmat = ksphasmat;
139: st->state = ST_STATE_INITIAL;
140: }
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /*@
145: STPrecondSetKSPHasMat - Sets a flag indicating that during STSetUp the coefficient
146: matrix of the KSP linear solver (A) must be set to be the same matrix as the
147: preconditioner (P).
149: Collective
151: Input Parameters:
152: + st - the spectral transformation context
153: - ksphasmat - the flag
155: Notes:
156: Often, the preconditioner matrix is used only in the PC object, but
157: in some solvers this matrix must be provided also as the A-matrix in
158: the KSP object.
160: Level: developer
162: .seealso: STPrecondGetKSPHasMat(), STSetShift()
163: @*/
164: PetscErrorCode STPrecondSetKSPHasMat(ST st,PetscBool ksphasmat)
165: {
166: PetscFunctionBegin;
169: PetscTryMethod(st,"STPrecondSetKSPHasMat_C",(ST,PetscBool),(st,ksphasmat));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode STPrecondGetKSPHasMat_Precond(ST st,PetscBool *ksphasmat)
174: {
175: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
177: PetscFunctionBegin;
178: *ksphasmat = ctx->ksphasmat;
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: /*@
183: STPrecondGetKSPHasMat - Returns the flag indicating if the coefficient
184: matrix of the KSP linear system (A) is set to be the same matrix as the
185: preconditioner (P).
187: Not Collective
189: Input Parameter:
190: . st - the spectral transformation context
192: Output Parameter:
193: . ksphasmat - the flag
195: Level: developer
197: .seealso: STPrecondSetKSPHasMat(), STSetShift()
198: @*/
199: PetscErrorCode STPrecondGetKSPHasMat(ST st,PetscBool *ksphasmat)
200: {
201: PetscFunctionBegin;
203: PetscAssertPointer(ksphasmat,2);
204: PetscUseMethod(st,"STPrecondGetKSPHasMat_C",(ST,PetscBool*),(st,ksphasmat));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: static PetscErrorCode STDestroy_Precond(ST st)
209: {
210: PetscFunctionBegin;
211: PetscCall(PetscFree(st->data));
212: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",NULL));
213: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",NULL));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: SLEPC_EXTERN PetscErrorCode STCreate_Precond(ST st)
218: {
219: ST_PRECOND *ctx;
221: PetscFunctionBegin;
222: PetscCall(PetscNew(&ctx));
223: st->data = (void*)ctx;
225: st->usesksp = PETSC_TRUE;
227: st->ops->apply = STApply_Generic;
228: st->ops->applymat = STApplyMat_Generic;
229: st->ops->applytrans = STApplyTranspose_Generic;
230: st->ops->applyhermtrans = STApplyHermitianTranspose_Generic;
231: st->ops->setshift = STSetShift_Precond;
232: st->ops->getbilinearform = STGetBilinearForm_Default;
233: st->ops->setup = STSetUp_Precond;
234: st->ops->computeoperator = STComputeOperator_Precond;
235: st->ops->postsolve = STPostSolve_Precond;
236: st->ops->destroy = STDestroy_Precond;
237: st->ops->setdefaultksp = STSetDefaultKSP_Precond;
239: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",STPrecondGetKSPHasMat_Precond));
240: PetscCall(PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",STPrecondSetKSPHasMat_Precond));
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }