Actual source code: stsles.c
slepc-3.21.0 2024-03-30
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: {
23: PetscFunctionBegin;
26: if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
27: PetscTryTypeMethod(st,setdefaultksp);
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: /*
32: This is done by all ST types except PRECOND.
33: The default is an LU direct solver, or GMRES+Jacobi if matmode=shell.
34: */
35: PetscErrorCode STSetDefaultKSP_Default(ST st)
36: {
37: PC pc;
38: PCType pctype;
39: KSPType ksptype;
41: PetscFunctionBegin;
42: PetscCall(KSPGetPC(st->ksp,&pc));
43: PetscCall(KSPGetType(st->ksp,&ksptype));
44: PetscCall(PCGetType(pc,&pctype));
45: if (!pctype && !ksptype) {
46: if (st->Pmat || st->Psplit) {
47: PetscCall(KSPSetType(st->ksp,KSPBCGS));
48: PetscCall(PCSetType(pc,PCBJACOBI));
49: } else if (st->matmode == ST_MATMODE_SHELL) {
50: PetscCall(KSPSetType(st->ksp,KSPGMRES));
51: PetscCall(PCSetType(pc,PCJACOBI));
52: } else {
53: PetscCall(KSPSetType(st->ksp,KSPPREONLY));
54: PetscCall(PCSetType(pc,st->asymm?PCCHOLESKY:PCLU));
55: }
56: }
57: PetscCall(KSPSetErrorIfNotConverged(st->ksp,PETSC_TRUE));
58: PetscFunctionReturn(PETSC_SUCCESS);
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
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(), STMatMultHermitianTranspose()
78: @*/
79: PetscErrorCode STMatMult(ST st,PetscInt k,Vec x,Vec y)
80: {
81: PetscFunctionBegin;
86: STCheckMatrices(st,1);
87: PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
88: PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
89: PetscCall(VecSetErrorIfLocked(y,3));
91: if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
92: PetscCall(VecLockReadPush(x));
93: PetscCall(PetscLogEventBegin(ST_MatMult,st,x,y,0));
94: if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
95: else PetscCall(MatMult(st->T[k],x,y));
96: PetscCall(PetscLogEventEnd(ST_MatMult,st,x,y,0));
97: PetscCall(VecLockReadPop(x));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /*@
102: STMatMultTranspose - Computes the matrix-vector product y = T[k]^T x, where T[k] is
103: the k-th matrix of the spectral transformation.
105: Neighbor-wise Collective
107: Input Parameters:
108: + st - the spectral transformation context
109: . k - index of matrix to use
110: - x - the vector to be multiplied
112: Output Parameter:
113: . y - the result
115: Level: developer
117: .seealso: STMatMult(), STMatMultHermitianTranspose()
118: @*/
119: PetscErrorCode STMatMultTranspose(ST st,PetscInt k,Vec x,Vec y)
120: {
121: PetscFunctionBegin;
126: STCheckMatrices(st,1);
127: PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
128: PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
129: PetscCall(VecSetErrorIfLocked(y,3));
131: if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
132: PetscCall(VecLockReadPush(x));
133: PetscCall(PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0));
134: if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
135: else PetscCall(MatMultTranspose(st->T[k],x,y));
136: PetscCall(PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0));
137: PetscCall(VecLockReadPop(x));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*@
142: STMatMultHermitianTranspose - Computes the matrix-vector product y = T[k]^H x, where T[k] is
143: the k-th matrix of the spectral transformation.
145: Neighbor-wise Collective
147: Input Parameters:
148: + st - the spectral transformation context
149: . k - index of matrix to use
150: - x - the vector to be multiplied
152: Output Parameter:
153: . y - the result
155: Level: developer
157: .seealso: STMatMult(), STMatMultTranspose()
158: @*/
159: PetscErrorCode STMatMultHermitianTranspose(ST st,PetscInt k,Vec x,Vec y)
160: {
161: PetscFunctionBegin;
166: STCheckMatrices(st,1);
167: PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
168: PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
169: PetscCall(VecSetErrorIfLocked(y,3));
171: if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
172: PetscCall(VecLockReadPush(x));
173: PetscCall(PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0));
174: if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
175: else PetscCall(MatMultHermitianTranspose(st->T[k],x,y));
176: PetscCall(PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0));
177: PetscCall(VecLockReadPop(x));
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /*@
182: STMatSolve - Solves P x = b, where P is the preconditioner matrix of
183: the spectral transformation, using a KSP object stored internally.
185: Collective
187: Input Parameters:
188: + st - the spectral transformation context
189: - b - right hand side vector
191: Output Parameter:
192: . x - computed solution
194: Level: developer
196: .seealso: STMatSolveTranspose(), STMatSolveHermitianTranspose(), STMatMatSolve()
197: @*/
198: PetscErrorCode STMatSolve(ST st,Vec b,Vec x)
199: {
200: PetscFunctionBegin;
204: STCheckMatrices(st,1);
205: PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
206: PetscCall(VecSetErrorIfLocked(x,3));
208: if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
209: PetscCall(VecLockReadPush(b));
210: PetscCall(PetscLogEventBegin(ST_MatSolve,st,b,x,0));
211: if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
212: else PetscCall(KSPSolve(st->ksp,b,x));
213: PetscCall(PetscLogEventEnd(ST_MatSolve,st,b,x,0));
214: PetscCall(VecLockReadPop(b));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*@
219: STMatMatSolve - Solves P X = B, where P is the preconditioner matrix of
220: the spectral transformation, using a KSP object stored internally.
222: Collective
224: Input Parameters:
225: + st - the spectral transformation context
226: - B - right hand side vectors
228: Output Parameter:
229: . X - computed solutions
231: Level: developer
233: .seealso: STMatSolve()
234: @*/
235: PetscErrorCode STMatMatSolve(ST st,Mat B,Mat X)
236: {
237: PetscFunctionBegin;
241: STCheckMatrices(st,1);
243: if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
244: PetscCall(PetscLogEventBegin(ST_MatSolve,st,B,X,0));
245: if (!st->P) PetscCall(MatCopy(B,X,SAME_NONZERO_PATTERN)); /* P=NULL means identity matrix */
246: else PetscCall(KSPMatSolve(st->ksp,B,X));
247: PetscCall(PetscLogEventEnd(ST_MatSolve,st,B,X,0));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: /*@
252: STMatSolveTranspose - Solves P^T x = b, where P is the preconditioner matrix of
253: the spectral transformation, using a KSP object stored internally.
255: Collective
257: Input Parameters:
258: + st - the spectral transformation context
259: - b - right hand side vector
261: Output Parameter:
262: . x - computed solution
264: Level: developer
266: .seealso: STMatSolve()
267: @*/
268: PetscErrorCode STMatSolveTranspose(ST st,Vec b,Vec x)
269: {
270: PetscFunctionBegin;
274: STCheckMatrices(st,1);
275: PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
276: PetscCall(VecSetErrorIfLocked(x,3));
278: if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
279: PetscCall(VecLockReadPush(b));
280: PetscCall(PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0));
281: if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
282: else PetscCall(KSPSolveTranspose(st->ksp,b,x));
283: PetscCall(PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0));
284: PetscCall(VecLockReadPop(b));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*@
289: STMatSolveHermitianTranspose - Solves P^H x = b, where P is the preconditioner matrix of
290: the spectral transformation, using a KSP object stored internally.
292: Collective
294: Input Parameters:
295: + st - the spectral transformation context
296: - b - right hand side vector
298: Output Parameter:
299: . x - computed solution
301: Level: developer
303: .seealso: STMatSolve()
304: @*/
305: PetscErrorCode STMatSolveHermitianTranspose(ST st,Vec b,Vec x)
306: {
307: Vec w;
309: PetscFunctionBegin;
313: STCheckMatrices(st,1);
314: PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
315: PetscCall(VecSetErrorIfLocked(x,3));
317: if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
318: PetscCall(VecLockReadPush(b));
319: PetscCall(PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0));
320: if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
321: else {
322: PetscCall(VecDuplicate(b,&w));
323: PetscCall(VecCopy(b,w));
324: PetscCall(VecConjugate(w));
325: PetscCall(KSPSolveTranspose(st->ksp,w,x));
326: PetscCall(VecDestroy(&w));
327: PetscCall(VecConjugate(x));
328: }
329: PetscCall(PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0));
330: PetscCall(VecLockReadPop(b));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: PetscErrorCode STCheckFactorPackage(ST st)
335: {
336: PC pc;
337: PetscMPIInt size;
338: PetscBool flg;
339: MatSolverType stype;
341: PetscFunctionBegin;
342: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)st),&size));
343: if (size==1) PetscFunctionReturn(PETSC_SUCCESS);
344: PetscCall(KSPGetPC(st->ksp,&pc));
345: PetscCall(PCFactorGetMatSolverType(pc,&stype));
346: if (stype) { /* currently selected PC is a factorization */
347: PetscCall(PetscStrcmp(stype,MATSOLVERPETSC,&flg));
348: PetscCheck(!flg,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");
349: }
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: /*@
354: STSetKSP - Sets the KSP object associated with the spectral
355: transformation.
357: Collective
359: Input Parameters:
360: + st - the spectral transformation context
361: - ksp - the linear system context
363: Level: advanced
365: .seealso: STGetKSP()
366: @*/
367: PetscErrorCode STSetKSP(ST st,KSP ksp)
368: {
369: PetscFunctionBegin;
372: PetscCheckSameComm(st,1,ksp,2);
373: STCheckNotSeized(st,1);
374: PetscCall(PetscObjectReference((PetscObject)ksp));
375: PetscCall(KSPDestroy(&st->ksp));
376: st->ksp = ksp;
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*@
381: STGetKSP - Gets the KSP object associated with the spectral
382: transformation.
384: Collective
386: Input Parameter:
387: . st - the spectral transformation context
389: Output Parameter:
390: . ksp - the linear system context
392: Level: intermediate
394: .seealso: STSetKSP()
395: @*/
396: PetscErrorCode STGetKSP(ST st,KSP* ksp)
397: {
398: PetscFunctionBegin;
400: PetscAssertPointer(ksp,2);
401: if (!st->ksp) {
402: PetscCall(KSPCreate(PetscObjectComm((PetscObject)st),&st->ksp));
403: PetscCall(PetscObjectIncrementTabLevel((PetscObject)st->ksp,(PetscObject)st,1));
404: PetscCall(KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix));
405: PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
406: PetscCall(PetscObjectSetOptions((PetscObject)st->ksp,((PetscObject)st)->options));
407: PetscCall(KSPSetTolerances(st->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
408: }
409: *ksp = st->ksp;
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: PetscErrorCode STCheckNullSpace_Default(ST st,BV V)
414: {
415: PetscInt nc,i,c;
416: PetscReal norm;
417: Vec *T,w,vi;
418: Mat A;
419: PC pc;
420: MatNullSpace nullsp;
422: PetscFunctionBegin;
423: PetscCall(BVGetNumConstraints(V,&nc));
424: PetscCall(PetscMalloc1(nc,&T));
425: if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
426: PetscCall(KSPGetPC(st->ksp,&pc));
427: PetscCall(PCGetOperators(pc,&A,NULL));
428: PetscCall(MatCreateVecs(A,NULL,&w));
429: c = 0;
430: for (i=0;i<nc;i++) {
431: PetscCall(BVGetColumn(V,-nc+i,&vi));
432: PetscCall(MatMult(A,vi,w));
433: PetscCall(VecNorm(w,NORM_2,&norm));
434: if (norm < 10.0*PETSC_SQRT_MACHINE_EPSILON) {
435: PetscCall(PetscInfo(st,"Vector %" PetscInt_FMT " included in the nullspace of OP, norm=%g\n",i,(double)norm));
436: PetscCall(BVCreateVec(V,T+c));
437: PetscCall(VecCopy(vi,T[c]));
438: PetscCall(VecNormalize(T[c],NULL));
439: c++;
440: } else PetscCall(PetscInfo(st,"Vector %" PetscInt_FMT " discarded as possible nullspace of OP, norm=%g\n",i,(double)norm));
441: PetscCall(BVRestoreColumn(V,-nc+i,&vi));
442: }
443: PetscCall(VecDestroy(&w));
444: if (c>0) {
445: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)st),PETSC_FALSE,c,T,&nullsp));
446: PetscCall(MatSetNullSpace(A,nullsp));
447: PetscCall(MatNullSpaceDestroy(&nullsp));
448: PetscCall(VecDestroyVecs(c,&T));
449: } else PetscCall(PetscFree(T));
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: /*@
454: STCheckNullSpace - Tests if constraint vectors are nullspace vectors.
456: Collective
458: Input Parameters:
459: + st - the spectral transformation context
460: - V - basis vectors to be checked
462: Notes:
463: Given a basis vectors object, this function tests each of its constraint
464: vectors to be a nullspace vector of the coefficient matrix of the
465: associated KSP object. All these nullspace vectors are passed to the KSP
466: object.
468: This function allows handling singular pencils and solving some problems
469: in which the nullspace is important (see the users guide for details).
471: Level: developer
473: .seealso: EPSSetDeflationSpace()
474: @*/
475: PetscErrorCode STCheckNullSpace(ST st,BV V)
476: {
477: PetscInt nc;
479: PetscFunctionBegin;
483: PetscCheckSameComm(st,1,V,2);
484: PetscCheck(st->state,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must call STSetUp() first");
486: PetscCall(BVGetNumConstraints(V,&nc));
487: if (nc) PetscTryTypeMethod(st,checknullspace,V);
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }