Actual source code: stsles.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcst.h" I*/
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->shift_matrix == ST_MATMODE_SHELL) {
50: KSPSetType(st->ksp,KSPGMRES);
51: PCSetType(pc,PCJACOBI);
52: } else {
53: KSPSetType(st->ksp,KSPPREONLY);
54: PCSetType(pc,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: 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: 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: {
173: PetscInt its;
174: PetscBool flg;
180: STCheckMatrices(st,1);
181: if (x == b) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
182: VecSetErrorIfLocked(x,3);
184: if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
185: VecLockReadPush(b);
186: PetscLogEventBegin(ST_MatSolve,st,b,x,0);
187: PetscObjectTypeCompareAny((PetscObject)st,&flg,STPRECOND,STSHELL,"");
188: if (!flg && !st->P) {
189: /* P=NULL means identity matrix */
190: VecCopy(b,x);
191: } else {
192: if (!st->ksp) { STGetKSP(st,&st->ksp); }
193: KSPSolve(st->ksp,b,x);
194: KSPGetIterationNumber(st->ksp,&its);
195: PetscInfo1(st,"Linear solve iterations=%D\n",its);
196: }
197: PetscLogEventEnd(ST_MatSolve,st,b,x,0);
198: VecLockReadPop(b);
199: return(0);
200: }
202: /*@
203: STMatSolveTranspose - Solves P' x = b, where P is the preconditioner matrix of
204: the spectral transformation, using a KSP object stored internally.
206: Collective on ST
208: Input Parameters:
209: . st - the spectral transformation context
210: . b - right hand side vector
212: Output Parameter:
213: . x - computed solution
215: Level: developer
217: .seealso: STMatSolve()
218: @*/
219: PetscErrorCode STMatSolveTranspose(ST st,Vec b,Vec x)
220: {
222: PetscInt its;
223: PetscBool flg;
229: STCheckMatrices(st,1);
230: if (x == b) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
231: VecSetErrorIfLocked(x,3);
233: if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
234: VecLockReadPush(b);
235: PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0);
236: PetscObjectTypeCompareAny((PetscObject)st,&flg,STPRECOND,STSHELL,"");
237: if (!flg && !st->P) {
238: /* P=NULL means identity matrix */
239: VecCopy(b,x);
240: } else {
241: if (!st->ksp) { STGetKSP(st,&st->ksp); }
242: KSPSolveTranspose(st->ksp,b,x);
243: KSPGetIterationNumber(st->ksp,&its);
244: PetscInfo1(st,"Linear solve iterations=%D\n",its);
245: }
246: PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0);
247: VecLockReadPop(b);
248: return(0);
249: }
251: /*
252: STMatSetHermitian - Sets the Hermitian flag to the ST matrix.
254: Input Parameters:
255: . st - the spectral transformation context
256: . M - matrix
257: */
258: PetscErrorCode STMatSetHermitian(ST st,Mat M)
259: {
261: PetscBool set,aherm,mherm;
262: PetscInt i;
265: #if defined(PETSC_USE_COMPLEX)
266: if (PetscImaginaryPart(st->sigma)!=0.0) mherm = PETSC_FALSE;
267: else {
268: mherm = PETSC_TRUE;
269: for (i=0;i<st->nmat;i++) {
270: MatIsHermitianKnown(st->A[i],&set,&aherm);
271: if (!set) aherm = PETSC_FALSE;
272: if (!aherm) { mherm = PETSC_FALSE; break; }
273: if (PetscRealPart(st->sigma)==0.0) break;
274: }
275: }
276: MatSetOption(M,MAT_HERMITIAN,mherm);
277: #else
278: mherm = PETSC_TRUE;
279: for (i=0;i<st->nmat;i++) {
280: MatIsSymmetricKnown(st->A[i],&set,&aherm);
281: if (!set) aherm = PETSC_FALSE;
282: if (!aherm) { mherm = PETSC_FALSE; break; }
283: }
284: MatSetOption(M,MAT_SYMMETRIC,mherm);
285: #endif
286: return(0);
287: }
289: PetscErrorCode STCheckFactorPackage(ST st)
290: {
292: PC pc;
293: PetscMPIInt size;
294: PetscBool flg;
295: MatSolverType stype;
298: MPI_Comm_size(PetscObjectComm((PetscObject)st),&size);
299: if (size==1) return(0);
300: KSPGetPC(st->ksp,&pc);
301: PCFactorGetMatSolverType(pc,&stype);
302: if (stype) { /* currently selected PC is a factorization */
303: PetscStrcmp(stype,MATSOLVERPETSC,&flg);
304: 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");
305: }
306: return(0);
307: }
309: /*@
310: STSetKSP - Sets the KSP object associated with the spectral
311: transformation.
313: Collective on ST
315: Input Parameters:
316: + st - the spectral transformation context
317: - ksp - the linear system context
319: Level: advanced
320: @*/
321: PetscErrorCode STSetKSP(ST st,KSP ksp)
322: {
329: PetscObjectReference((PetscObject)ksp);
330: KSPDestroy(&st->ksp);
331: st->ksp = ksp;
332: PetscLogObjectParent((PetscObject)st,(PetscObject)st->ksp);
333: return(0);
334: }
336: /*@
337: STGetKSP - Gets the KSP object associated with the spectral
338: transformation.
340: Not Collective
342: Input Parameter:
343: . st - the spectral transformation context
345: Output Parameter:
346: . ksp - the linear system context
348: Level: intermediate
349: @*/
350: PetscErrorCode STGetKSP(ST st,KSP* ksp)
351: {
357: if (!st->ksp) {
358: KSPCreate(PetscObjectComm((PetscObject)st),&st->ksp);
359: PetscObjectIncrementTabLevel((PetscObject)st->ksp,(PetscObject)st,1);
360: KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix);
361: KSPAppendOptionsPrefix(st->ksp,"st_");
362: PetscLogObjectParent((PetscObject)st,(PetscObject)st->ksp);
363: PetscObjectSetOptions((PetscObject)st->ksp,((PetscObject)st)->options);
364: KSPSetTolerances(st->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
365: }
366: *ksp = st->ksp;
367: return(0);
368: }
370: PetscErrorCode STCheckNullSpace_Default(ST st,BV V)
371: {
373: PetscInt nc,i,c;
374: PetscReal norm;
375: Vec *T,w,vi;
376: Mat A;
377: PC pc;
378: MatNullSpace nullsp;
381: BVGetNumConstraints(V,&nc);
382: PetscMalloc1(nc,&T);
383: if (!st->ksp) { STGetKSP(st,&st->ksp); }
384: KSPGetPC(st->ksp,&pc);
385: PCGetOperators(pc,&A,NULL);
386: MatCreateVecs(A,NULL,&w);
387: c = 0;
388: for (i=0;i<nc;i++) {
389: BVGetColumn(V,-nc+i,&vi);
390: MatMult(A,vi,w);
391: VecNorm(w,NORM_2,&norm);
392: if (norm < 1e-8) {
393: PetscInfo2(st,"Vector %D norm=%g\n",i,(double)norm);
394: BVCreateVec(V,T+c);
395: VecCopy(vi,T[c]);
396: c++;
397: }
398: BVRestoreColumn(V,-nc+i,&vi);
399: }
400: VecDestroy(&w);
401: if (c>0) {
402: MatNullSpaceCreate(PetscObjectComm((PetscObject)st),PETSC_FALSE,c,T,&nullsp);
403: MatSetNullSpace(A,nullsp);
404: MatNullSpaceDestroy(&nullsp);
405: VecDestroyVecs(c,&T);
406: } else {
407: PetscFree(T);
408: }
409: return(0);
410: }
412: /*@
413: STCheckNullSpace - Given a basis vectors object, this function tests each
414: of its constraint vectors to be a nullspace vector of the coefficient
415: matrix of the associated KSP object. All these nullspace vectors are passed
416: to the KSP object.
418: Collective on ST
420: Input Parameters:
421: + st - the spectral transformation context
422: - V - basis vectors to be checked
424: Note:
425: This function allows to handle singular pencils and to solve some problems
426: in which the nullspace is important (see the users guide for details).
428: Level: developer
430: .seealso: EPSSetDeflationSpace()
431: @*/
432: PetscErrorCode STCheckNullSpace(ST st,BV V)
433: {
435: PetscInt nc;
442: if (!st->state) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must call STSolve() first");
444: BVGetNumConstraints(V,&nc);
445: if (nc && st->ops->checknullspace) {
446: (*st->ops->checknullspace)(st,V);
447: }
448: return(0);
449: }