Actual source code: stfunc.c

slepc-3.21.1 2024-04-26
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:    The ST interface routines, callable by users
 12: */

 14: #include <slepc/private/stimpl.h>

 16: PetscClassId     ST_CLASSID = 0;
 17: PetscLogEvent    ST_SetUp = 0,ST_ComputeOperator = 0,ST_Apply = 0,ST_ApplyTranspose = 0,ST_ApplyHermitianTranspose = 0,ST_MatSetUp = 0,ST_MatMult = 0,ST_MatMultTranspose = 0,ST_MatSolve = 0,ST_MatSolveTranspose = 0;
 18: static PetscBool STPackageInitialized = PETSC_FALSE;

 20: const char *STMatModes[] = {"COPY","INPLACE","SHELL","STMatMode","ST_MATMODE_",NULL};

 22: /*@C
 23:    STFinalizePackage - This function destroys everything in the Slepc interface
 24:    to the ST package. It is called from SlepcFinalize().

 26:    Level: developer

 28: .seealso: SlepcFinalize()
 29: @*/
 30: PetscErrorCode STFinalizePackage(void)
 31: {
 32:   PetscFunctionBegin;
 33:   PetscCall(PetscFunctionListDestroy(&STList));
 34:   STPackageInitialized = PETSC_FALSE;
 35:   STRegisterAllCalled  = PETSC_FALSE;
 36:   PetscFunctionReturn(PETSC_SUCCESS);
 37: }

 39: /*@C
 40:    STInitializePackage - This function initializes everything in the ST package.
 41:    It is called from PetscDLLibraryRegister() when using dynamic libraries, and
 42:    on the first call to STCreate() when using static libraries.

 44:    Level: developer

 46: .seealso: SlepcInitialize()
 47: @*/
 48: PetscErrorCode STInitializePackage(void)
 49: {
 50:   char           logList[256];
 51:   PetscBool      opt,pkg;
 52:   PetscClassId   classids[1];

 54:   PetscFunctionBegin;
 55:   if (STPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
 56:   STPackageInitialized = PETSC_TRUE;
 57:   /* Register Classes */
 58:   PetscCall(PetscClassIdRegister("Spectral Transform",&ST_CLASSID));
 59:   /* Register Constructors */
 60:   PetscCall(STRegisterAll());
 61:   /* Register Events */
 62:   PetscCall(PetscLogEventRegister("STSetUp",ST_CLASSID,&ST_SetUp));
 63:   PetscCall(PetscLogEventRegister("STComputeOperatr",ST_CLASSID,&ST_ComputeOperator));
 64:   PetscCall(PetscLogEventRegister("STApply",ST_CLASSID,&ST_Apply));
 65:   PetscCall(PetscLogEventRegister("STApplyTranspose",ST_CLASSID,&ST_ApplyTranspose));
 66:   PetscCall(PetscLogEventRegister("STApplyHermTrans",ST_CLASSID,&ST_ApplyHermitianTranspose));
 67:   PetscCall(PetscLogEventRegister("STMatSetUp",ST_CLASSID,&ST_MatSetUp));
 68:   PetscCall(PetscLogEventRegister("STMatMult",ST_CLASSID,&ST_MatMult));
 69:   PetscCall(PetscLogEventRegister("STMatMultTranspose",ST_CLASSID,&ST_MatMultTranspose));
 70:   PetscCall(PetscLogEventRegister("STMatSolve",ST_CLASSID,&ST_MatSolve));
 71:   PetscCall(PetscLogEventRegister("STMatSolveTranspose",ST_CLASSID,&ST_MatSolveTranspose));
 72:   /* Process Info */
 73:   classids[0] = ST_CLASSID;
 74:   PetscCall(PetscInfoProcessClass("st",1,&classids[0]));
 75:   /* Process summary exclusions */
 76:   PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
 77:   if (opt) {
 78:     PetscCall(PetscStrInList("st",logList,',',&pkg));
 79:     if (pkg) PetscCall(PetscLogEventDeactivateClass(ST_CLASSID));
 80:   }
 81:   /* Register package finalizer */
 82:   PetscCall(PetscRegisterFinalize(STFinalizePackage));
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }

 86: /*@
 87:    STReset - Resets the ST context to the initial state (prior to setup)
 88:    and destroys any allocated Vecs and Mats.

 90:    Collective

 92:    Input Parameter:
 93: .  st - the spectral transformation context

 95:    Level: advanced

 97: .seealso: STDestroy()
 98: @*/
 99: PetscErrorCode STReset(ST st)
100: {
101:   PetscFunctionBegin;
103:   if (!st) PetscFunctionReturn(PETSC_SUCCESS);
104:   STCheckNotSeized(st,1);
105:   PetscTryTypeMethod(st,reset);
106:   if (st->ksp) PetscCall(KSPReset(st->ksp));
107:   PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->T));
108:   PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->A));
109:   st->nmat = 0;
110:   PetscCall(PetscFree(st->Astate));
111:   PetscCall(MatDestroy(&st->Op));
112:   PetscCall(MatDestroy(&st->P));
113:   PetscCall(MatDestroy(&st->Pmat));
114:   PetscCall(MatDestroyMatrices(st->nsplit,&st->Psplit));
115:   st->nsplit = 0;
116:   PetscCall(VecDestroyVecs(st->nwork,&st->work));
117:   st->nwork = 0;
118:   PetscCall(VecDestroy(&st->wb));
119:   PetscCall(VecDestroy(&st->wht));
120:   PetscCall(VecDestroy(&st->D));
121:   st->state   = ST_STATE_INITIAL;
122:   st->opready = PETSC_FALSE;
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: /*@C
127:    STDestroy - Destroys ST context that was created with STCreate().

129:    Collective

131:    Input Parameter:
132: .  st - the spectral transformation context

134:    Level: beginner

136: .seealso: STCreate(), STSetUp()
137: @*/
138: PetscErrorCode STDestroy(ST *st)
139: {
140:   PetscFunctionBegin;
141:   if (!*st) PetscFunctionReturn(PETSC_SUCCESS);
143:   if (--((PetscObject)*st)->refct > 0) { *st = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
144:   PetscCall(STReset(*st));
145:   PetscTryTypeMethod(*st,destroy);
146:   PetscCall(KSPDestroy(&(*st)->ksp));
147:   PetscCall(PetscHeaderDestroy(st));
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

151: /*@
152:    STCreate - Creates a spectral transformation context.

154:    Collective

156:    Input Parameter:
157: .  comm - MPI communicator

159:    Output Parameter:
160: .  newst - location to put the spectral transformation context

162:    Level: beginner

164: .seealso: STSetUp(), STApply(), STDestroy(), ST
165: @*/
166: PetscErrorCode STCreate(MPI_Comm comm,ST *newst)
167: {
168:   ST             st;

170:   PetscFunctionBegin;
171:   PetscAssertPointer(newst,2);
172:   *newst = NULL;
173:   PetscCall(STInitializePackage());
174:   PetscCall(SlepcHeaderCreate(st,ST_CLASSID,"ST","Spectral Transformation","ST",comm,STDestroy,STView));

176:   st->A            = NULL;
177:   st->nmat         = 0;
178:   st->sigma        = 0.0;
179:   st->defsigma     = 0.0;
180:   st->matmode      = ST_MATMODE_COPY;
181:   st->str          = UNKNOWN_NONZERO_PATTERN;
182:   st->transform    = PETSC_FALSE;
183:   st->D            = NULL;
184:   st->Pmat         = NULL;
185:   st->Pmat_set     = PETSC_FALSE;
186:   st->Psplit       = NULL;
187:   st->nsplit       = 0;
188:   st->strp         = UNKNOWN_NONZERO_PATTERN;

190:   st->ksp          = NULL;
191:   st->usesksp      = PETSC_FALSE;
192:   st->nwork        = 0;
193:   st->work         = NULL;
194:   st->wb           = NULL;
195:   st->wht          = NULL;
196:   st->state        = ST_STATE_INITIAL;
197:   st->Astate       = NULL;
198:   st->T            = NULL;
199:   st->Op           = NULL;
200:   st->opseized     = PETSC_FALSE;
201:   st->opready      = PETSC_FALSE;
202:   st->P            = NULL;
203:   st->M            = NULL;
204:   st->sigma_set    = PETSC_FALSE;
205:   st->asymm        = PETSC_FALSE;
206:   st->aherm        = PETSC_FALSE;
207:   st->data         = NULL;

209:   *newst = st;
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }

213: /*
214:    Checks whether the ST matrices are all symmetric or hermitian.
215: */
216: static inline PetscErrorCode STMatIsSymmetricKnown(ST st,PetscBool *symm,PetscBool *herm)
217: {
218:   PetscInt       i;
219:   PetscBool      sbaij=PETSC_FALSE,set,flg=PETSC_FALSE;

221:   PetscFunctionBegin;
222:   /* check if problem matrices are all sbaij */
223:   for (i=0;i<st->nmat;i++) {
224:     PetscCall(PetscObjectTypeCompareAny((PetscObject)st->A[i],&sbaij,MATSEQSBAIJ,MATMPISBAIJ,""));
225:     if (!sbaij) break;
226:   }
227:   /* check if user has set the symmetric flag */
228:   *symm = PETSC_TRUE;
229:   for (i=0;i<st->nmat;i++) {
230:     PetscCall(MatIsSymmetricKnown(st->A[i],&set,&flg));
231:     if (!set || !flg) { *symm = PETSC_FALSE; break; }
232:   }
233:   if (sbaij) *symm = PETSC_TRUE;
234: #if defined(PETSC_USE_COMPLEX)
235:   /* check if user has set the hermitian flag */
236:   *herm = PETSC_TRUE;
237:   for (i=0;i<st->nmat;i++) {
238:     PetscCall(MatIsHermitianKnown(st->A[i],&set,&flg));
239:     if (!set || !flg) { *herm = PETSC_FALSE; break; }
240:   }
241: #else
242:   *herm = *symm;
243: #endif
244:   PetscFunctionReturn(PETSC_SUCCESS);
245: }

247: /*@
248:    STSetMatrices - Sets the matrices associated with the eigenvalue problem.

250:    Collective

252:    Input Parameters:
253: +  st - the spectral transformation context
254: .  n  - number of matrices in array A
255: -  A  - the array of matrices associated with the eigensystem

257:    Notes:
258:    It must be called before STSetUp(). If it is called again after STSetUp() then
259:    the ST object is reset.

261:    Level: intermediate

263: .seealso: STGetMatrix(), STGetNumMatrices(), STSetUp(), STReset()
264: @*/
265: PetscErrorCode STSetMatrices(ST st,PetscInt n,Mat A[])
266: {
267:   PetscInt       i;
268:   PetscBool      same=PETSC_TRUE;

270:   PetscFunctionBegin;
273:   PetscCheck(n>0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more matrices, you have %" PetscInt_FMT,n);
274:   PetscAssertPointer(A,3);
275:   PetscCheckSameComm(st,1,*A,3);
276:   STCheckNotSeized(st,1);
277:   PetscCheck(!st->nsplit || st->nsplit==n,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"The number of matrices must be the same as in STSetSplitPreconditioner()");

279:   if (st->state) {
280:     if (n!=st->nmat) same = PETSC_FALSE;
281:     for (i=0;same&&i<n;i++) {
282:       if (A[i]!=st->A[i]) same = PETSC_FALSE;
283:     }
284:     if (!same) PetscCall(STReset(st));
285:   } else same = PETSC_FALSE;
286:   if (!same) {
287:     PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->A));
288:     PetscCall(PetscCalloc1(PetscMax(2,n),&st->A));
289:     PetscCall(PetscFree(st->Astate));
290:     PetscCall(PetscMalloc1(PetscMax(2,n),&st->Astate));
291:   }
292:   for (i=0;i<n;i++) {
294:     PetscCall(PetscObjectReference((PetscObject)A[i]));
295:     PetscCall(MatDestroy(&st->A[i]));
296:     st->A[i] = A[i];
297:     st->Astate[i] = ((PetscObject)A[i])->state;
298:   }
299:   if (n==1) {
300:     st->A[1] = NULL;
301:     st->Astate[1] = 0;
302:   }
303:   st->nmat = n;
304:   if (same) st->state = ST_STATE_UPDATED;
305:   else st->state = ST_STATE_INITIAL;
306:   PetscCheck(!same || !st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Support for changing the matrices while using a split preconditioner is not implemented yet");
307:   st->opready = PETSC_FALSE;
308:   if (!same) PetscCall(STMatIsSymmetricKnown(st,&st->asymm,&st->aherm));
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

312: /*@
313:    STGetMatrix - Gets the matrices associated with the original eigensystem.

315:    Not Collective

317:    Input Parameters:
318: +  st - the spectral transformation context
319: -  k  - the index of the requested matrix (starting in 0)

321:    Output Parameters:
322: .  A - the requested matrix

324:    Level: intermediate

326: .seealso: STSetMatrices(), STGetNumMatrices()
327: @*/
328: PetscErrorCode STGetMatrix(ST st,PetscInt k,Mat *A)
329: {
330:   PetscFunctionBegin;
333:   PetscAssertPointer(A,3);
334:   STCheckMatrices(st,1);
335:   PetscCheck(k>=0 && k<st->nmat,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat-1);
336:   PetscCheck(((PetscObject)st->A[k])->state==st->Astate[k],PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot retrieve original matrices (have been modified)");
337:   *A = st->A[k];
338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

341: /*@
342:    STGetMatrixTransformed - Gets the matrices associated with the transformed eigensystem.

344:    Not Collective

346:    Input Parameters:
347: +  st - the spectral transformation context
348: -  k  - the index of the requested matrix (starting in 0)

350:    Output Parameters:
351: .  T - the requested matrix

353:    Level: developer

355: .seealso: STGetMatrix(), STGetNumMatrices()
356: @*/
357: PetscErrorCode STGetMatrixTransformed(ST st,PetscInt k,Mat *T)
358: {
359:   PetscFunctionBegin;
362:   PetscAssertPointer(T,3);
363:   STCheckMatrices(st,1);
364:   PetscCheck(k>=0 && k<st->nmat,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat-1);
365:   PetscCheck(st->T,PetscObjectComm((PetscObject)st),PETSC_ERR_POINTER,"There are no transformed matrices");
366:   *T = st->T[k];
367:   PetscFunctionReturn(PETSC_SUCCESS);
368: }

370: /*@
371:    STGetNumMatrices - Returns the number of matrices stored in the ST.

373:    Not Collective

375:    Input Parameter:
376: .  st - the spectral transformation context

378:    Output Parameters:
379: .  n - the number of matrices passed in STSetMatrices()

381:    Level: intermediate

383: .seealso: STSetMatrices()
384: @*/
385: PetscErrorCode STGetNumMatrices(ST st,PetscInt *n)
386: {
387:   PetscFunctionBegin;
389:   PetscAssertPointer(n,2);
390:   *n = st->nmat;
391:   PetscFunctionReturn(PETSC_SUCCESS);
392: }

394: /*@
395:    STResetMatrixState - Resets the stored state of the matrices in the ST.

397:    Logically Collective

399:    Input Parameter:
400: .  st - the spectral transformation context

402:    Note:
403:    This is useful in solvers where the user matrices are modified during
404:    the computation, as in nonlinear inverse iteration. The effect is that
405:    STGetMatrix() will retrieve the modified matrices as if they were
406:    the matrices originally provided by the user.

408:    Level: developer

410: .seealso: STGetMatrix(), EPSPowerSetNonlinear()
411: @*/
412: PetscErrorCode STResetMatrixState(ST st)
413: {
414:   PetscInt i;

416:   PetscFunctionBegin;
418:   for (i=0;i<st->nmat;i++) st->Astate[i] = ((PetscObject)st->A[i])->state;
419:   PetscFunctionReturn(PETSC_SUCCESS);
420: }

422: /*@
423:    STSetPreconditionerMat - Sets the matrix to be used to build the preconditioner.

425:    Collective

427:    Input Parameters:
428: +  st  - the spectral transformation context
429: -  mat - the matrix that will be used in constructing the preconditioner

431:    Notes:
432:    This matrix will be passed to the internal KSP object (via the last argument
433:    of KSPSetOperators()) as the matrix to be used when constructing the preconditioner.
434:    If no matrix is set or mat is set to NULL, A-sigma*B will be used
435:    to build the preconditioner, being sigma the value set by STSetShift().

437:    More precisely, this is relevant for spectral transformations that represent
438:    a rational matrix function, and use a KSP object for the denominator, called
439:    K in the description of STGetOperator(). It includes also the STPRECOND case.
440:    If the user has a good approximation to matrix K that can be used to build a
441:    cheap preconditioner, it can be passed with this function. Note that it affects
442:    only the Pmat argument of KSPSetOperators(), not the Amat argument.

444:    If a preconditioner matrix is set, the default is to use an iterative KSP
445:    rather than a direct method.

447:    An alternative to pass an approximation of A-sigma*B with this function is
448:    to provide approximations of A and B via STSetSplitPreconditioner(). The
449:    difference is that when sigma changes the preconditioner is recomputed.

451:    Use NULL to remove a previously set matrix.

453:    Level: advanced

455: .seealso: STGetPreconditionerMat(), STSetShift(), STGetOperator(), STSetSplitPreconditioner()
456: @*/
457: PetscErrorCode STSetPreconditionerMat(ST st,Mat mat)
458: {
459:   PetscFunctionBegin;
461:   if (mat) {
463:     PetscCheckSameComm(st,1,mat,2);
464:   }
465:   STCheckNotSeized(st,1);
466:   PetscCheck(!mat || !st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot call both STSetPreconditionerMat and STSetSplitPreconditioner");
467:   if (mat) PetscCall(PetscObjectReference((PetscObject)mat));
468:   PetscCall(MatDestroy(&st->Pmat));
469:   st->Pmat     = mat;
470:   st->Pmat_set = mat? PETSC_TRUE: PETSC_FALSE;
471:   st->state    = ST_STATE_INITIAL;
472:   st->opready  = PETSC_FALSE;
473:   PetscFunctionReturn(PETSC_SUCCESS);
474: }

476: /*@
477:    STGetPreconditionerMat - Returns the matrix previously set by STSetPreconditionerMat().

479:    Not Collective

481:    Input Parameter:
482: .  st - the spectral transformation context

484:    Output Parameter:
485: .  mat - the matrix that will be used in constructing the preconditioner or
486:    NULL if no matrix was set by STSetPreconditionerMat().

488:    Level: advanced

490: .seealso: STSetPreconditionerMat()
491: @*/
492: PetscErrorCode STGetPreconditionerMat(ST st,Mat *mat)
493: {
494:   PetscFunctionBegin;
496:   PetscAssertPointer(mat,2);
497:   *mat = st->Pmat_set? st->Pmat: NULL;
498:   PetscFunctionReturn(PETSC_SUCCESS);
499: }

501: /*@
502:    STSetSplitPreconditioner - Sets the matrices from which to build the preconditioner
503:    in split form.

505:    Collective

507:    Input Parameters:
508: +  st     - the spectral transformation context
509: .  n      - number of matrices
510: .  Psplit - array of matrices
511: -  strp   - structure flag for Psplit matrices

513:    Notes:
514:    The number of matrices passed here must be the same as in STSetMatrices().

516:    For linear eigenproblems, the preconditioner matrix is computed as
517:    Pmat(sigma) = A0-sigma*B0, where A0 and B0 are approximations of A and B
518:    (the eigenproblem matrices) provided via the Psplit array in this function.
519:    Compared to STSetPreconditionerMat(), this function allows setting a preconditioner
520:    in a way that is independent of the shift sigma. Whenever the value of sigma
521:    changes the preconditioner is recomputed.

523:    Similarly, for polynomial eigenproblems the matrix for the preconditioner
524:    is expressed as Pmat(sigma) = sum_i Psplit_i*phi_i(sigma), for i=1,...,n, where
525:    the phi_i's are the polynomial basis functions.

527:    The structure flag provides information about the relative nonzero pattern of the
528:    Psplit_i matrices, in the same way as in STSetMatStructure().

530:    Use n=0 to reset a previously set split preconditioner.

532:    Level: advanced

534: .seealso: STGetSplitPreconditionerTerm(), STGetSplitPreconditionerInfo(), STSetPreconditionerMat(), STSetMatrices(), STSetMatStructure()
535: @*/
536: PetscErrorCode STSetSplitPreconditioner(ST st,PetscInt n,Mat Psplit[],MatStructure strp)
537: {
538:   PetscInt       i,N=0,M,M0=0,mloc,nloc,mloc0=0;

540:   PetscFunctionBegin;
543:   PetscCheck(n>=0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of n = %" PetscInt_FMT,n);
544:   PetscCheck(!n || !st->Pmat_set,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot call both STSetPreconditionerMat and STSetSplitPreconditioner");
545:   PetscCheck(!n || !st->nmat || st->nmat==n,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"The number of matrices must be the same as in STSetMatrices()");
546:   if (n) PetscAssertPointer(Psplit,3);
548:   STCheckNotSeized(st,1);

550:   for (i=0;i<n;i++) {
552:     PetscCheckSameComm(st,1,Psplit[i],3);
553:     PetscCall(MatGetSize(Psplit[i],&M,&N));
554:     PetscCall(MatGetLocalSize(Psplit[i],&mloc,&nloc));
555:     PetscCheck(M==N,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Psplit[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,M,N);
556:     PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Psplit[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
557:     if (!i) { M0 = M; mloc0 = mloc; }
558:     PetscCheck(M==M0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_INCOMP,"Dimensions of Psplit[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,M,M0);
559:     PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_INCOMP,"Local dimensions of Psplit[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
560:     PetscCall(PetscObjectReference((PetscObject)Psplit[i]));
561:   }

563:   if (st->Psplit) PetscCall(MatDestroyMatrices(st->nsplit,&st->Psplit));

565:   /* allocate space and copy matrices */
566:   if (n) {
567:     PetscCall(PetscMalloc1(n,&st->Psplit));
568:     for (i=0;i<n;i++) st->Psplit[i] = Psplit[i];
569:   }
570:   st->nsplit = n;
571:   st->strp   = strp;
572:   st->state  = ST_STATE_INITIAL;
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

576: /*@
577:    STGetSplitPreconditionerTerm - Gets the matrices associated with
578:    the split preconditioner.

580:    Not Collective

582:    Input Parameters:
583: +  st - the spectral transformation context
584: -  k  - the index of the requested matrix (starting in 0)

586:    Output Parameter:
587: .  Psplit - the returned matrix

589:    Level: advanced

591: .seealso: STSetSplitPreconditioner(), STGetSplitPreconditionerInfo()
592: @*/
593: PetscErrorCode STGetSplitPreconditionerTerm(ST st,PetscInt k,Mat *Psplit)
594: {
595:   PetscFunctionBegin;
598:   PetscAssertPointer(Psplit,3);
599:   PetscCheck(k>=0 && k<st->nsplit,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nsplit-1);
600:   PetscCheck(st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"You have not called STSetSplitPreconditioner()");
601:   *Psplit = st->Psplit[k];
602:   PetscFunctionReturn(PETSC_SUCCESS);
603: }

605: /*@
606:    STGetSplitPreconditionerInfo - Returns the number of matrices of the split
607:    preconditioner, as well as the structure flag.

609:    Not Collective

611:    Input Parameter:
612: .  st - the spectral transformation context

614:    Output Parameters:
615: +  n    - the number of matrices passed in STSetSplitPreconditioner()
616: -  strp - the matrix structure flag passed in STSetSplitPreconditioner()

618:    Level: advanced

620: .seealso: STSetSplitPreconditioner(), STGetSplitPreconditionerTerm()
621: @*/
622: PetscErrorCode STGetSplitPreconditionerInfo(ST st,PetscInt *n,MatStructure *strp)
623: {
624:   PetscFunctionBegin;
626:   if (n)    *n    = st->nsplit;
627:   if (strp) *strp = st->strp;
628:   PetscFunctionReturn(PETSC_SUCCESS);
629: }

631: /*@
632:    STSetShift - Sets the shift associated with the spectral transformation.

634:    Collective

636:    Input Parameters:
637: +  st - the spectral transformation context
638: -  shift - the value of the shift

640:    Notes:
641:    In some spectral transformations, changing the shift may have associated
642:    a lot of work, for example recomputing a factorization.

644:    This function is normally not directly called by users, since the shift is
645:    indirectly set by EPSSetTarget().

647:    Level: intermediate

649: .seealso: EPSSetTarget(), STGetShift(), STSetDefaultShift()
650: @*/
651: PetscErrorCode STSetShift(ST st,PetscScalar shift)
652: {
653:   PetscFunctionBegin;
657:   if (st->sigma != shift) {
658:     STCheckNotSeized(st,1);
659:     if (st->state==ST_STATE_SETUP) PetscTryTypeMethod(st,setshift,shift);
660:     st->sigma = shift;
661:   }
662:   st->sigma_set = PETSC_TRUE;
663:   PetscFunctionReturn(PETSC_SUCCESS);
664: }

666: /*@
667:    STGetShift - Gets the shift associated with the spectral transformation.

669:    Not Collective

671:    Input Parameter:
672: .  st - the spectral transformation context

674:    Output Parameter:
675: .  shift - the value of the shift

677:    Level: intermediate

679: .seealso: STSetShift()
680: @*/
681: PetscErrorCode STGetShift(ST st,PetscScalar* shift)
682: {
683:   PetscFunctionBegin;
685:   PetscAssertPointer(shift,2);
686:   *shift = st->sigma;
687:   PetscFunctionReturn(PETSC_SUCCESS);
688: }

690: /*@
691:    STSetDefaultShift - Sets the value of the shift that should be employed if
692:    the user did not specify one.

694:    Logically Collective

696:    Input Parameters:
697: +  st - the spectral transformation context
698: -  defaultshift - the default value of the shift

700:    Level: developer

702: .seealso: STSetShift()
703: @*/
704: PetscErrorCode STSetDefaultShift(ST st,PetscScalar defaultshift)
705: {
706:   PetscFunctionBegin;
709:   if (st->defsigma != defaultshift) {
710:     st->defsigma = defaultshift;
711:     st->state    = ST_STATE_INITIAL;
712:     st->opready  = PETSC_FALSE;
713:   }
714:   PetscFunctionReturn(PETSC_SUCCESS);
715: }

717: /*@
718:    STScaleShift - Multiply the shift with a given factor.

720:    Logically Collective

722:    Input Parameters:
723: +  st     - the spectral transformation context
724: -  factor - the scaling factor

726:    Note:
727:    This function does not update the transformation matrices, as opposed to
728:    STSetShift().

730:    Level: developer

732: .seealso: STSetShift()
733: @*/
734: PetscErrorCode STScaleShift(ST st,PetscScalar factor)
735: {
736:   PetscFunctionBegin;
739:   st->sigma *= factor;
740:   PetscFunctionReturn(PETSC_SUCCESS);
741: }

743: /*@
744:    STSetBalanceMatrix - Sets the diagonal matrix to be used for balancing.

746:    Collective

748:    Input Parameters:
749: +  st - the spectral transformation context
750: -  D  - the diagonal matrix (represented as a vector)

752:    Notes:
753:    If this matrix is set, STApply will effectively apply D*OP*D^{-1}. Use NULL
754:    to reset a previously passed D.

756:    Balancing is usually set via EPSSetBalance, but the advanced user may use
757:    this function to bypass the usual balancing methods.

759:    Level: developer

761: .seealso: EPSSetBalance(), STApply(), STGetBalanceMatrix()
762: @*/
763: PetscErrorCode STSetBalanceMatrix(ST st,Vec D)
764: {
765:   PetscFunctionBegin;
767:   if (st->D == D) PetscFunctionReturn(PETSC_SUCCESS);
768:   STCheckNotSeized(st,1);
769:   if (D) {
771:     PetscCheckSameComm(st,1,D,2);
772:     PetscCall(PetscObjectReference((PetscObject)D));
773:   }
774:   PetscCall(VecDestroy(&st->D));
775:   st->D = D;
776:   st->state   = ST_STATE_INITIAL;
777:   st->opready = PETSC_FALSE;
778:   PetscFunctionReturn(PETSC_SUCCESS);
779: }

781: /*@
782:    STGetBalanceMatrix - Gets the balance matrix used by the spectral transformation.

784:    Not Collective

786:    Input Parameter:
787: .  st - the spectral transformation context

789:    Output Parameter:
790: .  D  - the diagonal matrix (represented as a vector)

792:    Note:
793:    If the matrix was not set, a null pointer will be returned.

795:    Level: developer

797: .seealso: STSetBalanceMatrix()
798: @*/
799: PetscErrorCode STGetBalanceMatrix(ST st,Vec *D)
800: {
801:   PetscFunctionBegin;
803:   PetscAssertPointer(D,2);
804:   *D = st->D;
805:   PetscFunctionReturn(PETSC_SUCCESS);
806: }

808: /*@C
809:    STMatCreateVecs - Get vector(s) compatible with the ST matrices.

811:    Collective

813:    Input Parameter:
814: .  st - the spectral transformation context

816:    Output Parameters:
817: +  right - (optional) vector that the matrix can be multiplied against
818: -  left  - (optional) vector that the matrix vector product can be stored in

820:    Level: developer

822: .seealso: STMatCreateVecsEmpty()
823: @*/
824: PetscErrorCode STMatCreateVecs(ST st,Vec *right,Vec *left)
825: {
826:   PetscFunctionBegin;
827:   STCheckMatrices(st,1);
828:   PetscCall(MatCreateVecs(st->A[0],right,left));
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

832: /*@C
833:    STMatCreateVecsEmpty - Get vector(s) compatible with the ST matrices, i.e. with the same
834:    parallel layout, but without internal array.

836:    Collective

838:    Input Parameter:
839: .  st - the spectral transformation context

841:    Output Parameters:
842: +  right - (optional) vector that the matrix can be multiplied against
843: -  left  - (optional) vector that the matrix vector product can be stored in

845:    Level: developer

847: .seealso: STMatCreateVecs(), MatCreateVecsEmpty()
848: @*/
849: PetscErrorCode STMatCreateVecsEmpty(ST st,Vec *right,Vec *left)
850: {
851:   PetscFunctionBegin;
852:   STCheckMatrices(st,1);
853:   PetscCall(MatCreateVecsEmpty(st->A[0],right,left));
854:   PetscFunctionReturn(PETSC_SUCCESS);
855: }

857: /*@
858:    STMatGetSize - Returns the number of rows and columns of the ST matrices.

860:    Not Collective

862:    Input Parameter:
863: .  st - the spectral transformation context

865:    Output Parameters:
866: +  m - the number of global rows
867: -  n - the number of global columns

869:    Level: developer

871: .seealso: STMatGetLocalSize()
872: @*/
873: PetscErrorCode STMatGetSize(ST st,PetscInt *m,PetscInt *n)
874: {
875:   PetscFunctionBegin;
876:   STCheckMatrices(st,1);
877:   PetscCall(MatGetSize(st->A[0],m,n));
878:   PetscFunctionReturn(PETSC_SUCCESS);
879: }

881: /*@
882:    STMatGetLocalSize - Returns the number of local rows and columns of the ST matrices.

884:    Not Collective

886:    Input Parameter:
887: .  st - the spectral transformation context

889:    Output Parameters:
890: +  m - the number of local rows
891: -  n - the number of local columns

893:    Level: developer

895: .seealso: STMatGetSize()
896: @*/
897: PetscErrorCode STMatGetLocalSize(ST st,PetscInt *m,PetscInt *n)
898: {
899:   PetscFunctionBegin;
900:   STCheckMatrices(st,1);
901:   PetscCall(MatGetLocalSize(st->A[0],m,n));
902:   PetscFunctionReturn(PETSC_SUCCESS);
903: }

905: /*@C
906:    STSetOptionsPrefix - Sets the prefix used for searching for all
907:    ST options in the database.

909:    Logically Collective

911:    Input Parameters:
912: +  st     - the spectral transformation context
913: -  prefix - the prefix string to prepend to all ST option requests

915:    Notes:
916:    A hyphen (-) must NOT be given at the beginning of the prefix name.
917:    The first character of all runtime options is AUTOMATICALLY the
918:    hyphen.

920:    Level: advanced

922: .seealso: STAppendOptionsPrefix(), STGetOptionsPrefix()
923: @*/
924: PetscErrorCode STSetOptionsPrefix(ST st,const char *prefix)
925: {
926:   PetscFunctionBegin;
928:   if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
929:   PetscCall(KSPSetOptionsPrefix(st->ksp,prefix));
930:   PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
931:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)st,prefix));
932:   PetscFunctionReturn(PETSC_SUCCESS);
933: }

935: /*@C
936:    STAppendOptionsPrefix - Appends to the prefix used for searching for all
937:    ST options in the database.

939:    Logically Collective

941:    Input Parameters:
942: +  st     - the spectral transformation context
943: -  prefix - the prefix string to prepend to all ST option requests

945:    Notes:
946:    A hyphen (-) must NOT be given at the beginning of the prefix name.
947:    The first character of all runtime options is AUTOMATICALLY the
948:    hyphen.

950:    Level: advanced

952: .seealso: STSetOptionsPrefix(), STGetOptionsPrefix()
953: @*/
954: PetscErrorCode STAppendOptionsPrefix(ST st,const char *prefix)
955: {
956:   PetscFunctionBegin;
958:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)st,prefix));
959:   if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
960:   PetscCall(KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix));
961:   PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
962:   PetscFunctionReturn(PETSC_SUCCESS);
963: }

965: /*@C
966:    STGetOptionsPrefix - Gets the prefix used for searching for all
967:    ST options in the database.

969:    Not Collective

971:    Input Parameters:
972: .  st - the spectral transformation context

974:    Output Parameters:
975: .  prefix - pointer to the prefix string used, is returned

977:    Note:
978:    On the Fortran side, the user should pass in a string 'prefix' of
979:    sufficient length to hold the prefix.

981:    Level: advanced

983: .seealso: STSetOptionsPrefix(), STAppendOptionsPrefix()
984: @*/
985: PetscErrorCode STGetOptionsPrefix(ST st,const char *prefix[])
986: {
987:   PetscFunctionBegin;
989:   PetscAssertPointer(prefix,2);
990:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)st,prefix));
991:   PetscFunctionReturn(PETSC_SUCCESS);
992: }

994: /*@C
995:    STView - Prints the ST data structure.

997:    Collective

999:    Input Parameters:
1000: +  st - the ST context
1001: -  viewer - optional visualization context

1003:    Note:
1004:    The available visualization contexts include
1005: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1006: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1007:          output where only the first processor opens
1008:          the file.  All other processors send their
1009:          data to the first processor to print.

1011:    The user can open an alternative visualization contexts with
1012:    PetscViewerASCIIOpen() (output to a specified file).

1014:    Level: beginner

1016: .seealso: EPSView()
1017: @*/
1018: PetscErrorCode STView(ST st,PetscViewer viewer)
1019: {
1020:   STType         cstr;
1021:   char           str[50];
1022:   PetscBool      isascii,isstring;

1024:   PetscFunctionBegin;
1026:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)st),&viewer));
1028:   PetscCheckSameComm(st,1,viewer,2);

1030:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1031:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring));
1032:   if (isascii) {
1033:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)st,viewer));
1034:     PetscCall(PetscViewerASCIIPushTab(viewer));
1035:     PetscTryTypeMethod(st,view,viewer);
1036:     PetscCall(PetscViewerASCIIPopTab(viewer));
1037:     PetscCall(SlepcSNPrintfScalar(str,sizeof(str),st->sigma,PETSC_FALSE));
1038:     PetscCall(PetscViewerASCIIPrintf(viewer,"  shift: %s\n",str));
1039:     PetscCall(PetscViewerASCIIPrintf(viewer,"  number of matrices: %" PetscInt_FMT "\n",st->nmat));
1040:     switch (st->matmode) {
1041:     case ST_MATMODE_COPY:
1042:       break;
1043:     case ST_MATMODE_INPLACE:
1044:       PetscCall(PetscViewerASCIIPrintf(viewer,"  shifting the matrix and unshifting at exit\n"));
1045:       break;
1046:     case ST_MATMODE_SHELL:
1047:       PetscCall(PetscViewerASCIIPrintf(viewer,"  using a shell matrix\n"));
1048:       break;
1049:     }
1050:     if (st->nmat>1 && st->matmode != ST_MATMODE_SHELL) PetscCall(PetscViewerASCIIPrintf(viewer,"  nonzero pattern of the matrices: %s\n",MatStructures[st->str]));
1051:     if (st->Psplit) PetscCall(PetscViewerASCIIPrintf(viewer,"  using split preconditioner matrices with %s\n",MatStructures[st->strp]));
1052:     if (st->transform && st->nmat>2) PetscCall(PetscViewerASCIIPrintf(viewer,"  computing transformed matrices\n"));
1053:   } else if (isstring) {
1054:     PetscCall(STGetType(st,&cstr));
1055:     PetscCall(PetscViewerStringSPrintf(viewer," %-7.7s",cstr));
1056:     PetscTryTypeMethod(st,view,viewer);
1057:   }
1058:   if (st->usesksp) {
1059:     if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
1060:     PetscCall(PetscViewerASCIIPushTab(viewer));
1061:     PetscCall(KSPView(st->ksp,viewer));
1062:     PetscCall(PetscViewerASCIIPopTab(viewer));
1063:   }
1064:   PetscFunctionReturn(PETSC_SUCCESS);
1065: }

1067: /*@C
1068:    STViewFromOptions - View from options

1070:    Collective

1072:    Input Parameters:
1073: +  st   - the spectral transformation context
1074: .  obj  - optional object
1075: -  name - command line option

1077:    Level: intermediate

1079: .seealso: STView(), STCreate()
1080: @*/
1081: PetscErrorCode STViewFromOptions(ST st,PetscObject obj,const char name[])
1082: {
1083:   PetscFunctionBegin;
1085:   PetscCall(PetscObjectViewFromOptions((PetscObject)st,obj,name));
1086:   PetscFunctionReturn(PETSC_SUCCESS);
1087: }

1089: /*@C
1090:    STRegister - Adds a method to the spectral transformation package.

1092:    Not Collective

1094:    Input Parameters:
1095: +  name - name of a new user-defined transformation
1096: -  function - routine to create method context

1098:    Notes:
1099:    STRegister() may be called multiple times to add several user-defined
1100:    spectral transformations.

1102:    Example Usage:
1103: .vb
1104:     STRegister("my_transform",MyTransformCreate);
1105: .ve

1107:    Then, your spectral transform can be chosen with the procedural interface via
1108: $     STSetType(st,"my_transform")
1109:    or at runtime via the option
1110: $     -st_type my_transform

1112:    Level: advanced

1114: .seealso: STRegisterAll()
1115: @*/
1116: PetscErrorCode STRegister(const char *name,PetscErrorCode (*function)(ST))
1117: {
1118:   PetscFunctionBegin;
1119:   PetscCall(STInitializePackage());
1120:   PetscCall(PetscFunctionListAdd(&STList,name,function));
1121:   PetscFunctionReturn(PETSC_SUCCESS);
1122: }