Actual source code: stfunc.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: The ST interface routines, callable by users
12: */
14: #include <slepc/private/stimpl.h>
16: PetscClassId ST_CLASSID = 0;
17: PetscLogEvent ST_SetUp = 0,ST_ComputeOperator = 0,ST_Apply = 0,ST_ApplyTranspose = 0,ST_ApplyHermitianTranspose = 0,ST_MatSetUp = 0,ST_MatMult = 0,ST_MatMultTranspose = 0,ST_MatSolve = 0,ST_MatSolveTranspose = 0;
18: static PetscBool STPackageInitialized = PETSC_FALSE;
20: const char *STMatModes[] = {"COPY","INPLACE","SHELL","STMatMode","ST_MATMODE_",NULL};
22: /*@C
23: STFinalizePackage - This function destroys everything in the SLEPc interface
24: to the `ST` package. It is called from `SlepcFinalize()`.
26: Level: developer
28: .seealso: [](ch:st), `SlepcFinalize()`, `STInitializePackage()`
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_slepc()` when using dynamic libraries, and
42: on the first call to `STCreate()` when using shared or static libraries.
44: Note:
45: This function never needs to be called by SLEPc users.
47: Level: developer
49: .seealso: [](ch:st), `ST`, `SlepcInitialize()`, `STFinalizePackage()`
50: @*/
51: PetscErrorCode STInitializePackage(void)
52: {
53: char logList[256];
54: PetscBool opt,pkg;
55: PetscClassId classids[1];
57: PetscFunctionBegin;
58: if (STPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
59: STPackageInitialized = PETSC_TRUE;
60: /* Register Classes */
61: PetscCall(PetscClassIdRegister("Spectral Transform",&ST_CLASSID));
62: /* Register Constructors */
63: PetscCall(STRegisterAll());
64: /* Register Events */
65: PetscCall(PetscLogEventRegister("STSetUp",ST_CLASSID,&ST_SetUp));
66: PetscCall(PetscLogEventRegister("STComputeOperatr",ST_CLASSID,&ST_ComputeOperator));
67: PetscCall(PetscLogEventRegister("STApply",ST_CLASSID,&ST_Apply));
68: PetscCall(PetscLogEventRegister("STApplyTranspose",ST_CLASSID,&ST_ApplyTranspose));
69: PetscCall(PetscLogEventRegister("STApplyHermTrans",ST_CLASSID,&ST_ApplyHermitianTranspose));
70: PetscCall(PetscLogEventRegister("STMatSetUp",ST_CLASSID,&ST_MatSetUp));
71: PetscCall(PetscLogEventRegister("STMatMult",ST_CLASSID,&ST_MatMult));
72: PetscCall(PetscLogEventRegister("STMatMultTranspose",ST_CLASSID,&ST_MatMultTranspose));
73: PetscCall(PetscLogEventRegister("STMatSolve",ST_CLASSID,&ST_MatSolve));
74: PetscCall(PetscLogEventRegister("STMatSolveTranspose",ST_CLASSID,&ST_MatSolveTranspose));
75: /* Process Info */
76: classids[0] = ST_CLASSID;
77: PetscCall(PetscInfoProcessClass("st",1,&classids[0]));
78: /* Process summary exclusions */
79: PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
80: if (opt) {
81: PetscCall(PetscStrInList("st",logList,',',&pkg));
82: if (pkg) PetscCall(PetscLogEventDeactivateClass(ST_CLASSID));
83: }
84: /* Register package finalizer */
85: PetscCall(PetscRegisterFinalize(STFinalizePackage));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /*@
90: STReset - Resets the `ST` context to the initial state (prior to setup)
91: and destroys any allocated `Vec`s and `Mat`s.
93: Collective
95: Input Parameter:
96: . st - the spectral transformation context
98: Level: advanced
100: .seealso: [](ch:st), `STDestroy()`
101: @*/
102: PetscErrorCode STReset(ST st)
103: {
104: PetscFunctionBegin;
106: if (!st) PetscFunctionReturn(PETSC_SUCCESS);
107: STCheckNotSeized(st,1);
108: PetscTryTypeMethod(st,reset);
109: if (st->ksp) PetscCall(KSPReset(st->ksp));
110: PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->T));
111: PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->A));
112: st->nmat = 0;
113: PetscCall(PetscFree(st->Astate));
114: PetscCall(MatDestroy(&st->Op));
115: PetscCall(MatDestroy(&st->P));
116: PetscCall(MatDestroy(&st->Pmat));
117: PetscCall(MatDestroyMatrices(st->nsplit,&st->Psplit));
118: st->nsplit = 0;
119: PetscCall(VecDestroyVecs(st->nwork,&st->work));
120: st->nwork = 0;
121: PetscCall(VecDestroy(&st->wb));
122: PetscCall(VecDestroy(&st->wht));
123: PetscCall(VecDestroy(&st->D));
124: st->state = ST_STATE_INITIAL;
125: st->opready = PETSC_FALSE;
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*@
130: STDestroy - Destroys ST context that was created with `STCreate()`.
132: Collective
134: Input Parameter:
135: . st - the spectral transformation context
137: Level: beginner
139: .seealso: [](ch:st), `STCreate()`, `STSetUp()`
140: @*/
141: PetscErrorCode STDestroy(ST *st)
142: {
143: PetscFunctionBegin;
144: if (!*st) PetscFunctionReturn(PETSC_SUCCESS);
146: if (--((PetscObject)*st)->refct > 0) { *st = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
147: PetscCall(STReset(*st));
148: PetscTryTypeMethod(*st,destroy);
149: PetscCall(KSPDestroy(&(*st)->ksp));
150: PetscCall(PetscHeaderDestroy(st));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /*@
155: STCreate - Creates a spectral transformation context.
157: Collective
159: Input Parameter:
160: . comm - MPI communicator
162: Output Parameter:
163: . newst - location to put the spectral transformation context
165: Level: beginner
167: .seealso: [](ch:st), `STSetUp()`, `STApply()`, `STDestroy()`, `ST`
168: @*/
169: PetscErrorCode STCreate(MPI_Comm comm,ST *newst)
170: {
171: ST st;
173: PetscFunctionBegin;
174: PetscAssertPointer(newst,2);
175: PetscCall(STInitializePackage());
176: PetscCall(SlepcHeaderCreate(st,ST_CLASSID,"ST","Spectral Transformation","ST",comm,STDestroy,STView));
178: st->A = NULL;
179: st->nmat = 0;
180: st->sigma = 0.0;
181: st->defsigma = 0.0;
182: st->matmode = ST_MATMODE_COPY;
183: st->str = UNKNOWN_NONZERO_PATTERN;
184: st->transform = PETSC_FALSE;
185: st->structured = PETSC_FALSE;
186: st->D = NULL;
187: st->Pmat = NULL;
188: st->Pmat_set = PETSC_FALSE;
189: st->Psplit = NULL;
190: st->nsplit = 0;
191: st->strp = UNKNOWN_NONZERO_PATTERN;
193: st->ksp = NULL;
194: st->usesksp = PETSC_FALSE;
195: st->nwork = 0;
196: st->work = NULL;
197: st->wb = NULL;
198: st->wht = NULL;
199: st->state = ST_STATE_INITIAL;
200: st->Astate = NULL;
201: st->T = NULL;
202: st->Op = NULL;
203: st->opseized = PETSC_FALSE;
204: st->opready = PETSC_FALSE;
205: st->P = NULL;
206: st->M = NULL;
207: st->sigma_set = PETSC_FALSE;
208: st->asymm = PETSC_FALSE;
209: st->aherm = PETSC_FALSE;
210: st->data = NULL;
212: *newst = st;
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /*
217: Checks whether the ST matrices are all symmetric or hermitian.
218: */
219: static inline PetscErrorCode STMatIsSymmetricKnown(ST st,PetscBool *symm,PetscBool *herm)
220: {
221: PetscInt i;
222: PetscBool sbaij=PETSC_FALSE,set,flg=PETSC_FALSE;
224: PetscFunctionBegin;
225: /* check if problem matrices are all sbaij */
226: for (i=0;i<st->nmat;i++) {
227: PetscCall(PetscObjectTypeCompareAny((PetscObject)st->A[i],&sbaij,MATSEQSBAIJ,MATMPISBAIJ,""));
228: if (!sbaij) break;
229: }
230: /* check if user has set the symmetric flag */
231: *symm = PETSC_TRUE;
232: for (i=0;i<st->nmat;i++) {
233: PetscCall(MatIsSymmetricKnown(st->A[i],&set,&flg));
234: if (!set || !flg) { *symm = PETSC_FALSE; break; }
235: }
236: if (sbaij) *symm = PETSC_TRUE;
237: #if defined(PETSC_USE_COMPLEX)
238: /* check if user has set the hermitian flag */
239: *herm = PETSC_TRUE;
240: for (i=0;i<st->nmat;i++) {
241: PetscCall(MatIsHermitianKnown(st->A[i],&set,&flg));
242: if (!set || !flg) { *herm = PETSC_FALSE; break; }
243: }
244: #else
245: *herm = *symm;
246: #endif
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: /*@
251: STSetMatrices - Sets the matrices associated with the eigenvalue problem.
253: Collective
255: Input Parameters:
256: + st - the spectral transformation context
257: . n - number of matrices in array `A`
258: - A - the array of matrices associated with the eigensystem
260: Notes:
261: It must be called before `STSetUp()`. If it is called again after `STSetUp()` then
262: the `ST` object is reset.
264: In standard eigenproblems only one matrix is passed, while in generalized
265: problems two matrices are provided. The number of matrices is larger in
266: polynomial eigenproblems.
268: In normal usage, matrices are provided via the corresponding `EPS` of `PEP`
269: interface function.
271: Level: intermediate
273: .seealso: [](ch:st), `STGetMatrix()`, `STGetNumMatrices()`, `STSetUp()`, `STReset()`
274: @*/
275: PetscErrorCode STSetMatrices(ST st,PetscInt n,Mat A[])
276: {
277: PetscInt i;
278: PetscBool same=PETSC_TRUE;
280: PetscFunctionBegin;
283: PetscCheck(n>0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more matrices, you have %" PetscInt_FMT,n);
284: PetscAssertPointer(A,3);
285: PetscCheckSameComm(st,1,*A,3);
286: STCheckNotSeized(st,1);
287: PetscCheck(!st->nsplit || st->nsplit==n,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"The number of matrices must be the same as in STSetSplitPreconditioner()");
289: if (st->state) {
290: if (n!=st->nmat) same = PETSC_FALSE;
291: for (i=0;same&&i<n;i++) {
292: if (A[i]!=st->A[i]) same = PETSC_FALSE;
293: }
294: if (!same) PetscCall(STReset(st));
295: } else same = PETSC_FALSE;
296: if (!same) {
297: PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->A));
298: PetscCall(PetscCalloc1(PetscMax(2,n),&st->A));
299: PetscCall(PetscFree(st->Astate));
300: PetscCall(PetscMalloc1(PetscMax(2,n),&st->Astate));
301: }
302: for (i=0;i<n;i++) {
304: PetscCall(PetscObjectReference((PetscObject)A[i]));
305: PetscCall(MatDestroy(&st->A[i]));
306: st->A[i] = A[i];
307: st->Astate[i] = ((PetscObject)A[i])->state;
308: }
309: if (n==1) {
310: st->A[1] = NULL;
311: st->Astate[1] = 0;
312: }
313: st->nmat = n;
314: if (same) st->state = ST_STATE_UPDATED;
315: else st->state = ST_STATE_INITIAL;
316: PetscCheck(!same || !st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Support for changing the matrices while using a split preconditioner is not implemented yet");
317: st->opready = PETSC_FALSE;
318: if (!same) PetscCall(STMatIsSymmetricKnown(st,&st->asymm,&st->aherm));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /*@
323: STGetMatrix - Gets the matrices associated with the original eigensystem.
325: Collective
327: Input Parameters:
328: + st - the spectral transformation context
329: - k - the index of the requested matrix (starting in 0)
331: Output Parameter:
332: . A - the requested matrix
334: Level: intermediate
336: .seealso: [](ch:st), `STSetMatrices()`, `STGetNumMatrices()`
337: @*/
338: PetscErrorCode STGetMatrix(ST st,PetscInt k,Mat *A)
339: {
340: PetscFunctionBegin;
343: PetscAssertPointer(A,3);
344: STCheckMatrices(st,1);
345: PetscCheck(k>=0 && k<st->nmat,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat-1);
346: PetscCheck(((PetscObject)st->A[k])->state==st->Astate[k],PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot retrieve original matrices (have been modified)");
347: *A = st->A[k];
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*@
352: STGetMatrixTransformed - Gets the matrices associated with the transformed eigensystem.
354: Not Collective
356: Input Parameters:
357: + st - the spectral transformation context
358: - k - the index of the requested matrix (starting in 0)
360: Output Parameter:
361: . T - the requested matrix
363: Level: developer
365: .seealso: [](ch:st), `STGetMatrix()`, `STGetNumMatrices()`
366: @*/
367: PetscErrorCode STGetMatrixTransformed(ST st,PetscInt k,Mat *T)
368: {
369: PetscFunctionBegin;
372: PetscAssertPointer(T,3);
373: STCheckMatrices(st,1);
374: PetscCheck(k>=0 && k<st->nmat,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat-1);
375: PetscCheck(st->T,PetscObjectComm((PetscObject)st),PETSC_ERR_POINTER,"There are no transformed matrices");
376: *T = st->T[k];
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*@
381: STGetNumMatrices - Returns the number of matrices stored in the `ST`.
383: Not Collective
385: Input Parameter:
386: . st - the spectral transformation context
388: Output Parameter:
389: . n - the number of matrices passed in `STSetMatrices()`
391: Level: intermediate
393: .seealso: [](ch:st), `STSetMatrices()`
394: @*/
395: PetscErrorCode STGetNumMatrices(ST st,PetscInt *n)
396: {
397: PetscFunctionBegin;
399: PetscAssertPointer(n,2);
400: *n = st->nmat;
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@
405: STResetMatrixState - Resets the stored state of the matrices in the `ST`.
407: Logically Collective
409: Input Parameter:
410: . st - the spectral transformation context
412: Note:
413: This is useful in solvers where the user matrices are modified during
414: the computation, as in nonlinear inverse iteration. The effect is that
415: `STGetMatrix()` will retrieve the modified matrices as if they were
416: the matrices originally provided by the user.
418: Level: developer
420: .seealso: [](ch:st), `STGetMatrix()`, `EPSPowerSetNonlinear()`
421: @*/
422: PetscErrorCode STResetMatrixState(ST st)
423: {
424: PetscInt i;
426: PetscFunctionBegin;
428: for (i=0;i<st->nmat;i++) st->Astate[i] = ((PetscObject)st->A[i])->state;
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
432: /*@
433: STSetPreconditionerMat - Sets the matrix to be used to build the preconditioner.
435: Collective
437: Input Parameters:
438: + st - the spectral transformation context
439: - mat - the matrix that will be used in constructing the preconditioner
441: Notes:
442: This matrix will be passed to the internal `KSP` object (via the last argument
443: of `KSPSetOperators()`) as the matrix to be used when constructing the preconditioner.
444: If no matrix is set or `mat` is set to `NULL`, then $A-\sigma B$ will be used
445: to build the preconditioner, being $\sigma$ the value set by `STSetShift()`.
447: More precisely, this is relevant for spectral transformations that represent
448: a rational matrix function, and use a `KSP` object for the denominator, called
449: $K$ in the description of `STGetOperator()`. It includes also the `STPRECOND` case.
450: If the user has a good approximation to matrix $K$ that can be used to build a
451: cheap preconditioner, it can be passed with this function. Note that it affects
452: only the `Pmat` argument of `KSPSetOperators()`, not the `Amat` argument.
454: If a preconditioner matrix is set, the default is to use an iterative `KSP`
455: rather than a direct method.
457: An alternative to pass an approximation of $A-\sigma B$ with this function is
458: to provide approximations of $A$ and $B$ via `STSetSplitPreconditioner()`. The
459: difference is that when $\sigma$ changes the preconditioner is recomputed.
461: Use `NULL` to remove a previously set matrix.
463: Level: advanced
465: .seealso: [](ch:st), `STGetPreconditionerMat()`, `STSetShift()`, `STGetOperator()`, `STSetSplitPreconditioner()`
466: @*/
467: PetscErrorCode STSetPreconditionerMat(ST st,Mat mat)
468: {
469: PetscFunctionBegin;
471: if (mat) {
473: PetscCheckSameComm(st,1,mat,2);
474: }
475: STCheckNotSeized(st,1);
476: PetscCheck(!mat || !st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot call both STSetPreconditionerMat and STSetSplitPreconditioner");
477: if (mat) PetscCall(PetscObjectReference((PetscObject)mat));
478: PetscCall(MatDestroy(&st->Pmat));
479: st->Pmat = mat;
480: st->Pmat_set = mat? PETSC_TRUE: PETSC_FALSE;
481: st->state = ST_STATE_INITIAL;
482: st->opready = PETSC_FALSE;
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: /*@
487: STGetPreconditionerMat - Returns the matrix previously set by `STSetPreconditionerMat()`.
489: Not Collective
491: Input Parameter:
492: . st - the spectral transformation context
494: Output Parameter:
495: . mat - the matrix that will be used in constructing the preconditioner or
496: `NULL` if no matrix was set by `STSetPreconditionerMat()`.
498: Level: advanced
500: .seealso: [](ch:st), `STSetPreconditionerMat()`
501: @*/
502: PetscErrorCode STGetPreconditionerMat(ST st,Mat *mat)
503: {
504: PetscFunctionBegin;
506: PetscAssertPointer(mat,2);
507: *mat = st->Pmat_set? st->Pmat: NULL;
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: /*@
512: STSetSplitPreconditioner - Sets the matrices from which to build the preconditioner
513: in split form.
515: Collective
517: Input Parameters:
518: + st - the spectral transformation context
519: . n - number of matrices
520: . Psplit - array of matrices
521: - strp - structure flag for `Psplit` matrices
523: Notes:
524: The number of matrices passed here must be the same as in `STSetMatrices()`.
526: For linear eigenproblems, the preconditioner matrix is computed as
527: $P(\sigma) = A_0-\sigma B_0$, where $A_0$ and $B_0$ are approximations of $A$ and $B$
528: (the eigenproblem matrices) provided via the `Psplit` array in this function.
529: Compared to `STSetPreconditionerMat()`, this function allows setting a preconditioner
530: in a way that is independent of the shift $\sigma$. Whenever the value of $\sigma$
531: changes the preconditioner is recomputed.
533: Similarly, for polynomial eigenproblems the matrix for the preconditioner
534: is expressed as $P(\sigma) = \sum_i P_i \phi_i(\sigma)$, for $i=1,\dots,n$, where
535: $P_i$ are given in `Psplit` and the $\phi_i$'s are the polynomial basis functions.
537: The structure flag provides information about the relative nonzero pattern of the
538: `Psplit` matrices, in the same way as in `STSetMatStructure()`.
540: Use `n=0` to reset a previously set split preconditioner.
542: Level: advanced
544: .seealso: [](ch:st), `STGetSplitPreconditionerTerm()`, `STGetSplitPreconditionerInfo()`, `STSetPreconditionerMat()`, `STSetMatrices()`, `STSetMatStructure()`
545: @*/
546: PetscErrorCode STSetSplitPreconditioner(ST st,PetscInt n,Mat Psplit[],MatStructure strp)
547: {
548: PetscInt i,N=0,M,M0=0,mloc,nloc,mloc0=0;
550: PetscFunctionBegin;
553: PetscCheck(n>=0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of n = %" PetscInt_FMT,n);
554: PetscCheck(!n || !st->Pmat_set,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot call both STSetPreconditionerMat and STSetSplitPreconditioner");
555: PetscCheck(!n || !st->nmat || st->nmat==n,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"The number of matrices must be the same as in STSetMatrices()");
556: if (n) PetscAssertPointer(Psplit,3);
558: STCheckNotSeized(st,1);
560: for (i=0;i<n;i++) {
562: PetscCheckSameComm(st,1,Psplit[i],3);
563: PetscCall(MatGetSize(Psplit[i],&M,&N));
564: PetscCall(MatGetLocalSize(Psplit[i],&mloc,&nloc));
565: 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);
566: 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);
567: if (!i) { M0 = M; mloc0 = mloc; }
568: 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);
569: 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);
570: PetscCall(PetscObjectReference((PetscObject)Psplit[i]));
571: }
573: if (st->Psplit) PetscCall(MatDestroyMatrices(st->nsplit,&st->Psplit));
575: /* allocate space and copy matrices */
576: if (n) {
577: PetscCall(PetscMalloc1(n,&st->Psplit));
578: for (i=0;i<n;i++) st->Psplit[i] = Psplit[i];
579: }
580: st->nsplit = n;
581: st->strp = strp;
582: st->state = ST_STATE_INITIAL;
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: /*@
587: STGetSplitPreconditionerTerm - Gets the matrices associated with
588: the split preconditioner.
590: Not Collective
592: Input Parameters:
593: + st - the spectral transformation context
594: - k - the index of the requested matrix (starting in 0)
596: Output Parameter:
597: . Psplit - the returned matrix
599: Level: advanced
601: .seealso: [](ch:st), `STSetSplitPreconditioner()`, `STGetSplitPreconditionerInfo()`
602: @*/
603: PetscErrorCode STGetSplitPreconditionerTerm(ST st,PetscInt k,Mat *Psplit)
604: {
605: PetscFunctionBegin;
608: PetscAssertPointer(Psplit,3);
609: PetscCheck(k>=0 && k<st->nsplit,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nsplit-1);
610: PetscCheck(st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"You have not called STSetSplitPreconditioner()");
611: *Psplit = st->Psplit[k];
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*@
616: STGetSplitPreconditionerInfo - Returns the number of matrices of the split
617: preconditioner, as well as the structure flag.
619: Not Collective
621: Input Parameter:
622: . st - the spectral transformation context
624: Output Parameters:
625: + n - the number of matrices passed in `STSetSplitPreconditioner()`
626: - strp - the matrix structure flag passed in `STSetSplitPreconditioner()`
628: Level: advanced
630: .seealso: [](ch:st), `STSetSplitPreconditioner()`, `STGetSplitPreconditionerTerm()`
631: @*/
632: PetscErrorCode STGetSplitPreconditionerInfo(ST st,PetscInt *n,MatStructure *strp)
633: {
634: PetscFunctionBegin;
636: if (n) *n = st->nsplit;
637: if (strp) *strp = st->strp;
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: /*@
642: STSetShift - Sets the shift associated with the spectral transformation.
644: Collective
646: Input Parameters:
647: + st - the spectral transformation context
648: - shift - the value of the shift
650: Notes:
651: In some spectral transformations, changing the shift may have associated
652: a lot of work, for example recomputing a factorization.
654: This function is normally not directly called by users, since the shift is
655: indirectly set by `EPSSetTarget()`.
657: Level: intermediate
659: .seealso: [](ch:st), `EPSSetTarget()`, `STGetShift()`, `STSetDefaultShift()`
660: @*/
661: PetscErrorCode STSetShift(ST st,PetscScalar shift)
662: {
663: PetscFunctionBegin;
667: if (st->sigma != shift) {
668: STCheckNotSeized(st,1);
669: if (st->state==ST_STATE_SETUP) PetscTryTypeMethod(st,setshift,shift);
670: st->sigma = shift;
671: }
672: st->sigma_set = PETSC_TRUE;
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: /*@
677: STGetShift - Gets the shift associated with the spectral transformation.
679: Not Collective
681: Input Parameter:
682: . st - the spectral transformation context
684: Output Parameter:
685: . shift - the value of the shift
687: Level: intermediate
689: .seealso: [](ch:st), `STSetShift()`
690: @*/
691: PetscErrorCode STGetShift(ST st,PetscScalar* shift)
692: {
693: PetscFunctionBegin;
695: PetscAssertPointer(shift,2);
696: *shift = st->sigma;
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: /*@
701: STSetDefaultShift - Sets the value of the shift that should be employed if
702: the user did not specify one.
704: Logically Collective
706: Input Parameters:
707: + st - the spectral transformation context
708: - defaultshift - the default value of the shift
710: Level: developer
712: .seealso: [](ch:st), `STSetShift()`
713: @*/
714: PetscErrorCode STSetDefaultShift(ST st,PetscScalar defaultshift)
715: {
716: PetscFunctionBegin;
719: if (st->defsigma != defaultshift) {
720: st->defsigma = defaultshift;
721: st->state = ST_STATE_INITIAL;
722: st->opready = PETSC_FALSE;
723: }
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: /*@
728: STScaleShift - Multiply the shift by a given factor.
730: Logically Collective
732: Input Parameters:
733: + st - the spectral transformation context
734: - factor - the scaling factor
736: Note:
737: This function does not update the transformation matrices, as opposed to
738: `STSetShift()`.
740: Level: developer
742: .seealso: [](ch:st), `STSetShift()`
743: @*/
744: PetscErrorCode STScaleShift(ST st,PetscScalar factor)
745: {
746: PetscFunctionBegin;
749: st->sigma *= factor;
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: /*@
754: STSetBalanceMatrix - Sets the diagonal matrix to be used for balancing.
756: Collective
758: Input Parameters:
759: + st - the spectral transformation context
760: - D - the diagonal matrix (represented as a vector)
762: Notes:
763: If this matrix is set, `STApply()` will effectively apply $D K^{-1} M D^{-1}$,
764: see discussion at `STGetOperator()`. Use `NULL` to reset a previously passed `D`.
766: Balancing is usually set via `EPSSetBalance()`, but the advanced user may use
767: this function to bypass the usual balancing methods.
769: Level: developer
771: .seealso: [](ch:st), `EPSSetBalance()`, `STApply()`, `STGetBalanceMatrix()`, `STGetOperator()`
772: @*/
773: PetscErrorCode STSetBalanceMatrix(ST st,Vec D)
774: {
775: PetscFunctionBegin;
777: if (st->D == D) PetscFunctionReturn(PETSC_SUCCESS);
778: STCheckNotSeized(st,1);
779: if (D) {
781: PetscCheckSameComm(st,1,D,2);
782: PetscCall(PetscObjectReference((PetscObject)D));
783: }
784: PetscCall(VecDestroy(&st->D));
785: st->D = D;
786: st->state = ST_STATE_INITIAL;
787: st->opready = PETSC_FALSE;
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: /*@
792: STGetBalanceMatrix - Gets the balance matrix used by the spectral transformation.
794: Not Collective
796: Input Parameter:
797: . st - the spectral transformation context
799: Output Parameter:
800: . D - the diagonal matrix (represented as a vector)
802: Note:
803: If the matrix was not set, a `NULL` pointer will be returned.
805: Level: developer
807: .seealso: [](ch:st), `STSetBalanceMatrix()`
808: @*/
809: PetscErrorCode STGetBalanceMatrix(ST st,Vec *D)
810: {
811: PetscFunctionBegin;
813: PetscAssertPointer(D,2);
814: *D = st->D;
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: /*@
819: STMatCreateVecs - Get vector(s) compatible with the `ST` matrices.
821: Collective
823: Input Parameter:
824: . st - the spectral transformation context
826: Output Parameters:
827: + right - (optional) vector that the matrix can be multiplied against
828: - left - (optional) vector that the matrix vector product can be stored in
830: Level: developer
832: .seealso: [](ch:st), `STMatCreateVecsEmpty()`
833: @*/
834: PetscErrorCode STMatCreateVecs(ST st,Vec *right,Vec *left)
835: {
836: PetscFunctionBegin;
837: STCheckMatrices(st,1);
838: PetscCall(MatCreateVecs(st->A[0],right,left));
839: PetscFunctionReturn(PETSC_SUCCESS);
840: }
842: /*@
843: STMatCreateVecsEmpty - Get vector(s) compatible with the `ST` matrices, i.e.,
844: with the same parallel layout, but without internal array.
846: Collective
848: Input Parameter:
849: . st - the spectral transformation context
851: Output Parameters:
852: + right - (optional) vector that the matrix can be multiplied against
853: - left - (optional) vector that the matrix vector product can be stored in
855: Level: developer
857: .seealso: [](ch:st), `STMatCreateVecs()`, `MatCreateVecsEmpty()`
858: @*/
859: PetscErrorCode STMatCreateVecsEmpty(ST st,Vec *right,Vec *left)
860: {
861: PetscFunctionBegin;
862: STCheckMatrices(st,1);
863: PetscCall(MatCreateVecsEmpty(st->A[0],right,left));
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: STMatGetSize - Returns the number of rows and columns of the `ST` matrices.
870: Not Collective
872: Input Parameter:
873: . st - the spectral transformation context
875: Output Parameters:
876: + m - the number of global rows
877: - n - the number of global columns
879: Level: developer
881: .seealso: [](ch:st), `STMatGetLocalSize()`
882: @*/
883: PetscErrorCode STMatGetSize(ST st,PetscInt *m,PetscInt *n)
884: {
885: PetscFunctionBegin;
886: STCheckMatrices(st,1);
887: PetscCall(MatGetSize(st->A[0],m,n));
888: PetscFunctionReturn(PETSC_SUCCESS);
889: }
891: /*@
892: STMatGetLocalSize - Returns the number of local rows and columns of the `ST` matrices.
894: Not Collective
896: Input Parameter:
897: . st - the spectral transformation context
899: Output Parameters:
900: + m - the number of local rows
901: - n - the number of local columns
903: Level: developer
905: .seealso: [](ch:st), `STMatGetSize()`
906: @*/
907: PetscErrorCode STMatGetLocalSize(ST st,PetscInt *m,PetscInt *n)
908: {
909: PetscFunctionBegin;
910: STCheckMatrices(st,1);
911: PetscCall(MatGetLocalSize(st->A[0],m,n));
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: /*@
916: STSetOptionsPrefix - Sets the prefix used for searching for all
917: `ST` options in the database.
919: Logically Collective
921: Input Parameters:
922: + st - the spectral transformation context
923: - prefix - the prefix string to prepend to all `ST` option requests
925: Notes:
926: A hyphen (-) must NOT be given at the beginning of the prefix name.
927: The first character of all runtime options is AUTOMATICALLY the
928: hyphen.
930: Level: advanced
932: .seealso: [](ch:st), `STAppendOptionsPrefix()`, `STGetOptionsPrefix()`
933: @*/
934: PetscErrorCode STSetOptionsPrefix(ST st,const char prefix[])
935: {
936: PetscFunctionBegin;
938: if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
939: PetscCall(KSPSetOptionsPrefix(st->ksp,prefix));
940: PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
941: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)st,prefix));
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: /*@
946: STAppendOptionsPrefix - Appends to the prefix used for searching for all
947: `ST` options in the database.
949: Logically Collective
951: Input Parameters:
952: + st - the spectral transformation context
953: - prefix - the prefix string to prepend to all `ST` option requests
955: Notes:
956: A hyphen (-) must NOT be given at the beginning of the prefix name.
957: The first character of all runtime options is AUTOMATICALLY the
958: hyphen.
960: Level: advanced
962: .seealso: [](ch:st), `STSetOptionsPrefix()`, `STGetOptionsPrefix()`
963: @*/
964: PetscErrorCode STAppendOptionsPrefix(ST st,const char prefix[])
965: {
966: PetscFunctionBegin;
968: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)st,prefix));
969: if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
970: PetscCall(KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix));
971: PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: /*@
976: STGetOptionsPrefix - Gets the prefix used for searching for all
977: `ST` options in the database.
979: Not Collective
981: Input Parameter:
982: . st - the spectral transformation context
984: Output Parameter:
985: . prefix - pointer to the prefix string used, is returned
987: Level: advanced
989: .seealso: [](ch:st), `STSetOptionsPrefix()`, `STAppendOptionsPrefix()`
990: @*/
991: PetscErrorCode STGetOptionsPrefix(ST st,const char *prefix[])
992: {
993: PetscFunctionBegin;
995: PetscAssertPointer(prefix,2);
996: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)st,prefix));
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /*@
1001: STView - Prints the `ST` data structure.
1003: Collective
1005: Input Parameters:
1006: + st - the ST context
1007: - viewer - optional visualization context
1009: Note:
1010: The available visualization contexts include
1011: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1012: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the
1013: first process opens the file; all other processes send their data to the
1014: first one to print
1016: The user can open an alternative visualization context with `PetscViewerASCIIOpen()`
1017: to output to a specified file.
1019: Use `STViewFromOptions()` to allow the user to select many different `PetscViewerType`
1020: and formats from the options database.
1022: Level: beginner
1024: .seealso: [](ch:st), `STCreate()`, `STViewFromOptions()`
1025: @*/
1026: PetscErrorCode STView(ST st,PetscViewer viewer)
1027: {
1028: STType cstr;
1029: char str[50];
1030: PetscBool isascii,isstring;
1032: PetscFunctionBegin;
1034: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)st),&viewer));
1036: PetscCheckSameComm(st,1,viewer,2);
1038: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1039: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring));
1040: if (isascii) {
1041: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)st,viewer));
1042: PetscCall(PetscViewerASCIIPushTab(viewer));
1043: PetscTryTypeMethod(st,view,viewer);
1044: PetscCall(PetscViewerASCIIPopTab(viewer));
1045: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),st->sigma,PETSC_FALSE));
1046: PetscCall(PetscViewerASCIIPrintf(viewer," shift: %s\n",str));
1047: PetscCall(PetscViewerASCIIPrintf(viewer," number of matrices: %" PetscInt_FMT "\n",st->nmat));
1048: switch (st->matmode) {
1049: case ST_MATMODE_COPY:
1050: break;
1051: case ST_MATMODE_INPLACE:
1052: PetscCall(PetscViewerASCIIPrintf(viewer," shifting the matrix and unshifting at exit\n"));
1053: break;
1054: case ST_MATMODE_SHELL:
1055: PetscCall(PetscViewerASCIIPrintf(viewer," using a shell matrix\n"));
1056: break;
1057: }
1058: if (st->nmat>1 && st->matmode != ST_MATMODE_SHELL) PetscCall(PetscViewerASCIIPrintf(viewer," nonzero pattern of the matrices: %s\n",MatStructures[st->str]));
1059: if (st->Psplit) PetscCall(PetscViewerASCIIPrintf(viewer," using split preconditioner matrices with %s\n",MatStructures[st->strp]));
1060: if (st->transform && st->nmat>2) PetscCall(PetscViewerASCIIPrintf(viewer," computing transformed matrices\n"));
1061: if (st->structured) PetscCall(PetscViewerASCIIPrintf(viewer," exploiting structure in the application of the operator\n"));
1062: } else if (isstring) {
1063: PetscCall(STGetType(st,&cstr));
1064: PetscCall(PetscViewerStringSPrintf(viewer," %-7.7s",cstr));
1065: PetscTryTypeMethod(st,view,viewer);
1066: }
1067: if (st->usesksp) {
1068: if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
1069: PetscCall(PetscViewerASCIIPushTab(viewer));
1070: PetscCall(KSPView(st->ksp,viewer));
1071: PetscCall(PetscViewerASCIIPopTab(viewer));
1072: }
1073: PetscFunctionReturn(PETSC_SUCCESS);
1074: }
1076: /*@
1077: STViewFromOptions - View (print) an `ST` object based on values in the options database
1079: Collective
1081: Input Parameters:
1082: + st - the spectral transformation context
1083: . obj - optional object that provides the options prefix used to query the options database
1084: - name - command line option
1086: Level: intermediate
1088: .seealso: [](ch:st), `STView()`, `STCreate()`, `PetscObjectViewFromOptions()`
1089: @*/
1090: PetscErrorCode STViewFromOptions(ST st,PetscObject obj,const char name[])
1091: {
1092: PetscFunctionBegin;
1094: PetscCall(PetscObjectViewFromOptions((PetscObject)st,obj,name));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: /*@C
1099: STRegister - Adds a method to the spectral transformation package.
1101: Not Collective
1103: Input Parameters:
1104: + name - name of a new user-defined transformation
1105: - function - routine to create method
1107: Notes:
1108: `STRegister()` may be called multiple times to add several user-defined
1109: spectral transformations.
1111: Example Usage:
1112: .vb
1113: STRegister("my_transform",MyTransformCreate);
1114: .ve
1116: Then, your spectral transform can be chosen with the procedural interface via
1117: .vb
1118: STSetType(st,"my_transform")
1119: .ve
1120: or at runtime via the option `-st_type my_transform`
1122: Level: advanced
1124: .seealso: [](ch:st), `STSetType()`, `STRegisterAll()`
1125: @*/
1126: PetscErrorCode STRegister(const char *name,PetscErrorCode (*function)(ST))
1127: {
1128: PetscFunctionBegin;
1129: PetscCall(STInitializePackage());
1130: PetscCall(PetscFunctionListAdd(&STList,name,function));
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: }