Actual source code: stfunc.c

slepc-3.17.1 2022-04-11
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_",0};

 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:   PetscFunctionListDestroy(&STList);
 33:   STPackageInitialized = PETSC_FALSE;
 34:   STRegisterAllCalled  = PETSC_FALSE;
 35:   PetscFunctionReturn(0);
 36: }

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

 43:    Level: developer

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

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

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

 87:    Collective on st

 89:    Input Parameter:
 90: .  st - the spectral transformation context

 92:    Level: advanced

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

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

125:    Collective on st

127:    Input Parameter:
128: .  st - the spectral transformation context

130:    Level: beginner

132: .seealso: STCreate(), STSetUp()
133: @*/
134: PetscErrorCode STDestroy(ST *st)
135: {
136:   if (!*st) PetscFunctionReturn(0);
138:   if (--((PetscObject)(*st))->refct > 0) { *st = 0; PetscFunctionReturn(0); }
139:   STReset(*st);
140:   if ((*st)->ops->destroy) (*(*st)->ops->destroy)(*st);
141:   KSPDestroy(&(*st)->ksp);
142:   PetscHeaderDestroy(st);
143:   PetscFunctionReturn(0);
144: }

146: /*@
147:    STCreate - Creates a spectral transformation context.

149:    Collective

151:    Input Parameter:
152: .  comm - MPI communicator

154:    Output Parameter:
155: .  newst - location to put the spectral transformation context

157:    Level: beginner

159: .seealso: STSetUp(), STApply(), STDestroy(), ST
160: @*/
161: PetscErrorCode STCreate(MPI_Comm comm,ST *newst)
162: {
163:   ST             st;

166:   *newst = 0;
167:   STInitializePackage();
168:   SlepcHeaderCreate(st,ST_CLASSID,"ST","Spectral Transformation","ST",comm,STDestroy,STView);

170:   st->A            = NULL;
171:   st->nmat         = 0;
172:   st->sigma        = 0.0;
173:   st->defsigma     = 0.0;
174:   st->matmode      = ST_MATMODE_COPY;
175:   st->str          = UNKNOWN_NONZERO_PATTERN;
176:   st->transform    = PETSC_FALSE;
177:   st->D            = NULL;
178:   st->Pmat         = NULL;
179:   st->Pmat_set     = PETSC_FALSE;
180:   st->Psplit       = NULL;
181:   st->nsplit       = 0;
182:   st->strp         = UNKNOWN_NONZERO_PATTERN;

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

203:   *newst = st;
204:   PetscFunctionReturn(0);
205: }

207: /*
208:    Checks whether the ST matrices are all symmetric or hermitian.
209: */
210: static inline PetscErrorCode STMatIsSymmetricKnown(ST st,PetscBool *symm,PetscBool *herm)
211: {
212:   PetscInt       i;
213:   PetscBool      sbaij=PETSC_FALSE,set,flg=PETSC_FALSE;

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

240: /*@
241:    STSetMatrices - Sets the matrices associated with the eigenvalue problem.

243:    Collective on st

245:    Input Parameters:
246: +  st - the spectral transformation context
247: .  n  - number of matrices in array A
248: -  A  - the array of matrices associated with the eigensystem

250:    Notes:
251:    It must be called before STSetUp(). If it is called again after STSetUp() then
252:    the ST object is reset.

254:    Level: intermediate

256: .seealso: STGetMatrix(), STGetNumMatrices(), STSetUp(), STReset()
257:  @*/
258: PetscErrorCode STSetMatrices(ST st,PetscInt n,Mat A[])
259: {
260:   PetscInt       i;
261:   PetscBool      same=PETSC_TRUE;

268:   STCheckNotSeized(st,1);

271:   if (st->state) {
272:     if (n!=st->nmat) same = PETSC_FALSE;
273:     for (i=0;same&&i<n;i++) {
274:       if (A[i]!=st->A[i]) same = PETSC_FALSE;
275:     }
276:     if (!same) STReset(st);
277:   } else same = PETSC_FALSE;
278:   if (!same) {
279:     MatDestroyMatrices(PetscMax(2,st->nmat),&st->A);
280:     PetscCalloc1(PetscMax(2,n),&st->A);
281:     PetscLogObjectMemory((PetscObject)st,PetscMax(2,n)*sizeof(Mat));
282:     PetscFree(st->Astate);
283:     PetscMalloc1(PetscMax(2,n),&st->Astate);
284:     PetscLogObjectMemory((PetscObject)st,PetscMax(2,n)*sizeof(PetscInt));
285:   }
286:   for (i=0;i<n;i++) {
288:     PetscObjectReference((PetscObject)A[i]);
289:     MatDestroy(&st->A[i]);
290:     st->A[i] = A[i];
291:     st->Astate[i] = ((PetscObject)A[i])->state;
292:   }
293:   if (n==1) {
294:     st->A[1] = NULL;
295:     st->Astate[1] = 0;
296:   }
297:   st->nmat = n;
298:   if (same) st->state = ST_STATE_UPDATED;
299:   else st->state = ST_STATE_INITIAL;
301:   st->opready = PETSC_FALSE;
302:   if (!same) STMatIsSymmetricKnown(st,&st->asymm,&st->aherm);
303:   PetscFunctionReturn(0);
304: }

306: /*@
307:    STGetMatrix - Gets the matrices associated with the original eigensystem.

309:    Not collective, though parallel Mats are returned if the ST is parallel

311:    Input Parameters:
312: +  st - the spectral transformation context
313: -  k  - the index of the requested matrix (starting in 0)

315:    Output Parameters:
316: .  A - the requested matrix

318:    Level: intermediate

320: .seealso: STSetMatrices(), STGetNumMatrices()
321: @*/
322: PetscErrorCode STGetMatrix(ST st,PetscInt k,Mat *A)
323: {
327:   STCheckMatrices(st,1);
330:   *A = st->A[k];
331:   PetscFunctionReturn(0);
332: }

334: /*@
335:    STGetMatrixTransformed - Gets the matrices associated with the transformed eigensystem.

337:    Not collective, though parallel Mats are returned if the ST is parallel

339:    Input Parameters:
340: +  st - the spectral transformation context
341: -  k  - the index of the requested matrix (starting in 0)

343:    Output Parameters:
344: .  T - the requested matrix

346:    Level: developer

348: .seealso: STGetMatrix(), STGetNumMatrices()
349: @*/
350: PetscErrorCode STGetMatrixTransformed(ST st,PetscInt k,Mat *T)
351: {
355:   STCheckMatrices(st,1);
358:   *T = st->T[k];
359:   PetscFunctionReturn(0);
360: }

362: /*@
363:    STGetNumMatrices - Returns the number of matrices stored in the ST.

365:    Not collective

367:    Input Parameter:
368: .  st - the spectral transformation context

370:    Output Parameters:
371: .  n - the number of matrices passed in STSetMatrices()

373:    Level: intermediate

375: .seealso: STSetMatrices()
376: @*/
377: PetscErrorCode STGetNumMatrices(ST st,PetscInt *n)
378: {
381:   *n = st->nmat;
382:   PetscFunctionReturn(0);
383: }

385: /*@
386:    STResetMatrixState - Resets the stored state of the matrices in the ST.

388:    Logically Collective on st

390:    Input Parameter:
391: .  st - the spectral transformation context

393:    Note:
394:    This is useful in solvers where the user matrices are modified during
395:    the computation, as in nonlinear inverse iteration. The effect is that
396:    STGetMatrix() will retrieve the modified matrices as if they were
397:    the matrices originally provided by the user.

399:    Level: developer

401: .seealso: STGetMatrix(), EPSPowerSetNonlinear()
402: @*/
403: PetscErrorCode STResetMatrixState(ST st)
404: {
405:   PetscInt i;

408:   for (i=0;i<st->nmat;i++) st->Astate[i] = ((PetscObject)st->A[i])->state;
409:   PetscFunctionReturn(0);
410: }

412: /*@
413:    STSetPreconditionerMat - Sets the matrix to be used to build the preconditioner.

415:    Collective on st

417:    Input Parameters:
418: +  st  - the spectral transformation context
419: -  mat - the matrix that will be used in constructing the preconditioner

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

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

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

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

441:    Use NULL to remove a previously set matrix.

443:    Level: advanced

445: .seealso: STGetPreconditionerMat(), STSetShift(), STGetOperator(), STSetSplitPreconditioner()
446: @*/
447: PetscErrorCode STSetPreconditionerMat(ST st,Mat mat)
448: {
450:   if (mat) {
453:   }
454:   STCheckNotSeized(st,1);
456:   if (mat) PetscObjectReference((PetscObject)mat);
457:   MatDestroy(&st->Pmat);
458:   st->Pmat     = mat;
459:   st->Pmat_set = mat? PETSC_TRUE: PETSC_FALSE;
460:   st->state    = ST_STATE_INITIAL;
461:   st->opready  = PETSC_FALSE;
462:   PetscFunctionReturn(0);
463: }

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

468:    Not Collective

470:    Input Parameter:
471: .  st - the spectral transformation context

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

477:    Level: advanced

479: .seealso: STSetPreconditionerMat()
480: @*/
481: PetscErrorCode STGetPreconditionerMat(ST st,Mat *mat)
482: {
485:   *mat = st->Pmat_set? st->Pmat: NULL;
486:   PetscFunctionReturn(0);
487: }

489: /*@
490:    STSetSplitPreconditioner - Sets the matrices from which to build the preconditioner
491:    in split form.

493:    Collective on st

495:    Input Parameters:
496: +  st     - the spectral transformation context
497: .  n      - number of matrices
498: .  Psplit - array of matrices
499: -  strp   - structure flag for Psplit matrices

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

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

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

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

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

520:    Level: advanced

522: .seealso: STGetSplitPreconditionerTerm(), STGetSplitPreconditionerInfo(), STSetPreconditionerMat(), STSetMatrices(), STSetMatStructure()
523:  @*/
524: PetscErrorCode STSetSplitPreconditioner(ST st,PetscInt n,Mat Psplit[],MatStructure strp)
525: {
526:   PetscInt       i,N=0,M,M0=0,mloc,nloc,mloc0=0;

535:   STCheckNotSeized(st,1);

537:   for (i=0;i<n;i++) {
540:     MatGetSize(Psplit[i],&M,&N);
541:     MatGetLocalSize(Psplit[i],&mloc,&nloc);
544:     if (!i) { M0 = M; mloc0 = mloc; }
547:     PetscObjectReference((PetscObject)Psplit[i]);
548:   }

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

552:   /* allocate space and copy matrices */
553:   if (n) {
554:     PetscMalloc1(n,&st->Psplit);
555:     PetscLogObjectMemory((PetscObject)st,n*sizeof(Mat));
556:     for (i=0;i<n;i++) st->Psplit[i] = Psplit[i];
557:   }
558:   st->nsplit = n;
559:   st->strp   = strp;
560:   st->state  = ST_STATE_INITIAL;
561:   PetscFunctionReturn(0);
562: }

564: /*@
565:    STGetSplitPreconditionerTerm - Gets the matrices associated with
566:    the split preconditioner.

568:    Not collective, though parallel Mats are returned if the ST is parallel

570:    Input Parameters:
571: +  st - the spectral transformation context
572: -  k  - the index of the requested matrix (starting in 0)

574:    Output Parameter:
575: .  Psplit - the returned matrix

577:    Level: advanced

579: .seealso: STSetSplitPreconditioner(), STGetSplitPreconditionerInfo()
580: @*/
581: PetscErrorCode STGetSplitPreconditionerTerm(ST st,PetscInt k,Mat *Psplit)
582: {
588:   *Psplit = st->Psplit[k];
589:   PetscFunctionReturn(0);
590: }

592: /*@
593:    STGetSplitPreconditionerInfo - Returns the number of matrices of the split
594:    preconditioner, as well as the structure flag.

596:    Not collective

598:    Input Parameter:
599: .  st - the spectral transformation context

601:    Output Parameters:
602: +  n    - the number of matrices passed in STSetSplitPreconditioner()
603: -  strp - the matrix structure flag passed in STSetSplitPreconditioner()

605:    Level: advanced

607: .seealso: STSetSplitPreconditioner(), STGetSplitPreconditionerTerm()
608: @*/
609: PetscErrorCode STGetSplitPreconditionerInfo(ST st,PetscInt *n,MatStructure *strp)
610: {
612:   if (n)    *n    = st->nsplit;
613:   if (strp) *strp = st->strp;
614:   PetscFunctionReturn(0);
615: }

617: /*@
618:    STSetShift - Sets the shift associated with the spectral transformation.

620:    Logically Collective on st

622:    Input Parameters:
623: +  st - the spectral transformation context
624: -  shift - the value of the shift

626:    Notes:
627:    In some spectral transformations, changing the shift may have associated
628:    a lot of work, for example recomputing a factorization.

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

633:    Level: intermediate

635: .seealso: EPSSetTarget(), STGetShift(), STSetDefaultShift()
636: @*/
637: PetscErrorCode STSetShift(ST st,PetscScalar shift)
638: {
642:   if (st->sigma != shift) {
643:     STCheckNotSeized(st,1);
644:     if (st->state==ST_STATE_SETUP && st->ops->setshift) (*st->ops->setshift)(st,shift);
645:     st->sigma = shift;
646:   }
647:   st->sigma_set = PETSC_TRUE;
648:   PetscFunctionReturn(0);
649: }

651: /*@
652:    STGetShift - Gets the shift associated with the spectral transformation.

654:    Not Collective

656:    Input Parameter:
657: .  st - the spectral transformation context

659:    Output Parameter:
660: .  shift - the value of the shift

662:    Level: intermediate

664: .seealso: STSetShift()
665: @*/
666: PetscErrorCode STGetShift(ST st,PetscScalar* shift)
667: {
670:   *shift = st->sigma;
671:   PetscFunctionReturn(0);
672: }

674: /*@
675:    STSetDefaultShift - Sets the value of the shift that should be employed if
676:    the user did not specify one.

678:    Logically Collective on st

680:    Input Parameters:
681: +  st - the spectral transformation context
682: -  defaultshift - the default value of the shift

684:    Level: developer

686: .seealso: STSetShift()
687: @*/
688: PetscErrorCode STSetDefaultShift(ST st,PetscScalar defaultshift)
689: {
692:   if (st->defsigma != defaultshift) {
693:     st->defsigma = defaultshift;
694:     st->state    = ST_STATE_INITIAL;
695:     st->opready  = PETSC_FALSE;
696:   }
697:   PetscFunctionReturn(0);
698: }

700: /*@
701:    STScaleShift - Multiply the shift with a given factor.

703:    Logically Collective on st

705:    Input Parameters:
706: +  st     - the spectral transformation context
707: -  factor - the scaling factor

709:    Note:
710:    This function does not update the transformation matrices, as opposed to
711:    STSetShift().

713:    Level: developer

715: .seealso: STSetShift()
716: @*/
717: PetscErrorCode STScaleShift(ST st,PetscScalar factor)
718: {
721:   st->sigma *= factor;
722:   PetscFunctionReturn(0);
723: }

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

728:    Collective on st

730:    Input Parameters:
731: +  st - the spectral transformation context
732: -  D  - the diagonal matrix (represented as a vector)

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

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

741:    Level: developer

743: .seealso: EPSSetBalance(), STApply(), STGetBalanceMatrix()
744: @*/
745: PetscErrorCode STSetBalanceMatrix(ST st,Vec D)
746: {
748:   if (st->D == D) PetscFunctionReturn(0);
749:   STCheckNotSeized(st,1);
750:   if (D) {
753:     PetscObjectReference((PetscObject)D);
754:   }
755:   VecDestroy(&st->D);
756:   st->D = D;
757:   st->state   = ST_STATE_INITIAL;
758:   st->opready = PETSC_FALSE;
759:   PetscFunctionReturn(0);
760: }

762: /*@
763:    STGetBalanceMatrix - Gets the balance matrix used by the spectral transformation.

765:    Not collective, but vector is shared by all processors that share the ST

767:    Input Parameter:
768: .  st - the spectral transformation context

770:    Output Parameter:
771: .  D  - the diagonal matrix (represented as a vector)

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

776:    Level: developer

778: .seealso: STSetBalanceMatrix()
779: @*/
780: PetscErrorCode STGetBalanceMatrix(ST st,Vec *D)
781: {
784:   *D = st->D;
785:   PetscFunctionReturn(0);
786: }

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

791:    Collective on st

793:    Input Parameter:
794: .  st - the spectral transformation context

796:    Output Parameters:
797: +  right - (optional) vector that the matrix can be multiplied against
798: -  left  - (optional) vector that the matrix vector product can be stored in

800:    Level: developer

802: .seealso: STMatCreateVecsEmpty()
803: @*/
804: PetscErrorCode STMatCreateVecs(ST st,Vec *right,Vec *left)
805: {
806:   STCheckMatrices(st,1);
807:   MatCreateVecs(st->A[0],right,left);
808:   PetscFunctionReturn(0);
809: }

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

815:    Collective on st

817:    Input Parameter:
818: .  st - the spectral transformation context

820:    Output Parameters:
821: +  right - (optional) vector that the matrix can be multiplied against
822: -  left  - (optional) vector that the matrix vector product can be stored in

824:    Level: developer

826: .seealso: STMatCreateVecs(), MatCreateVecsEmpty()
827: @*/
828: PetscErrorCode STMatCreateVecsEmpty(ST st,Vec *right,Vec *left)
829: {
830:   STCheckMatrices(st,1);
831:   MatCreateVecsEmpty(st->A[0],right,left);
832:   PetscFunctionReturn(0);
833: }

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

838:    Not Collective

840:    Input Parameter:
841: .  st - the spectral transformation context

843:    Output Parameters:
844: +  m - the number of global rows
845: -  n - the number of global columns

847:    Level: developer

849: .seealso: STMatGetLocalSize()
850: @*/
851: PetscErrorCode STMatGetSize(ST st,PetscInt *m,PetscInt *n)
852: {
853:   STCheckMatrices(st,1);
854:   MatGetSize(st->A[0],m,n);
855:   PetscFunctionReturn(0);
856: }

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

861:    Not Collective

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

866:    Output Parameters:
867: +  m - the number of local rows
868: -  n - the number of local columns

870:    Level: developer

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

881: /*@C
882:    STSetOptionsPrefix - Sets the prefix used for searching for all
883:    ST options in the database.

885:    Logically Collective on st

887:    Input Parameters:
888: +  st     - the spectral transformation context
889: -  prefix - the prefix string to prepend to all ST option requests

891:    Notes:
892:    A hyphen (-) must NOT be given at the beginning of the prefix name.
893:    The first character of all runtime options is AUTOMATICALLY the
894:    hyphen.

896:    Level: advanced

898: .seealso: STAppendOptionsPrefix(), STGetOptionsPrefix()
899: @*/
900: PetscErrorCode STSetOptionsPrefix(ST st,const char *prefix)
901: {
903:   if (!st->ksp) STGetKSP(st,&st->ksp);
904:   KSPSetOptionsPrefix(st->ksp,prefix);
905:   KSPAppendOptionsPrefix(st->ksp,"st_");
906:   PetscObjectSetOptionsPrefix((PetscObject)st,prefix);
907:   PetscFunctionReturn(0);
908: }

910: /*@C
911:    STAppendOptionsPrefix - Appends to the prefix used for searching for all
912:    ST options in the database.

914:    Logically Collective on st

916:    Input Parameters:
917: +  st     - the spectral transformation context
918: -  prefix - the prefix string to prepend to all ST option requests

920:    Notes:
921:    A hyphen (-) must NOT be given at the beginning of the prefix name.
922:    The first character of all runtime options is AUTOMATICALLY the
923:    hyphen.

925:    Level: advanced

927: .seealso: STSetOptionsPrefix(), STGetOptionsPrefix()
928: @*/
929: PetscErrorCode STAppendOptionsPrefix(ST st,const char *prefix)
930: {
932:   PetscObjectAppendOptionsPrefix((PetscObject)st,prefix);
933:   if (!st->ksp) STGetKSP(st,&st->ksp);
934:   KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix);
935:   KSPAppendOptionsPrefix(st->ksp,"st_");
936:   PetscFunctionReturn(0);
937: }

939: /*@C
940:    STGetOptionsPrefix - Gets the prefix used for searching for all
941:    ST options in the database.

943:    Not Collective

945:    Input Parameters:
946: .  st - the spectral transformation context

948:    Output Parameters:
949: .  prefix - pointer to the prefix string used, is returned

951:    Note:
952:    On the Fortran side, the user should pass in a string 'prefix' of
953:    sufficient length to hold the prefix.

955:    Level: advanced

957: .seealso: STSetOptionsPrefix(), STAppendOptionsPrefix()
958: @*/
959: PetscErrorCode STGetOptionsPrefix(ST st,const char *prefix[])
960: {
963:   PetscObjectGetOptionsPrefix((PetscObject)st,prefix);
964:   PetscFunctionReturn(0);
965: }

967: /*@C
968:    STView - Prints the ST data structure.

970:    Collective on st

972:    Input Parameters:
973: +  st - the ST context
974: -  viewer - optional visualization context

976:    Note:
977:    The available visualization contexts include
978: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
979: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
980:          output where only the first processor opens
981:          the file.  All other processors send their
982:          data to the first processor to print.

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

987:    Level: beginner

989: .seealso: EPSView()
990: @*/
991: PetscErrorCode STView(ST st,PetscViewer viewer)
992: {
993:   STType         cstr;
994:   char           str[50];
995:   PetscBool      isascii,isstring;

998:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)st),&viewer);

1002:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1003:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1004:   if (isascii) {
1005:     PetscObjectPrintClassNamePrefixType((PetscObject)st,viewer);
1006:     if (st->ops->view) {
1007:       PetscViewerASCIIPushTab(viewer);
1008:       (*st->ops->view)(st,viewer);
1009:       PetscViewerASCIIPopTab(viewer);
1010:     }
1011:     SlepcSNPrintfScalar(str,sizeof(str),st->sigma,PETSC_FALSE);
1012:     PetscViewerASCIIPrintf(viewer,"  shift: %s\n",str);
1013:     PetscViewerASCIIPrintf(viewer,"  number of matrices: %" PetscInt_FMT "\n",st->nmat);
1014:     switch (st->matmode) {
1015:     case ST_MATMODE_COPY:
1016:       break;
1017:     case ST_MATMODE_INPLACE:
1018:       PetscViewerASCIIPrintf(viewer,"  shifting the matrix and unshifting at exit\n");
1019:       break;
1020:     case ST_MATMODE_SHELL:
1021:       PetscViewerASCIIPrintf(viewer,"  using a shell matrix\n");
1022:       break;
1023:     }
1024:     if (st->nmat>1 && st->matmode != ST_MATMODE_SHELL) PetscViewerASCIIPrintf(viewer,"  nonzero pattern of the matrices: %s\n",MatStructures[st->str]);
1025:     if (st->Psplit) PetscViewerASCIIPrintf(viewer,"  using split preconditioner matrices with %s\n",MatStructures[st->strp]);
1026:     if (st->transform && st->nmat>2) PetscViewerASCIIPrintf(viewer,"  computing transformed matrices\n");
1027:   } else if (isstring) {
1028:     STGetType(st,&cstr);
1029:     PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1030:     if (st->ops->view) (*st->ops->view)(st,viewer);
1031:   }
1032:   if (st->usesksp) {
1033:     if (!st->ksp) STGetKSP(st,&st->ksp);
1034:     PetscViewerASCIIPushTab(viewer);
1035:     KSPView(st->ksp,viewer);
1036:     PetscViewerASCIIPopTab(viewer);
1037:   }
1038:   PetscFunctionReturn(0);
1039: }

1041: /*@C
1042:    STViewFromOptions - View from options

1044:    Collective on ST

1046:    Input Parameters:
1047: +  st   - the spectral transformation context
1048: .  obj  - optional object
1049: -  name - command line option

1051:    Level: intermediate

1053: .seealso: STView(), STCreate()
1054: @*/
1055: PetscErrorCode STViewFromOptions(ST st,PetscObject obj,const char name[])
1056: {
1058:   PetscObjectViewFromOptions((PetscObject)st,obj,name);
1059:   PetscFunctionReturn(0);
1060: }

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

1065:    Not collective

1067:    Input Parameters:
1068: +  name - name of a new user-defined transformation
1069: -  function - routine to create method context

1071:    Notes:
1072:    STRegister() may be called multiple times to add several user-defined
1073:    spectral transformations.

1075:    Sample usage:
1076: .vb
1077:     STRegister("my_transform",MyTransformCreate);
1078: .ve

1080:    Then, your spectral transform can be chosen with the procedural interface via
1081: $     STSetType(st,"my_transform")
1082:    or at runtime via the option
1083: $     -st_type my_transform

1085:    Level: advanced

1087: .seealso: STRegisterAll()
1088: @*/
1089: PetscErrorCode STRegister(const char *name,PetscErrorCode (*function)(ST))
1090: {
1091:   STInitializePackage();
1092:   PetscFunctionListAdd(&STList,name,function);
1093:   PetscFunctionReturn(0);
1094: }