Actual source code: stfunc.c

slepc-3.20.2 2024-03-15
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_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("STMatSetUp",ST_CLASSID,&ST_MatSetUp));
 67:   PetscCall(PetscLogEventRegister("STMatMult",ST_CLASSID,&ST_MatMult));
 68:   PetscCall(PetscLogEventRegister("STMatMultTranspose",ST_CLASSID,&ST_MatMultTranspose));
 69:   PetscCall(PetscLogEventRegister("STMatSolve",ST_CLASSID,&ST_MatSolve));
 70:   PetscCall(PetscLogEventRegister("STMatSolveTranspose",ST_CLASSID,&ST_MatSolveTranspose));
 71:   /* Process Info */
 72:   classids[0] = ST_CLASSID;
 73:   PetscCall(PetscInfoProcessClass("st",1,&classids[0]));
 74:   /* Process summary exclusions */
 75:   PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
 76:   if (opt) {
 77:     PetscCall(PetscStrInList("st",logList,',',&pkg));
 78:     if (pkg) PetscCall(PetscLogEventDeactivateClass(ST_CLASSID));
 79:   }
 80:   /* Register package finalizer */
 81:   PetscCall(PetscRegisterFinalize(STFinalizePackage));
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

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

 89:    Collective

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

 94:    Level: advanced

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

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

128:    Collective

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

133:    Level: beginner

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

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

153:    Collective

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

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

161:    Level: beginner

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

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

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

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

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

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

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

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

249:    Collective

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

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

260:    Level: intermediate

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

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

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

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

314:    Not Collective

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

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

323:    Level: intermediate

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

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

343:    Not Collective

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

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

352:    Level: developer

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

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

372:    Not Collective

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

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

380:    Level: intermediate

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

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

396:    Logically Collective

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

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

407:    Level: developer

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

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

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

424:    Collective

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

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

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

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

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

450:    Use NULL to remove a previously set matrix.

452:    Level: advanced

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

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

478:    Not Collective

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

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

487:    Level: advanced

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

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

504:    Collective

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

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

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

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

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

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

531:    Level: advanced

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

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

549:   for (i=0;i<n;i++) {
551:     PetscCheckSameComm(st,1,Psplit[i],3);
552:     PetscCall(MatGetSize(Psplit[i],&M,&N));
553:     PetscCall(MatGetLocalSize(Psplit[i],&mloc,&nloc));
554:     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);
555:     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);
556:     if (!i) { M0 = M; mloc0 = mloc; }
557:     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);
558:     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);
559:     PetscCall(PetscObjectReference((PetscObject)Psplit[i]));
560:   }

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

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

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

579:    Not Collective

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

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

588:    Level: advanced

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

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

608:    Not Collective

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

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

617:    Level: advanced

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

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

633:    Collective

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

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

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

646:    Level: intermediate

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

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

668:    Not Collective

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

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

676:    Level: intermediate

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

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

693:    Logically Collective

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

699:    Level: developer

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

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

719:    Logically Collective

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

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

729:    Level: developer

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

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

745:    Collective

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

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

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

758:    Level: developer

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

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

783:    Not Collective

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

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

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

794:    Level: developer

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

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

810:    Collective

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

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

819:    Level: developer

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

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

835:    Collective

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

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

844:    Level: developer

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

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

859:    Not Collective

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

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

868:    Level: developer

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

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

883:    Not Collective

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

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

892:    Level: developer

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

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

908:    Logically Collective

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

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

919:    Level: advanced

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

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

938:    Logically Collective

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

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

949:    Level: advanced

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

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

968:    Not Collective

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

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

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

980:    Level: advanced

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

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

996:    Collective

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

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

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

1013:    Level: beginner

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

1023:   PetscFunctionBegin;
1025:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)st),&viewer));
1027:   PetscCheckSameComm(st,1,viewer,2);

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

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

1069:    Collective

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

1076:    Level: intermediate

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

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

1091:    Not Collective

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

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

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

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

1111:    Level: advanced

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