1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: {
28: if (!st->ksp) { STGetKSP(st,&st->ksp); }
29: if (st->ops->setdefaultksp) { (*st->ops->setdefaultksp)(st); }
30: return(0);
31: }
33: /*
34: This is done by all ST types except PRECOND.
35: The default is an LU direct solver, or GMRES+Jacobi if matmode=shell.
36: */
37: PetscErrorCode STSetDefaultKSP_Default(ST st) 38: {
40: PC pc;
41: PCType pctype;
42: KSPType ksptype;
45: KSPGetPC(st->ksp,&pc);
46: KSPGetType(st->ksp,&ksptype);
47: PCGetType(pc,&pctype);
48: if (!pctype && !ksptype) {
49: if (st->matmode == ST_MATMODE_SHELL) {
50: KSPSetType(st->ksp,KSPGMRES);
51: PCSetType(pc,PCJACOBI);
52: } else {
53: KSPSetType(st->ksp,KSPPREONLY);
54: PCSetType(pc,st->asymm?PCCHOLESKY:PCLU);
55: }
56: }
57: KSPSetErrorIfNotConverged(st->ksp,PETSC_TRUE);
58: return(0);
59: }
61: /*@
62: STMatMult - Computes the matrix-vector product y = T[k] x, where T[k] is
63: the k-th matrix of the spectral transformation.
65: Neighbor-wise Collective on st
67: Input Parameters:
68: + st - the spectral transformation context
69: . k - index of matrix to use
70: - x - the vector to be multiplied
72: Output Parameter:
73: . y - the result
75: Level: developer
77: .seealso: STMatMultTranspose()
78: @*/
79: PetscErrorCode STMatMult(ST st,PetscInt k,Vec x,Vec y) 80: {
88: STCheckMatrices(st,1);
89: if (k<0 || k>=PetscMax(2,st->nmat)) SETERRQ1(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",st->nmat);
90: if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
91: VecSetErrorIfLocked(y,3);
93: if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
94: VecLockReadPush(x);
95: PetscLogEventBegin(ST_MatMult,st,x,y,0);
96: if (!st->T[k]) {
97: /* T[k]=NULL means identity matrix */
98: VecCopy(x,y);
99: } else {
100: MatMult(st->T[k],x,y);
101: }
102: PetscLogEventEnd(ST_MatMult,st,x,y,0);
103: VecLockReadPop(x);
104: return(0);
105: }
107: /*@
108: STMatMultTranspose - Computes the matrix-vector product y = T[k]' x, where T[k] is
109: the k-th matrix of the spectral transformation.
111: Neighbor-wise Collective on st
113: Input Parameters:
114: + st - the spectral transformation context
115: . k - index of matrix to use
116: - x - the vector to be multiplied
118: Output Parameter:
119: . y - the result
121: Level: developer
123: .seealso: STMatMult()
124: @*/
125: PetscErrorCode STMatMultTranspose(ST st,PetscInt k,Vec x,Vec y)126: {
134: STCheckMatrices(st,1);
135: if (k<0 || k>=PetscMax(2,st->nmat)) SETERRQ1(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",st->nmat);
136: if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
137: VecSetErrorIfLocked(y,3);
139: if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
140: VecLockReadPush(x);
141: PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0);
142: if (!st->T[k]) {
143: /* T[k]=NULL means identity matrix */
144: VecCopy(x,y);
145: } else {
146: MatMultTranspose(st->T[k],x,y);
147: }
148: PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0);
149: VecLockReadPop(x);
150: return(0);
151: }
153: /*@
154: STMatSolve - Solves P x = b, where P is the preconditioner matrix of
155: the spectral transformation, using a KSP object stored internally.
157: Collective on st
159: Input Parameters:
160: + st - the spectral transformation context
161: - b - right hand side vector
163: Output Parameter:
164: . x - computed solution
166: Level: developer
168: .seealso: STMatSolveTranspose()
169: @*/
170: PetscErrorCode STMatSolve(ST st,Vec b,Vec x)171: {
178: STCheckMatrices(st,1);
179: if (x == b) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
180: VecSetErrorIfLocked(x,3);
182: if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
183: VecLockReadPush(b);
184: PetscLogEventBegin(ST_MatSolve,st,b,x,0);
185: if (!st->P) {
186: /* P=NULL means identity matrix */
187: VecCopy(b,x);
188: } else {
189: KSPSolve(st->ksp,b,x);
190: }
191: PetscLogEventEnd(ST_MatSolve,st,b,x,0);
192: VecLockReadPop(b);
193: return(0);
194: }
196: /*@
197: STMatSolveTranspose - Solves P' x = b, where P is the preconditioner matrix of
198: the spectral transformation, using a KSP object stored internally.
200: Collective on st
202: Input Parameters:
203: . st - the spectral transformation context
204: . b - right hand side vector
206: Output Parameter:
207: . x - computed solution
209: Level: developer
211: .seealso: STMatSolve()
212: @*/
213: PetscErrorCode STMatSolveTranspose(ST st,Vec b,Vec x)214: {
221: STCheckMatrices(st,1);
222: if (x == b) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
223: VecSetErrorIfLocked(x,3);
225: if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
226: VecLockReadPush(b);
227: PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0);
228: if (!st->P) {
229: /* P=NULL means identity matrix */
230: VecCopy(b,x);
231: } else {
232: KSPSolveTranspose(st->ksp,b,x);
233: }
234: PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0);
235: VecLockReadPop(b);
236: return(0);
237: }
239: PetscErrorCode STCheckFactorPackage(ST st)240: {
242: PC pc;
243: PetscMPIInt size;
244: PetscBool flg;
245: 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);
254: if (flg) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"You chose to solve linear systems with a factorization, but in parallel runs you need to select an external package; see the users guide for details");
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
270: @*/
271: PetscErrorCode STSetKSP(ST st,KSP ksp)272: {
279: STCheckNotSeized(st,1);
280: PetscObjectReference((PetscObject)ksp);
281: KSPDestroy(&st->ksp);
282: st->ksp = ksp;
283: PetscLogObjectParent((PetscObject)st,(PetscObject)st->ksp);
284: return(0);
285: }
287: /*@
288: STGetKSP - Gets the KSP object associated with the spectral
289: transformation.
291: Not Collective
293: Input Parameter:
294: . st - the spectral transformation context
296: Output Parameter:
297: . ksp - the linear system context
299: Level: intermediate
300: @*/
301: PetscErrorCode STGetKSP(ST st,KSP* ksp)302: {
308: if (!st->ksp) {
309: KSPCreate(PetscObjectComm((PetscObject)st),&st->ksp);
310: PetscObjectIncrementTabLevel((PetscObject)st->ksp,(PetscObject)st,1);
311: KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix);
312: KSPAppendOptionsPrefix(st->ksp,"st_");
313: PetscLogObjectParent((PetscObject)st,(PetscObject)st->ksp);
314: PetscObjectSetOptions((PetscObject)st->ksp,((PetscObject)st)->options);
315: KSPSetTolerances(st->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
316: }
317: *ksp = st->ksp;
318: return(0);
319: }
321: PetscErrorCode STCheckNullSpace_Default(ST st,BV V)322: {
324: PetscInt nc,i,c;
325: PetscReal norm;
326: Vec *T,w,vi;
327: Mat A;
328: PC pc;
329: MatNullSpace nullsp;
332: BVGetNumConstraints(V,&nc);
333: PetscMalloc1(nc,&T);
334: if (!st->ksp) { STGetKSP(st,&st->ksp); }
335: KSPGetPC(st->ksp,&pc);
336: PCGetOperators(pc,&A,NULL);
337: MatCreateVecs(A,NULL,&w);
338: c = 0;
339: for (i=0;i<nc;i++) {
340: BVGetColumn(V,-nc+i,&vi);
341: MatMult(A,vi,w);
342: VecNorm(w,NORM_2,&norm);
343: if (norm < 1e-8) {
344: PetscInfo2(st,"Vector %D norm=%g\n",i,(double)norm);
345: BVCreateVec(V,T+c);
346: VecCopy(vi,T[c]);
347: VecNormalize(T[c],NULL);
348: c++;
349: }
350: BVRestoreColumn(V,-nc+i,&vi);
351: }
352: VecDestroy(&w);
353: if (c>0) {
354: MatNullSpaceCreate(PetscObjectComm((PetscObject)st),PETSC_FALSE,c,T,&nullsp);
355: MatSetNullSpace(A,nullsp);
356: MatNullSpaceDestroy(&nullsp);
357: VecDestroyVecs(c,&T);
358: } else {
359: PetscFree(T);
360: }
361: return(0);
362: }
364: /*@
365: STCheckNullSpace - Given a basis vectors object, this function tests each
366: of its constraint vectors to be a nullspace vector of the coefficient
367: matrix of the associated KSP object. All these nullspace vectors are passed
368: to the KSP object.
370: Collective on st
372: Input Parameters:
373: + st - the spectral transformation context
374: - V - basis vectors to be checked
376: Note:
377: This function allows to handle singular pencils and to solve some problems
378: in which the nullspace is important (see the users guide for details).
380: Level: developer
382: .seealso: EPSSetDeflationSpace()
383: @*/
384: PetscErrorCode STCheckNullSpace(ST st,BV V)385: {
387: PetscInt nc;
394: if (!st->state) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must call STSolve() first");
396: BVGetNumConstraints(V,&nc);
397: if (nc && st->ops->checknullspace) {
398: (*st->ops->checknullspace)(st,V);
399: }
400: return(0);
401: }