Actual source code: stfunc.c
slepc-3.20.2 2024-03-15
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: }