Actual source code: stsles.c

slepc-3.22.1 2024-10-28
Report Typos and Errors
  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_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
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: }