1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: ST interface routines, callable by users
12: */
14: #include <slepc/private/stimpl.h> 16: PetscErrorCode STApply_Generic(ST st,Vec x,Vec y) 17: {
21: if (st->M && st->P) {
22: MatMult(st->M,x,st->work[0]);
23: STMatSolve(st,st->work[0],y);
24: } else if (st->M) {
25: MatMult(st->M,x,y);
26: } else {
27: STMatSolve(st,x,y);
28: }
29: return(0);
30: }
32: /*@
33: STApply - Applies the spectral transformation operator to a vector, for
34: instance (A - sB)^-1 B in the case of the shift-and-invert transformation
35: and generalized eigenproblem.
37: Collective on st
39: Input Parameters:
40: + st - the spectral transformation context
41: - x - input vector
43: Output Parameter:
44: . y - output vector
46: Level: developer
48: .seealso: STApplyTranspose(), STApplyHermitianTranspose()
49: @*/
50: PetscErrorCode STApply(ST st,Vec x,Vec y) 51: {
53: Mat Op;
60: STCheckMatrices(st,1);
61: if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
62: VecSetErrorIfLocked(y,3);
63: if (!st->ops->apply) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST does not have apply");
64: STGetOperator_Private(st,&Op);
65: MatMult(Op,x,y);
66: return(0);
67: }
69: PetscErrorCode STApplyMat_Generic(ST st,Mat B,Mat C) 70: {
72: Mat work;
75: if (st->M && st->P) {
76: MatMatMult(st->M,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&work);
77: STMatMatSolve(st,work,C);
78: MatDestroy(&work);
79: } else if (st->M) {
80: MatMatMult(st->M,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
81: } else {
82: STMatMatSolve(st,B,C);
83: }
84: return(0);
85: }
87: /*@
88: STApplyMat - Applies the spectral transformation operator to a matrix, for
89: instance (A - sB)^-1 B in the case of the shift-and-invert transformation
90: and generalized eigenproblem.
92: Collective on st
94: Input Parameters:
95: + st - the spectral transformation context
96: - X - input matrix
98: Output Parameter:
99: . y - output matrix
101: Level: developer
103: .seealso: STApply()
104: @*/
105: PetscErrorCode STApplyMat(ST st,Mat X,Mat Y)106: {
114: STCheckMatrices(st,1);
115: if (X == Y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"X and Y must be different matrices");
116: if (!st->ops->applymat) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST does not have applymat");
117: (*st->ops->applymat)(st,X,Y);
118: return(0);
119: }
121: PetscErrorCode STApplyTranspose_Generic(ST st,Vec x,Vec y)122: {
126: if (st->M && st->P) {
127: STMatSolveTranspose(st,x,st->work[0]);
128: MatMultTranspose(st->M,st->work[0],y);
129: } else if (st->M) {
130: MatMultTranspose(st->M,x,y);
131: } else {
132: STMatSolveTranspose(st,x,y);
133: }
134: return(0);
135: }
137: /*@
138: STApplyTranspose - Applies the transpose of the operator to a vector, for
139: instance B^T(A - sB)^-T in the case of the shift-and-invert transformation
140: and generalized eigenproblem.
142: Collective on st
144: Input Parameters:
145: + st - the spectral transformation context
146: - x - input vector
148: Output Parameter:
149: . y - output vector
151: Level: developer
153: .seealso: STApply(), STApplyHermitianTranspose()
154: @*/
155: PetscErrorCode STApplyTranspose(ST st,Vec x,Vec y)156: {
158: Mat Op;
165: STCheckMatrices(st,1);
166: if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
167: VecSetErrorIfLocked(y,3);
168: if (!st->ops->applytrans) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST does not have applytrans");
169: STGetOperator_Private(st,&Op);
170: MatMultTranspose(Op,x,y);
171: return(0);
172: }
174: /*@
175: STApplyHermitianTranspose - Applies the hermitian-transpose of the operator
176: to a vector, for instance B^H(A - sB)^-H in the case of the shift-and-invert
177: transformation and generalized eigenproblem.
179: Collective on st
181: Input Parameters:
182: + st - the spectral transformation context
183: - x - input vector
185: Output Parameter:
186: . y - output vector
188: Note:
189: Currently implemented via STApplyTranspose() with appropriate conjugation.
191: Level: developer
193: .seealso: STApply(), STApplyTranspose()
194: @*/
195: PetscErrorCode STApplyHermitianTranspose(ST st,Vec x,Vec y)196: {
198: Mat Op;
205: STCheckMatrices(st,1);
206: if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
207: VecSetErrorIfLocked(y,3);
208: if (!st->ops->applytrans) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST does not have applytrans");
209: STGetOperator_Private(st,&Op);
210: MatMultHermitianTranspose(Op,x,y);
211: return(0);
212: }
214: /*@
215: STGetBilinearForm - Returns the matrix used in the bilinear form with a
216: generalized problem with semi-definite B.
218: Not collective, though a parallel Mat may be returned
220: Input Parameters:
221: . st - the spectral transformation context
223: Output Parameter:
224: . B - output matrix
226: Notes:
227: The output matrix B must be destroyed after use. It will be NULL in
228: case of standard eigenproblems.
230: Level: developer
231: @*/
232: PetscErrorCode STGetBilinearForm(ST st,Mat *B)233: {
240: STCheckMatrices(st,1);
241: (*st->ops->getbilinearform)(st,B);
242: return(0);
243: }
245: PetscErrorCode STGetBilinearForm_Default(ST st,Mat *B)246: {
250: if (st->nmat==1) *B = NULL;
251: else {
252: *B = st->A[1];
253: PetscObjectReference((PetscObject)*B);
254: }
255: return(0);
256: }
258: static PetscErrorCode MatMult_STOperator(Mat Op,Vec x,Vec y)259: {
261: ST st;
264: MatShellGetContext(Op,(void**)&st);
265: STSetUp(st);
266: PetscLogEventBegin(ST_Apply,st,x,y,0);
267: if (st->D) { /* with balancing */
268: VecPointwiseDivide(st->wb,x,st->D);
269: (*st->ops->apply)(st,st->wb,y);
270: VecPointwiseMult(y,y,st->D);
271: } else {
272: (*st->ops->apply)(st,x,y);
273: }
274: PetscLogEventEnd(ST_Apply,st,x,y,0);
275: return(0);
276: }
278: static PetscErrorCode MatMultTranspose_STOperator(Mat Op,Vec x,Vec y)279: {
281: ST st;
284: MatShellGetContext(Op,(void**)&st);
285: STSetUp(st);
286: PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);
287: if (st->D) { /* with balancing */
288: VecPointwiseMult(st->wb,x,st->D);
289: (*st->ops->applytrans)(st,st->wb,y);
290: VecPointwiseDivide(y,y,st->D);
291: } else {
292: (*st->ops->applytrans)(st,x,y);
293: }
294: PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);
295: return(0);
296: }
298: #if defined(PETSC_USE_COMPLEX)
299: static PetscErrorCode MatMultHermitianTranspose_STOperator(Mat Op,Vec x,Vec y)300: {
302: ST st;
305: MatShellGetContext(Op,(void**)&st);
306: STSetUp(st);
307: PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);
308: if (!st->wht) {
309: MatCreateVecs(st->A[0],&st->wht,NULL);
310: PetscLogObjectParent((PetscObject)st,(PetscObject)st->wht);
311: }
312: VecCopy(x,st->wht);
313: VecConjugate(st->wht);
314: if (st->D) { /* with balancing */
315: VecPointwiseMult(st->wb,st->wht,st->D);
316: (*st->ops->applytrans)(st,st->wb,y);
317: VecPointwiseDivide(y,y,st->D);
318: } else {
319: (*st->ops->applytrans)(st,st->wht,y);
320: }
321: VecConjugate(y);
322: PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);
323: return(0);
324: }
325: #endif
327: static PetscErrorCode MatMatMult_STOperator(Mat Op,Mat B,Mat C,void *ctx)328: {
330: ST st;
333: MatShellGetContext(Op,(void**)&st);
334: STSetUp(st);
335: PetscLogEventBegin(ST_Apply,st,B,C,0);
336: STApplyMat_Generic(st,B,C);
337: PetscLogEventEnd(ST_Apply,st,B,C,0);
338: return(0);
339: }
341: PetscErrorCode STGetOperator_Private(ST st,Mat *Op)342: {
344: PetscInt m,n,M,N;
345: Vec v;
346: VecType vtype;
349: if (!st->Op) {
350: if (Op) *Op = NULL;
351: /* create the shell matrix */
352: MatGetLocalSize(st->A[0],&m,&n);
353: MatGetSize(st->A[0],&M,&N);
354: MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,st,&st->Op);
355: MatShellSetOperation(st->Op,MATOP_MULT,(void(*)(void))MatMult_STOperator);
356: MatShellSetOperation(st->Op,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator);
357: #if defined(PETSC_USE_COMPLEX)
358: MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_STOperator);
359: #else
360: MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator);
361: #endif
362: if (!st->D && st->ops->apply==STApply_Generic) {
363: MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSE,MATDENSE);
364: }
365: /* make sure the shell matrix generates a vector of the same type as the problem matrices */
366: MatCreateVecs(st->A[0],&v,NULL);
367: VecGetType(v,&vtype);
368: MatShellSetVecType(st->Op,vtype);
369: VecDestroy(&v);
370: /* build the operator matrices */
371: STComputeOperator(st);
372: }
373: if (Op) *Op = st->Op;
374: return(0);
375: }
377: /*@
378: STGetOperator - Returns a shell matrix that represents the operator of the
379: spectral transformation.
381: Collective on st
383: Input Parameter:
384: . st - the spectral transformation context
386: Output Parameter:
387: . Op - operator matrix
389: Notes:
390: The operator is defined in linear eigenproblems only, not in polynomial ones,
391: so the call will fail if more than 2 matrices were passed in STSetMatrices().
393: The returned shell matrix is essentially a wrapper to the STApply() and
394: STApplyTranspose() operations. The operator can often be expressed as
396: .vb
397: Op = D*inv(K)*M*inv(D)
398: .ve
400: where D is the balancing matrix, and M and K are two matrices corresponding
401: to the numerator and denominator for spectral transformations that represent
402: a rational matrix function. In the case of STSHELL, the inner part inv(K)*M
403: is replaced by the user-provided operation from STShellSetApply().
405: The preconditioner matrix K typically depends on the value of the shift, and
406: its inverse is handled via an internal KSP object. Normal usage does not
407: require explicitly calling STGetOperator(), but it can be used to force the
408: creation of K and M, and then K is passed to the KSP. This is useful for
409: setting options associated with the PCFactor (to set MUMPS options, for instance).
411: The returned matrix must NOT be destroyed by the user. Instead, when no
412: longer needed it must be returned with STRestoreOperator(). In particular,
413: this is required before modifying the ST matrices or the shift.
415: A NULL pointer can be passed in Op in case the matrix is not required but we
416: want to force its creation. In this case, STRestoreOperator() should not be
417: called.
419: Level: advanced
421: .seealso: STApply(), STApplyTranspose(), STSetBalanceMatrix(), STShellSetApply(),
422: STGetKSP(), STSetShift(), STRestoreOperator(), STSetMatrices()
423: @*/
424: PetscErrorCode STGetOperator(ST st,Mat *Op)425: {
431: STCheckMatrices(st,1);
432: STCheckNotSeized(st,1);
433: if (st->nmat>2) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"The operator is not defined in polynomial eigenproblems");
434: STGetOperator_Private(st,Op);
435: if (Op) st->opseized = PETSC_TRUE;
436: return(0);
437: }
439: /*@
440: STRestoreOperator - Restore the previously seized operator matrix.
442: Collective on st
444: Input Parameters:
445: + st - the spectral transformation context
446: - Op - operator matrix
448: Notes:
449: The arguments must match the corresponding call to STGetOperator().
451: Level: advanced
453: .seealso: STGetOperator()
454: @*/
455: PetscErrorCode STRestoreOperator(ST st,Mat *Op)456: {
461: if (!st->opseized) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must be called after STGetOperator()");
462: *Op = NULL;
463: st->opseized = PETSC_FALSE;
464: return(0);
465: }
467: /*
468: STComputeOperator - Computes the matrices that constitute the operator
470: Op = D*inv(K)*M*inv(D).
472: K and M are computed here (D is user-provided) from the system matrices
473: and the shift sigma (whenever these are changed, this function recomputes
474: K and M). This is used only in linear eigenproblems (nmat<3).
476: K is the "preconditioner matrix": it is the denominator in rational operators,
477: e.g. (A-sigma*B) in shift-and-invert. In non-rational transformations such
478: as STFILTER, K=NULL which means identity. After computing K, it is passed to
479: the internal KSP object via KSPSetOperators.
481: M is the numerator in rational operators. If unused it is set to NULL (e.g.
482: in STPRECOND).
484: STSHELL does not compute anything here, but sets the flag as if it was ready.
485: */
486: PetscErrorCode STComputeOperator(ST st)487: {
489: PC pc;
494: if (!st->opready && st->ops->computeoperator) {
495: PetscInfo(st,"Building the operator matrices\n");
496: STCheckMatrices(st,1);
497: if (!st->T) {
498: PetscCalloc1(PetscMax(2,st->nmat),&st->T);
499: PetscLogObjectMemory((PetscObject)st,PetscMax(2,st->nmat)*sizeof(Mat));
500: }
501: PetscLogEventBegin(ST_ComputeOperator,st,0,0,0);
502: (*st->ops->computeoperator)(st);
503: PetscLogEventEnd(ST_ComputeOperator,st,0,0,0);
504: if (st->usesksp) {
505: if (!st->ksp) { STGetKSP(st,&st->ksp); }
506: if (st->P) {
507: STSetDefaultKSP(st);
508: STCheckFactorPackage(st);
509: KSPSetOperators(st->ksp,st->P,st->P);
510: } else {
511: /* STPRECOND defaults to PCNONE if st->P is empty */
512: KSPGetPC(st->ksp,&pc);
513: PCSetType(pc,PCNONE);
514: }
515: }
516: }
517: st->opready = PETSC_TRUE;
518: return(0);
519: }
521: /*@
522: STSetUp - Prepares for the use of a spectral transformation.
524: Collective on st
526: Input Parameter:
527: . st - the spectral transformation context
529: Level: advanced
531: .seealso: STCreate(), STApply(), STDestroy()
532: @*/
533: PetscErrorCode STSetUp(ST st)534: {
535: PetscInt i,n,k;
541: STCheckMatrices(st,1);
542: switch (st->state) {
543: case ST_STATE_INITIAL:
544: PetscInfo(st,"Setting up new ST\n");
545: if (!((PetscObject)st)->type_name) {
546: STSetType(st,STSHIFT);
547: }
548: break;
549: case ST_STATE_SETUP:
550: return(0);
551: case ST_STATE_UPDATED:
552: PetscInfo(st,"Setting up updated ST\n");
553: break;
554: }
555: PetscLogEventBegin(ST_SetUp,st,0,0,0);
556: if (st->state!=ST_STATE_UPDATED) {
557: if (!(st->nmat<3 && st->opready)) {
558: if (st->T) {
559: for (i=0;i<PetscMax(2,st->nmat);i++) {
560: MatDestroy(&st->T[i]);
561: }
562: }
563: MatDestroy(&st->P);
564: }
565: }
566: if (st->D) {
567: MatGetLocalSize(st->A[0],NULL,&n);
568: VecGetLocalSize(st->D,&k);
569: if (n != k) SETERRQ2(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Balance matrix has wrong dimension %D (should be %D)",k,n);
570: if (!st->wb) {
571: VecDuplicate(st->D,&st->wb);
572: PetscLogObjectParent((PetscObject)st,(PetscObject)st->wb);
573: }
574: }
575: if (st->nmat<3 && st->transform) {
576: STComputeOperator(st);
577: } else {
578: if (!st->T) {
579: PetscCalloc1(PetscMax(2,st->nmat),&st->T);
580: PetscLogObjectMemory((PetscObject)st,PetscMax(2,st->nmat)*sizeof(Mat));
581: }
582: }
583: if (st->ops->setup) { (*st->ops->setup)(st); }
584: st->state = ST_STATE_SETUP;
585: PetscLogEventEnd(ST_SetUp,st,0,0,0);
586: return(0);
587: }
589: /*
590: Computes coefficients for the transformed polynomial,
591: and stores the result in argument S.
593: alpha - value of the parameter of the transformed polynomial
594: beta - value of the previous shift (only used in inplace mode)
595: k - index of first matrix included in the computation
596: coeffs - coefficients of the expansion
597: initial - true if this is the first time (only relevant for shell mode)
598: */
599: PetscErrorCode STMatMAXPY_Private(ST st,PetscScalar alpha,PetscScalar beta,PetscInt k,PetscScalar *coeffs,PetscBool initial,Mat *S)600: {
602: PetscInt *matIdx=NULL,nmat,i,ini=-1;
603: PetscScalar t=1.0,ta,gamma;
604: PetscBool nz=PETSC_FALSE;
607: nmat = st->nmat-k;
608: switch (st->matmode) {
609: case ST_MATMODE_INPLACE:
610: if (st->nmat>2) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for polynomial eigenproblems");
611: if (initial) {
612: PetscObjectReference((PetscObject)st->A[0]);
613: *S = st->A[0];
614: gamma = alpha;
615: } else gamma = alpha-beta;
616: if (gamma != 0.0) {
617: if (st->nmat>1) {
618: MatAXPY(*S,gamma,st->A[1],st->str);
619: } else {
620: MatShift(*S,gamma);
621: }
622: }
623: break;
624: case ST_MATMODE_SHELL:
625: if (initial) {
626: if (st->nmat>2) {
627: PetscMalloc1(nmat,&matIdx);
628: for (i=0;i<nmat;i++) matIdx[i] = k+i;
629: }
630: STMatShellCreate(st,alpha,nmat,matIdx,coeffs,S);
631: PetscLogObjectParent((PetscObject)st,(PetscObject)*S);
632: if (st->nmat>2) { PetscFree(matIdx); }
633: } else {
634: STMatShellShift(*S,alpha);
635: }
636: break;
637: case ST_MATMODE_COPY:
638: if (coeffs) {
639: for (i=0;i<nmat && ini==-1;i++) {
640: if (coeffs[i]!=0.0) ini = i;
641: else t *= alpha;
642: }
643: if (coeffs[ini] != 1.0) nz = PETSC_TRUE;
644: for (i=ini+1;i<nmat&&!nz;i++) if (coeffs[i]!=0.0) nz = PETSC_TRUE;
645: } else { nz = PETSC_TRUE; ini = 0; }
646: if ((alpha == 0.0 || !nz) && t==1.0) {
647: PetscObjectReference((PetscObject)st->A[k+ini]);
648: MatDestroy(S);
649: *S = st->A[k+ini];
650: } else {
651: if (*S && *S!=st->A[k+ini]) {
652: MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
653: MatCopy(st->A[k+ini],*S,DIFFERENT_NONZERO_PATTERN);
654: } else {
655: MatDestroy(S);
656: MatDuplicate(st->A[k+ini],MAT_COPY_VALUES,S);
657: MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
658: PetscLogObjectParent((PetscObject)st,(PetscObject)*S);
659: }
660: if (coeffs && coeffs[ini]!=1.0) {
661: MatScale(*S,coeffs[ini]);
662: }
663: for (i=ini+k+1;i<PetscMax(2,st->nmat);i++) {
664: t *= alpha;
665: ta = t;
666: if (coeffs) ta *= coeffs[i-k];
667: if (ta!=0.0) {
668: if (st->nmat>1) {
669: MatAXPY(*S,ta,st->A[i],st->str);
670: } else {
671: MatShift(*S,ta);
672: }
673: }
674: }
675: }
676: }
677: MatSetOption(*S,MAT_SYMMETRIC,st->asymm);
678: MatSetOption(*S,MAT_HERMITIAN,(PetscImaginaryPart(st->sigma)==0.0)?st->aherm:PETSC_FALSE);
679: return(0);
680: }
682: /*
683: Computes the values of the coefficients required by STMatMAXPY_Private
684: for the case of monomial basis.
685: */
686: PetscErrorCode STCoeffs_Monomial(ST st, PetscScalar *coeffs)687: {
688: PetscInt k,i,ini,inip;
691: /* Compute binomial coefficients */
692: ini = (st->nmat*(st->nmat-1))/2;
693: for (i=0;i<st->nmat;i++) coeffs[ini+i]=1.0;
694: for (k=st->nmat-1;k>=1;k--) {
695: inip = ini+1;
696: ini = (k*(k-1))/2;
697: coeffs[ini] = 1.0;
698: for (i=1;i<k;i++) coeffs[ini+i] = coeffs[ini+i-1]+coeffs[inip+i-1];
699: }
700: return(0);
701: }
703: /*@
704: STPostSolve - Optional post-solve phase, intended for any actions that must
705: be performed on the ST object after the eigensolver has finished.
707: Collective on st
709: Input Parameters:
710: . st - the spectral transformation context
712: Level: developer
714: .seealso: EPSSolve()
715: @*/
716: PetscErrorCode STPostSolve(ST st)717: {
723: if (st->ops->postsolve) {
724: (*st->ops->postsolve)(st);
725: }
726: return(0);
727: }
729: /*@
730: STBackTransform - Back-transformation phase, intended for
731: spectral transformations which require to transform the computed
732: eigenvalues back to the original eigenvalue problem.
734: Not Collective
736: Input Parameters:
737: + st - the spectral transformation context
738: eigr - real part of a computed eigenvalue
739: - eigi - imaginary part of a computed eigenvalue
741: Level: developer
743: .seealso: STIsInjective()
744: @*/
745: PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)746: {
752: if (st->ops->backtransform) {
753: (*st->ops->backtransform)(st,n,eigr,eigi);
754: }
755: return(0);
756: }
758: /*@
759: STIsInjective - Ask if this spectral transformation is injective or not
760: (that is, if it corresponds to a one-to-one mapping). If not, then it
761: does not make sense to call STBackTransform().
763: Not collective
765: Input Parameter:
766: . st - the spectral transformation context
768: Output Parameter:
769: . is - the answer
771: Level: developer
773: .seealso: STBackTransform()
774: @*/
775: PetscErrorCode STIsInjective(ST st,PetscBool* is)776: {
778: PetscBool shell;
785: PetscObjectTypeCompare((PetscObject)st,STSHELL,&shell);
786: if (shell) {
787: STIsInjective_Shell(st,is);
788: } else *is = st->ops->backtransform? PETSC_TRUE: PETSC_FALSE;
789: return(0);
790: }
792: /*@
793: STMatSetUp - Build the preconditioner matrix used in STMatSolve().
795: Collective on st
797: Input Parameters:
798: + st - the spectral transformation context
799: . sigma - the shift
800: - coeffs - the coefficients (may be NULL)
802: Note:
803: This function is not intended to be called by end users, but by SLEPc
804: solvers that use ST. It builds matrix st->P as follows, then calls KSPSetUp().
805: .vb
806: If (coeffs): st->P = Sum_{i=0:nmat-1} coeffs[i]*sigma^i*A_i.
807: else st->P = Sum_{i=0:nmat-1} sigma^i*A_i
808: .ve
810: Level: developer
812: .seealso: STMatSolve()
813: @*/
814: PetscErrorCode STMatSetUp(ST st,PetscScalar sigma,PetscScalar *coeffs)815: {
821: STCheckMatrices(st,1);
823: PetscLogEventBegin(ST_MatSetUp,st,0,0,0);
824: STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,&st->P);
825: if (!st->ksp) { STGetKSP(st,&st->ksp); }
826: STCheckFactorPackage(st);
827: KSPSetOperators(st->ksp,st->P,st->P);
828: KSPSetUp(st->ksp);
829: PetscLogEventEnd(ST_MatSetUp,st,0,0,0);
830: return(0);
831: }
833: /*@
834: STSetWorkVecs - Sets a number of work vectors into the ST object.
836: Collective on st
838: Input Parameters:
839: + st - the spectral transformation context
840: - nw - number of work vectors to allocate
842: Developers Note:
843: This is SLEPC_EXTERN because it may be required by shell STs.
845: Level: developer
846: @*/
847: PetscErrorCode STSetWorkVecs(ST st,PetscInt nw)848: {
850: PetscInt i;
855: if (nw <= 0) SETERRQ1(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %D",nw);
856: if (st->nwork < nw) {
857: VecDestroyVecs(st->nwork,&st->work);
858: st->nwork = nw;
859: PetscMalloc1(nw,&st->work);
860: for (i=0;i<nw;i++) { STMatCreateVecs(st,&st->work[i],NULL); }
861: PetscLogObjectParents(st,nw,st->work);
862: }
863: return(0);
864: }