Actual source code: stfunc.c
slepc-3.21.0 2024-03-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: The ST interface routines, callable by users
12: */
14: #include <slepc/private/stimpl.h>
16: PetscClassId ST_CLASSID = 0;
17: PetscLogEvent ST_SetUp = 0,ST_ComputeOperator = 0,ST_Apply = 0,ST_ApplyTranspose = 0,ST_ApplyHermitianTranspose = 0,ST_MatSetUp = 0,ST_MatMult = 0,ST_MatMultTranspose = 0,ST_MatSolve = 0,ST_MatSolveTranspose = 0;
18: static PetscBool STPackageInitialized = PETSC_FALSE;
20: const char *STMatModes[] = {"COPY","INPLACE","SHELL","STMatMode","ST_MATMODE_",NULL};
22: /*@C
23: STFinalizePackage - This function destroys everything in the Slepc interface
24: to the ST package. It is called from SlepcFinalize().
26: Level: developer
28: .seealso: SlepcFinalize()
29: @*/
30: PetscErrorCode STFinalizePackage(void)
31: {
32: PetscFunctionBegin;
33: PetscCall(PetscFunctionListDestroy(&STList));
34: STPackageInitialized = PETSC_FALSE;
35: STRegisterAllCalled = PETSC_FALSE;
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: /*@C
40: STInitializePackage - This function initializes everything in the ST package.
41: It is called from PetscDLLibraryRegister() when using dynamic libraries, and
42: on the first call to STCreate() when using static libraries.
44: Level: developer
46: .seealso: SlepcInitialize()
47: @*/
48: PetscErrorCode STInitializePackage(void)
49: {
50: char logList[256];
51: PetscBool opt,pkg;
52: PetscClassId classids[1];
54: PetscFunctionBegin;
55: if (STPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
56: STPackageInitialized = PETSC_TRUE;
57: /* Register Classes */
58: PetscCall(PetscClassIdRegister("Spectral Transform",&ST_CLASSID));
59: /* Register Constructors */
60: PetscCall(STRegisterAll());
61: /* Register Events */
62: PetscCall(PetscLogEventRegister("STSetUp",ST_CLASSID,&ST_SetUp));
63: PetscCall(PetscLogEventRegister("STComputeOperatr",ST_CLASSID,&ST_ComputeOperator));
64: PetscCall(PetscLogEventRegister("STApply",ST_CLASSID,&ST_Apply));
65: PetscCall(PetscLogEventRegister("STApplyTranspose",ST_CLASSID,&ST_ApplyTranspose));
66: PetscCall(PetscLogEventRegister("STApplyHermTrans",ST_CLASSID,&ST_ApplyHermitianTranspose));
67: PetscCall(PetscLogEventRegister("STMatSetUp",ST_CLASSID,&ST_MatSetUp));
68: PetscCall(PetscLogEventRegister("STMatMult",ST_CLASSID,&ST_MatMult));
69: PetscCall(PetscLogEventRegister("STMatMultTranspose",ST_CLASSID,&ST_MatMultTranspose));
70: PetscCall(PetscLogEventRegister("STMatSolve",ST_CLASSID,&ST_MatSolve));
71: PetscCall(PetscLogEventRegister("STMatSolveTranspose",ST_CLASSID,&ST_MatSolveTranspose));
72: /* Process Info */
73: classids[0] = ST_CLASSID;
74: PetscCall(PetscInfoProcessClass("st",1,&classids[0]));
75: /* Process summary exclusions */
76: PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
77: if (opt) {
78: PetscCall(PetscStrInList("st",logList,',',&pkg));
79: if (pkg) PetscCall(PetscLogEventDeactivateClass(ST_CLASSID));
80: }
81: /* Register package finalizer */
82: PetscCall(PetscRegisterFinalize(STFinalizePackage));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: /*@
87: STReset - Resets the ST context to the initial state (prior to setup)
88: and destroys any allocated Vecs and Mats.
90: Collective
92: Input Parameter:
93: . st - the spectral transformation context
95: Level: advanced
97: .seealso: STDestroy()
98: @*/
99: PetscErrorCode STReset(ST st)
100: {
101: PetscFunctionBegin;
103: if (!st) PetscFunctionReturn(PETSC_SUCCESS);
104: STCheckNotSeized(st,1);
105: PetscTryTypeMethod(st,reset);
106: if (st->ksp) PetscCall(KSPReset(st->ksp));
107: PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->T));
108: PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->A));
109: st->nmat = 0;
110: PetscCall(PetscFree(st->Astate));
111: PetscCall(MatDestroy(&st->Op));
112: PetscCall(MatDestroy(&st->P));
113: PetscCall(MatDestroy(&st->Pmat));
114: PetscCall(MatDestroyMatrices(st->nsplit,&st->Psplit));
115: st->nsplit = 0;
116: PetscCall(VecDestroyVecs(st->nwork,&st->work));
117: st->nwork = 0;
118: PetscCall(VecDestroy(&st->wb));
119: PetscCall(VecDestroy(&st->wht));
120: PetscCall(VecDestroy(&st->D));
121: st->state = ST_STATE_INITIAL;
122: st->opready = PETSC_FALSE;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /*@C
127: STDestroy - Destroys ST context that was created with STCreate().
129: Collective
131: Input Parameter:
132: . st - the spectral transformation context
134: Level: beginner
136: .seealso: STCreate(), STSetUp()
137: @*/
138: PetscErrorCode STDestroy(ST *st)
139: {
140: PetscFunctionBegin;
141: if (!*st) PetscFunctionReturn(PETSC_SUCCESS);
143: if (--((PetscObject)*st)->refct > 0) { *st = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
144: PetscCall(STReset(*st));
145: PetscTryTypeMethod(*st,destroy);
146: PetscCall(KSPDestroy(&(*st)->ksp));
147: PetscCall(PetscHeaderDestroy(st));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: STCreate - Creates a spectral transformation context.
154: Collective
156: Input Parameter:
157: . comm - MPI communicator
159: Output Parameter:
160: . newst - location to put the spectral transformation context
162: Level: beginner
164: .seealso: STSetUp(), STApply(), STDestroy(), ST
165: @*/
166: PetscErrorCode STCreate(MPI_Comm comm,ST *newst)
167: {
168: ST st;
170: PetscFunctionBegin;
171: PetscAssertPointer(newst,2);
172: *newst = NULL;
173: PetscCall(STInitializePackage());
174: PetscCall(SlepcHeaderCreate(st,ST_CLASSID,"ST","Spectral Transformation","ST",comm,STDestroy,STView));
176: st->A = NULL;
177: st->nmat = 0;
178: st->sigma = 0.0;
179: st->defsigma = 0.0;
180: st->matmode = ST_MATMODE_COPY;
181: st->str = UNKNOWN_NONZERO_PATTERN;
182: st->transform = PETSC_FALSE;
183: st->D = NULL;
184: st->Pmat = NULL;
185: st->Pmat_set = PETSC_FALSE;
186: st->Psplit = NULL;
187: st->nsplit = 0;
188: st->strp = UNKNOWN_NONZERO_PATTERN;
190: st->ksp = NULL;
191: st->usesksp = PETSC_FALSE;
192: st->nwork = 0;
193: st->work = NULL;
194: st->wb = NULL;
195: st->wht = NULL;
196: st->state = ST_STATE_INITIAL;
197: st->Astate = NULL;
198: st->T = NULL;
199: st->Op = NULL;
200: st->opseized = PETSC_FALSE;
201: st->opready = PETSC_FALSE;
202: st->P = NULL;
203: st->M = NULL;
204: st->sigma_set = PETSC_FALSE;
205: st->asymm = PETSC_FALSE;
206: st->aherm = PETSC_FALSE;
207: st->data = NULL;
209: *newst = st;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*
214: Checks whether the ST matrices are all symmetric or hermitian.
215: */
216: static inline PetscErrorCode STMatIsSymmetricKnown(ST st,PetscBool *symm,PetscBool *herm)
217: {
218: PetscInt i;
219: PetscBool sbaij=PETSC_FALSE,set,flg=PETSC_FALSE;
221: PetscFunctionBegin;
222: /* check if problem matrices are all sbaij */
223: for (i=0;i<st->nmat;i++) {
224: PetscCall(PetscObjectTypeCompareAny((PetscObject)st->A[i],&sbaij,MATSEQSBAIJ,MATMPISBAIJ,""));
225: if (!sbaij) break;
226: }
227: /* check if user has set the symmetric flag */
228: *symm = PETSC_TRUE;
229: for (i=0;i<st->nmat;i++) {
230: PetscCall(MatIsSymmetricKnown(st->A[i],&set,&flg));
231: if (!set || !flg) { *symm = PETSC_FALSE; break; }
232: }
233: if (sbaij) *symm = PETSC_TRUE;
234: #if defined(PETSC_USE_COMPLEX)
235: /* check if user has set the hermitian flag */
236: *herm = PETSC_TRUE;
237: for (i=0;i<st->nmat;i++) {
238: PetscCall(MatIsHermitianKnown(st->A[i],&set,&flg));
239: if (!set || !flg) { *herm = PETSC_FALSE; break; }
240: }
241: #else
242: *herm = *symm;
243: #endif
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*@
248: STSetMatrices - Sets the matrices associated with the eigenvalue problem.
250: Collective
252: Input Parameters:
253: + st - the spectral transformation context
254: . n - number of matrices in array A
255: - A - the array of matrices associated with the eigensystem
257: Notes:
258: It must be called before STSetUp(). If it is called again after STSetUp() then
259: the ST object is reset.
261: Level: intermediate
263: .seealso: STGetMatrix(), STGetNumMatrices(), STSetUp(), STReset()
264: @*/
265: PetscErrorCode STSetMatrices(ST st,PetscInt n,Mat A[])
266: {
267: PetscInt i;
268: PetscBool same=PETSC_TRUE;
270: PetscFunctionBegin;
273: PetscCheck(n>0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more matrices, you have %" PetscInt_FMT,n);
274: PetscAssertPointer(A,3);
275: PetscCheckSameComm(st,1,*A,3);
276: STCheckNotSeized(st,1);
277: PetscCheck(!st->nsplit || st->nsplit==n,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"The number of matrices must be the same as in STSetSplitPreconditioner()");
279: if (st->state) {
280: if (n!=st->nmat) same = PETSC_FALSE;
281: for (i=0;same&&i<n;i++) {
282: if (A[i]!=st->A[i]) same = PETSC_FALSE;
283: }
284: if (!same) PetscCall(STReset(st));
285: } else same = PETSC_FALSE;
286: if (!same) {
287: PetscCall(MatDestroyMatrices(PetscMax(2,st->nmat),&st->A));
288: PetscCall(PetscCalloc1(PetscMax(2,n),&st->A));
289: PetscCall(PetscFree(st->Astate));
290: PetscCall(PetscMalloc1(PetscMax(2,n),&st->Astate));
291: }
292: for (i=0;i<n;i++) {
294: PetscCall(PetscObjectReference((PetscObject)A[i]));
295: PetscCall(MatDestroy(&st->A[i]));
296: st->A[i] = A[i];
297: st->Astate[i] = ((PetscObject)A[i])->state;
298: }
299: if (n==1) {
300: st->A[1] = NULL;
301: st->Astate[1] = 0;
302: }
303: st->nmat = n;
304: if (same) st->state = ST_STATE_UPDATED;
305: else st->state = ST_STATE_INITIAL;
306: PetscCheck(!same || !st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Support for changing the matrices while using a split preconditioner is not implemented yet");
307: st->opready = PETSC_FALSE;
308: if (!same) PetscCall(STMatIsSymmetricKnown(st,&st->asymm,&st->aherm));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@
313: STGetMatrix - Gets the matrices associated with the original eigensystem.
315: Not Collective
317: Input Parameters:
318: + st - the spectral transformation context
319: - k - the index of the requested matrix (starting in 0)
321: Output Parameters:
322: . A - the requested matrix
324: Level: intermediate
326: .seealso: STSetMatrices(), STGetNumMatrices()
327: @*/
328: PetscErrorCode STGetMatrix(ST st,PetscInt k,Mat *A)
329: {
330: PetscFunctionBegin;
333: PetscAssertPointer(A,3);
334: STCheckMatrices(st,1);
335: PetscCheck(k>=0 && k<st->nmat,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat-1);
336: PetscCheck(((PetscObject)st->A[k])->state==st->Astate[k],PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot retrieve original matrices (have been modified)");
337: *A = st->A[k];
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /*@
342: STGetMatrixTransformed - Gets the matrices associated with the transformed eigensystem.
344: Not Collective
346: Input Parameters:
347: + st - the spectral transformation context
348: - k - the index of the requested matrix (starting in 0)
350: Output Parameters:
351: . T - the requested matrix
353: Level: developer
355: .seealso: STGetMatrix(), STGetNumMatrices()
356: @*/
357: PetscErrorCode STGetMatrixTransformed(ST st,PetscInt k,Mat *T)
358: {
359: PetscFunctionBegin;
362: PetscAssertPointer(T,3);
363: STCheckMatrices(st,1);
364: PetscCheck(k>=0 && k<st->nmat,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat-1);
365: PetscCheck(st->T,PetscObjectComm((PetscObject)st),PETSC_ERR_POINTER,"There are no transformed matrices");
366: *T = st->T[k];
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: /*@
371: STGetNumMatrices - Returns the number of matrices stored in the ST.
373: Not Collective
375: Input Parameter:
376: . st - the spectral transformation context
378: Output Parameters:
379: . n - the number of matrices passed in STSetMatrices()
381: Level: intermediate
383: .seealso: STSetMatrices()
384: @*/
385: PetscErrorCode STGetNumMatrices(ST st,PetscInt *n)
386: {
387: PetscFunctionBegin;
389: PetscAssertPointer(n,2);
390: *n = st->nmat;
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*@
395: STResetMatrixState - Resets the stored state of the matrices in the ST.
397: Logically Collective
399: Input Parameter:
400: . st - the spectral transformation context
402: Note:
403: This is useful in solvers where the user matrices are modified during
404: the computation, as in nonlinear inverse iteration. The effect is that
405: STGetMatrix() will retrieve the modified matrices as if they were
406: the matrices originally provided by the user.
408: Level: developer
410: .seealso: STGetMatrix(), EPSPowerSetNonlinear()
411: @*/
412: PetscErrorCode STResetMatrixState(ST st)
413: {
414: PetscInt i;
416: PetscFunctionBegin;
418: for (i=0;i<st->nmat;i++) st->Astate[i] = ((PetscObject)st->A[i])->state;
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@
423: STSetPreconditionerMat - Sets the matrix to be used to build the preconditioner.
425: Collective
427: Input Parameters:
428: + st - the spectral transformation context
429: - mat - the matrix that will be used in constructing the preconditioner
431: Notes:
432: This matrix will be passed to the internal KSP object (via the last argument
433: of KSPSetOperators()) as the matrix to be used when constructing the preconditioner.
434: If no matrix is set or mat is set to NULL, A-sigma*B will be used
435: to build the preconditioner, being sigma the value set by STSetShift().
437: More precisely, this is relevant for spectral transformations that represent
438: a rational matrix function, and use a KSP object for the denominator, called
439: K in the description of STGetOperator(). It includes also the STPRECOND case.
440: If the user has a good approximation to matrix K that can be used to build a
441: cheap preconditioner, it can be passed with this function. Note that it affects
442: only the Pmat argument of KSPSetOperators(), not the Amat argument.
444: If a preconditioner matrix is set, the default is to use an iterative KSP
445: rather than a direct method.
447: An alternative to pass an approximation of A-sigma*B with this function is
448: to provide approximations of A and B via STSetSplitPreconditioner(). The
449: difference is that when sigma changes the preconditioner is recomputed.
451: Use NULL to remove a previously set matrix.
453: Level: advanced
455: .seealso: STGetPreconditionerMat(), STSetShift(), STGetOperator(), STSetSplitPreconditioner()
456: @*/
457: PetscErrorCode STSetPreconditionerMat(ST st,Mat mat)
458: {
459: PetscFunctionBegin;
461: if (mat) {
463: PetscCheckSameComm(st,1,mat,2);
464: }
465: STCheckNotSeized(st,1);
466: PetscCheck(!mat || !st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot call both STSetPreconditionerMat and STSetSplitPreconditioner");
467: if (mat) PetscCall(PetscObjectReference((PetscObject)mat));
468: PetscCall(MatDestroy(&st->Pmat));
469: st->Pmat = mat;
470: st->Pmat_set = mat? PETSC_TRUE: PETSC_FALSE;
471: st->state = ST_STATE_INITIAL;
472: st->opready = PETSC_FALSE;
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*@
477: STGetPreconditionerMat - Returns the matrix previously set by STSetPreconditionerMat().
479: Not Collective
481: Input Parameter:
482: . st - the spectral transformation context
484: Output Parameter:
485: . mat - the matrix that will be used in constructing the preconditioner or
486: NULL if no matrix was set by STSetPreconditionerMat().
488: Level: advanced
490: .seealso: STSetPreconditionerMat()
491: @*/
492: PetscErrorCode STGetPreconditionerMat(ST st,Mat *mat)
493: {
494: PetscFunctionBegin;
496: PetscAssertPointer(mat,2);
497: *mat = st->Pmat_set? st->Pmat: NULL;
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: /*@
502: STSetSplitPreconditioner - Sets the matrices from which to build the preconditioner
503: in split form.
505: Collective
507: Input Parameters:
508: + st - the spectral transformation context
509: . n - number of matrices
510: . Psplit - array of matrices
511: - strp - structure flag for Psplit matrices
513: Notes:
514: The number of matrices passed here must be the same as in STSetMatrices().
516: For linear eigenproblems, the preconditioner matrix is computed as
517: Pmat(sigma) = A0-sigma*B0, where A0 and B0 are approximations of A and B
518: (the eigenproblem matrices) provided via the Psplit array in this function.
519: Compared to STSetPreconditionerMat(), this function allows setting a preconditioner
520: in a way that is independent of the shift sigma. Whenever the value of sigma
521: changes the preconditioner is recomputed.
523: Similarly, for polynomial eigenproblems the matrix for the preconditioner
524: is expressed as Pmat(sigma) = sum_i Psplit_i*phi_i(sigma), for i=1,...,n, where
525: the phi_i's are the polynomial basis functions.
527: The structure flag provides information about the relative nonzero pattern of the
528: Psplit_i matrices, in the same way as in STSetMatStructure().
530: Use n=0 to reset a previously set split preconditioner.
532: Level: advanced
534: .seealso: STGetSplitPreconditionerTerm(), STGetSplitPreconditionerInfo(), STSetPreconditionerMat(), STSetMatrices(), STSetMatStructure()
535: @*/
536: PetscErrorCode STSetSplitPreconditioner(ST st,PetscInt n,Mat Psplit[],MatStructure strp)
537: {
538: PetscInt i,N=0,M,M0=0,mloc,nloc,mloc0=0;
540: PetscFunctionBegin;
543: PetscCheck(n>=0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of n = %" PetscInt_FMT,n);
544: PetscCheck(!n || !st->Pmat_set,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cannot call both STSetPreconditionerMat and STSetSplitPreconditioner");
545: PetscCheck(!n || !st->nmat || st->nmat==n,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"The number of matrices must be the same as in STSetMatrices()");
546: if (n) PetscAssertPointer(Psplit,3);
548: STCheckNotSeized(st,1);
550: for (i=0;i<n;i++) {
552: PetscCheckSameComm(st,1,Psplit[i],3);
553: PetscCall(MatGetSize(Psplit[i],&M,&N));
554: PetscCall(MatGetLocalSize(Psplit[i],&mloc,&nloc));
555: PetscCheck(M==N,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Psplit[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,M,N);
556: PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONG,"Psplit[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
557: if (!i) { M0 = M; mloc0 = mloc; }
558: PetscCheck(M==M0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_INCOMP,"Dimensions of Psplit[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,M,M0);
559: PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_INCOMP,"Local dimensions of Psplit[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
560: PetscCall(PetscObjectReference((PetscObject)Psplit[i]));
561: }
563: if (st->Psplit) PetscCall(MatDestroyMatrices(st->nsplit,&st->Psplit));
565: /* allocate space and copy matrices */
566: if (n) {
567: PetscCall(PetscMalloc1(n,&st->Psplit));
568: for (i=0;i<n;i++) st->Psplit[i] = Psplit[i];
569: }
570: st->nsplit = n;
571: st->strp = strp;
572: st->state = ST_STATE_INITIAL;
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*@
577: STGetSplitPreconditionerTerm - Gets the matrices associated with
578: the split preconditioner.
580: Not Collective
582: Input Parameters:
583: + st - the spectral transformation context
584: - k - the index of the requested matrix (starting in 0)
586: Output Parameter:
587: . Psplit - the returned matrix
589: Level: advanced
591: .seealso: STSetSplitPreconditioner(), STGetSplitPreconditionerInfo()
592: @*/
593: PetscErrorCode STGetSplitPreconditionerTerm(ST st,PetscInt k,Mat *Psplit)
594: {
595: PetscFunctionBegin;
598: PetscAssertPointer(Psplit,3);
599: PetscCheck(k>=0 && k<st->nsplit,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nsplit-1);
600: PetscCheck(st->Psplit,PetscObjectComm((PetscObject)st),PETSC_ERR_ORDER,"You have not called STSetSplitPreconditioner()");
601: *Psplit = st->Psplit[k];
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: /*@
606: STGetSplitPreconditionerInfo - Returns the number of matrices of the split
607: preconditioner, as well as the structure flag.
609: Not Collective
611: Input Parameter:
612: . st - the spectral transformation context
614: Output Parameters:
615: + n - the number of matrices passed in STSetSplitPreconditioner()
616: - strp - the matrix structure flag passed in STSetSplitPreconditioner()
618: Level: advanced
620: .seealso: STSetSplitPreconditioner(), STGetSplitPreconditionerTerm()
621: @*/
622: PetscErrorCode STGetSplitPreconditionerInfo(ST st,PetscInt *n,MatStructure *strp)
623: {
624: PetscFunctionBegin;
626: if (n) *n = st->nsplit;
627: if (strp) *strp = st->strp;
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: /*@
632: STSetShift - Sets the shift associated with the spectral transformation.
634: Collective
636: Input Parameters:
637: + st - the spectral transformation context
638: - shift - the value of the shift
640: Notes:
641: In some spectral transformations, changing the shift may have associated
642: a lot of work, for example recomputing a factorization.
644: This function is normally not directly called by users, since the shift is
645: indirectly set by EPSSetTarget().
647: Level: intermediate
649: .seealso: EPSSetTarget(), STGetShift(), STSetDefaultShift()
650: @*/
651: PetscErrorCode STSetShift(ST st,PetscScalar shift)
652: {
653: PetscFunctionBegin;
657: if (st->sigma != shift) {
658: STCheckNotSeized(st,1);
659: if (st->state==ST_STATE_SETUP) PetscTryTypeMethod(st,setshift,shift);
660: st->sigma = shift;
661: }
662: st->sigma_set = PETSC_TRUE;
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: /*@
667: STGetShift - Gets the shift associated with the spectral transformation.
669: Not Collective
671: Input Parameter:
672: . st - the spectral transformation context
674: Output Parameter:
675: . shift - the value of the shift
677: Level: intermediate
679: .seealso: STSetShift()
680: @*/
681: PetscErrorCode STGetShift(ST st,PetscScalar* shift)
682: {
683: PetscFunctionBegin;
685: PetscAssertPointer(shift,2);
686: *shift = st->sigma;
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: /*@
691: STSetDefaultShift - Sets the value of the shift that should be employed if
692: the user did not specify one.
694: Logically Collective
696: Input Parameters:
697: + st - the spectral transformation context
698: - defaultshift - the default value of the shift
700: Level: developer
702: .seealso: STSetShift()
703: @*/
704: PetscErrorCode STSetDefaultShift(ST st,PetscScalar defaultshift)
705: {
706: PetscFunctionBegin;
709: if (st->defsigma != defaultshift) {
710: st->defsigma = defaultshift;
711: st->state = ST_STATE_INITIAL;
712: st->opready = PETSC_FALSE;
713: }
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }
717: /*@
718: STScaleShift - Multiply the shift with a given factor.
720: Logically Collective
722: Input Parameters:
723: + st - the spectral transformation context
724: - factor - the scaling factor
726: Note:
727: This function does not update the transformation matrices, as opposed to
728: STSetShift().
730: Level: developer
732: .seealso: STSetShift()
733: @*/
734: PetscErrorCode STScaleShift(ST st,PetscScalar factor)
735: {
736: PetscFunctionBegin;
739: st->sigma *= factor;
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: STSetBalanceMatrix - Sets the diagonal matrix to be used for balancing.
746: Collective
748: Input Parameters:
749: + st - the spectral transformation context
750: - D - the diagonal matrix (represented as a vector)
752: Notes:
753: If this matrix is set, STApply will effectively apply D*OP*D^{-1}. Use NULL
754: to reset a previously passed D.
756: Balancing is usually set via EPSSetBalance, but the advanced user may use
757: this function to bypass the usual balancing methods.
759: Level: developer
761: .seealso: EPSSetBalance(), STApply(), STGetBalanceMatrix()
762: @*/
763: PetscErrorCode STSetBalanceMatrix(ST st,Vec D)
764: {
765: PetscFunctionBegin;
767: if (st->D == D) PetscFunctionReturn(PETSC_SUCCESS);
768: STCheckNotSeized(st,1);
769: if (D) {
771: PetscCheckSameComm(st,1,D,2);
772: PetscCall(PetscObjectReference((PetscObject)D));
773: }
774: PetscCall(VecDestroy(&st->D));
775: st->D = D;
776: st->state = ST_STATE_INITIAL;
777: st->opready = PETSC_FALSE;
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: /*@
782: STGetBalanceMatrix - Gets the balance matrix used by the spectral transformation.
784: Not Collective
786: Input Parameter:
787: . st - the spectral transformation context
789: Output Parameter:
790: . D - the diagonal matrix (represented as a vector)
792: Note:
793: If the matrix was not set, a null pointer will be returned.
795: Level: developer
797: .seealso: STSetBalanceMatrix()
798: @*/
799: PetscErrorCode STGetBalanceMatrix(ST st,Vec *D)
800: {
801: PetscFunctionBegin;
803: PetscAssertPointer(D,2);
804: *D = st->D;
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: /*@C
809: STMatCreateVecs - Get vector(s) compatible with the ST matrices.
811: Collective
813: Input Parameter:
814: . st - the spectral transformation context
816: Output Parameters:
817: + right - (optional) vector that the matrix can be multiplied against
818: - left - (optional) vector that the matrix vector product can be stored in
820: Level: developer
822: .seealso: STMatCreateVecsEmpty()
823: @*/
824: PetscErrorCode STMatCreateVecs(ST st,Vec *right,Vec *left)
825: {
826: PetscFunctionBegin;
827: STCheckMatrices(st,1);
828: PetscCall(MatCreateVecs(st->A[0],right,left));
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: /*@C
833: STMatCreateVecsEmpty - Get vector(s) compatible with the ST matrices, i.e. with the same
834: parallel layout, but without internal array.
836: Collective
838: Input Parameter:
839: . st - the spectral transformation context
841: Output Parameters:
842: + right - (optional) vector that the matrix can be multiplied against
843: - left - (optional) vector that the matrix vector product can be stored in
845: Level: developer
847: .seealso: STMatCreateVecs(), MatCreateVecsEmpty()
848: @*/
849: PetscErrorCode STMatCreateVecsEmpty(ST st,Vec *right,Vec *left)
850: {
851: PetscFunctionBegin;
852: STCheckMatrices(st,1);
853: PetscCall(MatCreateVecsEmpty(st->A[0],right,left));
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: /*@
858: STMatGetSize - Returns the number of rows and columns of the ST matrices.
860: Not Collective
862: Input Parameter:
863: . st - the spectral transformation context
865: Output Parameters:
866: + m - the number of global rows
867: - n - the number of global columns
869: Level: developer
871: .seealso: STMatGetLocalSize()
872: @*/
873: PetscErrorCode STMatGetSize(ST st,PetscInt *m,PetscInt *n)
874: {
875: PetscFunctionBegin;
876: STCheckMatrices(st,1);
877: PetscCall(MatGetSize(st->A[0],m,n));
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: /*@
882: STMatGetLocalSize - Returns the number of local rows and columns of the ST matrices.
884: Not Collective
886: Input Parameter:
887: . st - the spectral transformation context
889: Output Parameters:
890: + m - the number of local rows
891: - n - the number of local columns
893: Level: developer
895: .seealso: STMatGetSize()
896: @*/
897: PetscErrorCode STMatGetLocalSize(ST st,PetscInt *m,PetscInt *n)
898: {
899: PetscFunctionBegin;
900: STCheckMatrices(st,1);
901: PetscCall(MatGetLocalSize(st->A[0],m,n));
902: PetscFunctionReturn(PETSC_SUCCESS);
903: }
905: /*@C
906: STSetOptionsPrefix - Sets the prefix used for searching for all
907: ST options in the database.
909: Logically Collective
911: Input Parameters:
912: + st - the spectral transformation context
913: - prefix - the prefix string to prepend to all ST option requests
915: Notes:
916: A hyphen (-) must NOT be given at the beginning of the prefix name.
917: The first character of all runtime options is AUTOMATICALLY the
918: hyphen.
920: Level: advanced
922: .seealso: STAppendOptionsPrefix(), STGetOptionsPrefix()
923: @*/
924: PetscErrorCode STSetOptionsPrefix(ST st,const char *prefix)
925: {
926: PetscFunctionBegin;
928: if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
929: PetscCall(KSPSetOptionsPrefix(st->ksp,prefix));
930: PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
931: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)st,prefix));
932: PetscFunctionReturn(PETSC_SUCCESS);
933: }
935: /*@C
936: STAppendOptionsPrefix - Appends to the prefix used for searching for all
937: ST options in the database.
939: Logically Collective
941: Input Parameters:
942: + st - the spectral transformation context
943: - prefix - the prefix string to prepend to all ST option requests
945: Notes:
946: A hyphen (-) must NOT be given at the beginning of the prefix name.
947: The first character of all runtime options is AUTOMATICALLY the
948: hyphen.
950: Level: advanced
952: .seealso: STSetOptionsPrefix(), STGetOptionsPrefix()
953: @*/
954: PetscErrorCode STAppendOptionsPrefix(ST st,const char *prefix)
955: {
956: PetscFunctionBegin;
958: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)st,prefix));
959: if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
960: PetscCall(KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix));
961: PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
962: PetscFunctionReturn(PETSC_SUCCESS);
963: }
965: /*@C
966: STGetOptionsPrefix - Gets the prefix used for searching for all
967: ST options in the database.
969: Not Collective
971: Input Parameters:
972: . st - the spectral transformation context
974: Output Parameters:
975: . prefix - pointer to the prefix string used, is returned
977: Note:
978: On the Fortran side, the user should pass in a string 'prefix' of
979: sufficient length to hold the prefix.
981: Level: advanced
983: .seealso: STSetOptionsPrefix(), STAppendOptionsPrefix()
984: @*/
985: PetscErrorCode STGetOptionsPrefix(ST st,const char *prefix[])
986: {
987: PetscFunctionBegin;
989: PetscAssertPointer(prefix,2);
990: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)st,prefix));
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: /*@C
995: STView - Prints the ST data structure.
997: Collective
999: Input Parameters:
1000: + st - the ST context
1001: - viewer - optional visualization context
1003: Note:
1004: The available visualization contexts include
1005: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1006: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1007: output where only the first processor opens
1008: the file. All other processors send their
1009: data to the first processor to print.
1011: The user can open an alternative visualization contexts with
1012: PetscViewerASCIIOpen() (output to a specified file).
1014: Level: beginner
1016: .seealso: EPSView()
1017: @*/
1018: PetscErrorCode STView(ST st,PetscViewer viewer)
1019: {
1020: STType cstr;
1021: char str[50];
1022: PetscBool isascii,isstring;
1024: PetscFunctionBegin;
1026: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)st),&viewer));
1028: PetscCheckSameComm(st,1,viewer,2);
1030: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1031: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring));
1032: if (isascii) {
1033: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)st,viewer));
1034: PetscCall(PetscViewerASCIIPushTab(viewer));
1035: PetscTryTypeMethod(st,view,viewer);
1036: PetscCall(PetscViewerASCIIPopTab(viewer));
1037: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),st->sigma,PETSC_FALSE));
1038: PetscCall(PetscViewerASCIIPrintf(viewer," shift: %s\n",str));
1039: PetscCall(PetscViewerASCIIPrintf(viewer," number of matrices: %" PetscInt_FMT "\n",st->nmat));
1040: switch (st->matmode) {
1041: case ST_MATMODE_COPY:
1042: break;
1043: case ST_MATMODE_INPLACE:
1044: PetscCall(PetscViewerASCIIPrintf(viewer," shifting the matrix and unshifting at exit\n"));
1045: break;
1046: case ST_MATMODE_SHELL:
1047: PetscCall(PetscViewerASCIIPrintf(viewer," using a shell matrix\n"));
1048: break;
1049: }
1050: if (st->nmat>1 && st->matmode != ST_MATMODE_SHELL) PetscCall(PetscViewerASCIIPrintf(viewer," nonzero pattern of the matrices: %s\n",MatStructures[st->str]));
1051: if (st->Psplit) PetscCall(PetscViewerASCIIPrintf(viewer," using split preconditioner matrices with %s\n",MatStructures[st->strp]));
1052: if (st->transform && st->nmat>2) PetscCall(PetscViewerASCIIPrintf(viewer," computing transformed matrices\n"));
1053: } else if (isstring) {
1054: PetscCall(STGetType(st,&cstr));
1055: PetscCall(PetscViewerStringSPrintf(viewer," %-7.7s",cstr));
1056: PetscTryTypeMethod(st,view,viewer);
1057: }
1058: if (st->usesksp) {
1059: if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
1060: PetscCall(PetscViewerASCIIPushTab(viewer));
1061: PetscCall(KSPView(st->ksp,viewer));
1062: PetscCall(PetscViewerASCIIPopTab(viewer));
1063: }
1064: PetscFunctionReturn(PETSC_SUCCESS);
1065: }
1067: /*@C
1068: STViewFromOptions - View from options
1070: Collective
1072: Input Parameters:
1073: + st - the spectral transformation context
1074: . obj - optional object
1075: - name - command line option
1077: Level: intermediate
1079: .seealso: STView(), STCreate()
1080: @*/
1081: PetscErrorCode STViewFromOptions(ST st,PetscObject obj,const char name[])
1082: {
1083: PetscFunctionBegin;
1085: PetscCall(PetscObjectViewFromOptions((PetscObject)st,obj,name));
1086: PetscFunctionReturn(PETSC_SUCCESS);
1087: }
1089: /*@C
1090: STRegister - Adds a method to the spectral transformation package.
1092: Not Collective
1094: Input Parameters:
1095: + name - name of a new user-defined transformation
1096: - function - routine to create method context
1098: Notes:
1099: STRegister() may be called multiple times to add several user-defined
1100: spectral transformations.
1102: Example Usage:
1103: .vb
1104: STRegister("my_transform",MyTransformCreate);
1105: .ve
1107: Then, your spectral transform can be chosen with the procedural interface via
1108: $ STSetType(st,"my_transform")
1109: or at runtime via the option
1110: $ -st_type my_transform
1112: Level: advanced
1114: .seealso: STRegisterAll()
1115: @*/
1116: PetscErrorCode STRegister(const char *name,PetscErrorCode (*function)(ST))
1117: {
1118: PetscFunctionBegin;
1119: PetscCall(STInitializePackage());
1120: PetscCall(PetscFunctionListAdd(&STList,name,function));
1121: PetscFunctionReturn(PETSC_SUCCESS);
1122: }