Actual source code: stsles.c
slepc-3.18.1 2022-11-02
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: ST interface routines related to the KSP object associated with it
12: */
14: #include <slepc/private/stimpl.h>
16: /*
17: This is used to set a default type for the KSP and PC objects.
18: It is called at STSetFromOptions (before KSPSetFromOptions)
19: and also at STSetUp (in case STSetFromOptions was not called).
20: */
21: PetscErrorCode STSetDefaultKSP(ST st)
22: {
25: if (!st->ksp) STGetKSP(st,&st->ksp);
26: PetscTryTypeMethod(st,setdefaultksp);
27: return 0;
28: }
30: /*
31: This is done by all ST types except PRECOND.
32: The default is an LU direct solver, or GMRES+Jacobi if matmode=shell.
33: */
34: PetscErrorCode STSetDefaultKSP_Default(ST st)
35: {
36: PC pc;
37: PCType pctype;
38: KSPType ksptype;
40: KSPGetPC(st->ksp,&pc);
41: KSPGetType(st->ksp,&ksptype);
42: PCGetType(pc,&pctype);
43: if (!pctype && !ksptype) {
44: if (st->Pmat || st->Psplit) {
45: KSPSetType(st->ksp,KSPBCGS);
46: PCSetType(pc,PCBJACOBI);
47: } else if (st->matmode == ST_MATMODE_SHELL) {
48: KSPSetType(st->ksp,KSPGMRES);
49: PCSetType(pc,PCJACOBI);
50: } else {
51: KSPSetType(st->ksp,KSPPREONLY);
52: PCSetType(pc,st->asymm?PCCHOLESKY:PCLU);
53: }
54: }
55: KSPSetErrorIfNotConverged(st->ksp,PETSC_TRUE);
56: return 0;
57: }
59: /*@
60: STMatMult - Computes the matrix-vector product y = T[k] x, where T[k] is
61: the k-th matrix of the spectral transformation.
63: Neighbor-wise Collective on st
65: Input Parameters:
66: + st - the spectral transformation context
67: . k - index of matrix to use
68: - x - the vector to be multiplied
70: Output Parameter:
71: . y - the result
73: Level: developer
75: .seealso: STMatMultTranspose()
76: @*/
77: PetscErrorCode STMatMult(ST st,PetscInt k,Vec x,Vec y)
78: {
83: STCheckMatrices(st,1);
86: VecSetErrorIfLocked(y,3);
88: if (st->state!=ST_STATE_SETUP) STSetUp(st);
89: VecLockReadPush(x);
90: PetscLogEventBegin(ST_MatMult,st,x,y,0);
91: if (!st->T[k]) VecCopy(x,y); /* T[k]=NULL means identity matrix */
92: else MatMult(st->T[k],x,y);
93: PetscLogEventEnd(ST_MatMult,st,x,y,0);
94: VecLockReadPop(x);
95: return 0;
96: }
98: /*@
99: STMatMultTranspose - Computes the matrix-vector product y = T[k]' x, where T[k] is
100: the k-th matrix of the spectral transformation.
102: Neighbor-wise Collective on st
104: Input Parameters:
105: + st - the spectral transformation context
106: . k - index of matrix to use
107: - x - the vector to be multiplied
109: Output Parameter:
110: . y - the result
112: Level: developer
114: .seealso: STMatMult()
115: @*/
116: PetscErrorCode STMatMultTranspose(ST st,PetscInt k,Vec x,Vec y)
117: {
122: STCheckMatrices(st,1);
125: VecSetErrorIfLocked(y,3);
127: if (st->state!=ST_STATE_SETUP) STSetUp(st);
128: VecLockReadPush(x);
129: PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0);
130: if (!st->T[k]) VecCopy(x,y); /* T[k]=NULL means identity matrix */
131: else MatMultTranspose(st->T[k],x,y);
132: PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0);
133: VecLockReadPop(x);
134: return 0;
135: }
137: /*@
138: STMatSolve - Solves P x = b, where P is the preconditioner matrix of
139: the spectral transformation, using a KSP object stored internally.
141: Collective on st
143: Input Parameters:
144: + st - the spectral transformation context
145: - b - right hand side vector
147: Output Parameter:
148: . x - computed solution
150: Level: developer
152: .seealso: STMatSolveTranspose()
153: @*/
154: PetscErrorCode STMatSolve(ST st,Vec b,Vec x)
155: {
159: STCheckMatrices(st,1);
161: VecSetErrorIfLocked(x,3);
163: if (st->state!=ST_STATE_SETUP) STSetUp(st);
164: VecLockReadPush(b);
165: PetscLogEventBegin(ST_MatSolve,st,b,x,0);
166: if (!st->P) VecCopy(b,x); /* P=NULL means identity matrix */
167: else KSPSolve(st->ksp,b,x);
168: PetscLogEventEnd(ST_MatSolve,st,b,x,0);
169: VecLockReadPop(b);
170: return 0;
171: }
173: /*@
174: STMatMatSolve - Solves P X = B, where P is the preconditioner matrix of
175: the spectral transformation, using a KSP object stored internally.
177: Collective on st
179: Input Parameters:
180: + st - the spectral transformation context
181: - B - right hand side vectors
183: Output Parameter:
184: . X - computed solutions
186: Level: developer
188: .seealso: STMatSolve()
189: @*/
190: PetscErrorCode STMatMatSolve(ST st,Mat B,Mat X)
191: {
195: STCheckMatrices(st,1);
197: if (st->state!=ST_STATE_SETUP) STSetUp(st);
198: PetscLogEventBegin(ST_MatSolve,st,B,X,0);
199: if (!st->P) MatCopy(B,X,SAME_NONZERO_PATTERN); /* P=NULL means identity matrix */
200: else KSPMatSolve(st->ksp,B,X);
201: PetscLogEventEnd(ST_MatSolve,st,B,X,0);
202: return 0;
203: }
205: /*@
206: STMatSolveTranspose - Solves P' x = b, where P is the preconditioner matrix of
207: the spectral transformation, using a KSP object stored internally.
209: Collective on st
211: Input Parameters:
212: + st - the spectral transformation context
213: - b - right hand side vector
215: Output Parameter:
216: . x - computed solution
218: Level: developer
220: .seealso: STMatSolve()
221: @*/
222: PetscErrorCode STMatSolveTranspose(ST st,Vec b,Vec x)
223: {
227: STCheckMatrices(st,1);
229: VecSetErrorIfLocked(x,3);
231: if (st->state!=ST_STATE_SETUP) STSetUp(st);
232: VecLockReadPush(b);
233: PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0);
234: if (!st->P) VecCopy(b,x); /* P=NULL means identity matrix */
235: else KSPSolveTranspose(st->ksp,b,x);
236: PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0);
237: VecLockReadPop(b);
238: return 0;
239: }
241: PetscErrorCode STCheckFactorPackage(ST st)
242: {
243: PC pc;
244: PetscMPIInt size;
245: PetscBool flg;
246: MatSolverType stype;
248: MPI_Comm_size(PetscObjectComm((PetscObject)st),&size);
249: if (size==1) return 0;
250: KSPGetPC(st->ksp,&pc);
251: PCFactorGetMatSolverType(pc,&stype);
252: if (stype) { /* currently selected PC is a factorization */
253: PetscStrcmp(stype,MATSOLVERPETSC,&flg);
255: }
256: return 0;
257: }
259: /*@
260: STSetKSP - Sets the KSP object associated with the spectral
261: transformation.
263: Collective on st
265: Input Parameters:
266: + st - the spectral transformation context
267: - ksp - the linear system context
269: Level: advanced
271: .seealso: STGetKSP()
272: @*/
273: PetscErrorCode STSetKSP(ST st,KSP ksp)
274: {
278: STCheckNotSeized(st,1);
279: PetscObjectReference((PetscObject)ksp);
280: KSPDestroy(&st->ksp);
281: st->ksp = ksp;
282: return 0;
283: }
285: /*@
286: STGetKSP - Gets the KSP object associated with the spectral
287: transformation.
289: Not Collective
291: Input Parameter:
292: . st - the spectral transformation context
294: Output Parameter:
295: . ksp - the linear system context
297: Level: intermediate
299: .seealso: STSetKSP()
300: @*/
301: PetscErrorCode STGetKSP(ST st,KSP* ksp)
302: {
305: if (!st->ksp) {
306: KSPCreate(PetscObjectComm((PetscObject)st),&st->ksp);
307: PetscObjectIncrementTabLevel((PetscObject)st->ksp,(PetscObject)st,1);
308: KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix);
309: KSPAppendOptionsPrefix(st->ksp,"st_");
310: PetscObjectSetOptions((PetscObject)st->ksp,((PetscObject)st)->options);
311: KSPSetTolerances(st->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
312: }
313: *ksp = st->ksp;
314: return 0;
315: }
317: PetscErrorCode STCheckNullSpace_Default(ST st,BV V)
318: {
319: PetscInt nc,i,c;
320: PetscReal norm;
321: Vec *T,w,vi;
322: Mat A;
323: PC pc;
324: MatNullSpace nullsp;
326: BVGetNumConstraints(V,&nc);
327: PetscMalloc1(nc,&T);
328: if (!st->ksp) STGetKSP(st,&st->ksp);
329: KSPGetPC(st->ksp,&pc);
330: PCGetOperators(pc,&A,NULL);
331: MatCreateVecs(A,NULL,&w);
332: c = 0;
333: for (i=0;i<nc;i++) {
334: BVGetColumn(V,-nc+i,&vi);
335: MatMult(A,vi,w);
336: VecNorm(w,NORM_2,&norm);
337: if (norm < 10.0*PETSC_SQRT_MACHINE_EPSILON) {
338: PetscInfo(st,"Vector %" PetscInt_FMT " included in the nullspace of OP, norm=%g\n",i,(double)norm);
339: BVCreateVec(V,T+c);
340: VecCopy(vi,T[c]);
341: VecNormalize(T[c],NULL);
342: c++;
343: } else PetscInfo(st,"Vector %" PetscInt_FMT " discarded as possible nullspace of OP, norm=%g\n",i,(double)norm);
344: BVRestoreColumn(V,-nc+i,&vi);
345: }
346: VecDestroy(&w);
347: if (c>0) {
348: MatNullSpaceCreate(PetscObjectComm((PetscObject)st),PETSC_FALSE,c,T,&nullsp);
349: MatSetNullSpace(A,nullsp);
350: MatNullSpaceDestroy(&nullsp);
351: VecDestroyVecs(c,&T);
352: } else PetscFree(T);
353: return 0;
354: }
356: /*@
357: STCheckNullSpace - Tests if constraint vectors are nullspace vectors.
359: Collective on st
361: Input Parameters:
362: + st - the spectral transformation context
363: - V - basis vectors to be checked
365: Notes:
366: Given a basis vectors object, this function tests each of its constraint
367: vectors to be a nullspace vector of the coefficient matrix of the
368: associated KSP object. All these nullspace vectors are passed to the KSP
369: object.
371: This function allows handling singular pencils and solving some problems
372: in which the nullspace is important (see the users guide for details).
374: Level: developer
376: .seealso: EPSSetDeflationSpace()
377: @*/
378: PetscErrorCode STCheckNullSpace(ST st,BV V)
379: {
380: PetscInt nc;
388: BVGetNumConstraints(V,&nc);
389: if (nc) PetscTryTypeMethod(st,checknullspace,V);
390: return 0;
391: }