Actual source code: stsolve.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    ST interface routines, callable by users
 12: */

 14: #include <slepc/private/stimpl.h>

 16: PetscErrorCode STApply_Generic(ST st,Vec x,Vec y)
 17: {
 18:   if (st->M && st->P) {
 19:     MatMult(st->M,x,st->work[0]);
 20:     STMatSolve(st,st->work[0],y);
 21:   } else if (st->M) MatMult(st->M,x,y);
 22:   else STMatSolve(st,x,y);
 23:   PetscFunctionReturn(0);
 24: }

 26: /*@
 27:    STApply - Applies the spectral transformation operator to a vector, for
 28:    instance (A - sB)^-1 B in the case of the shift-and-invert transformation
 29:    and generalized eigenproblem.

 31:    Collective on st

 33:    Input Parameters:
 34: +  st - the spectral transformation context
 35: -  x  - input vector

 37:    Output Parameter:
 38: .  y - output vector

 40:    Level: developer

 42: .seealso: STApplyTranspose(), STApplyHermitianTranspose()
 43: @*/
 44: PetscErrorCode STApply(ST st,Vec x,Vec y)
 45: {
 46:   Mat            Op;

 52:   STCheckMatrices(st,1);
 54:   VecSetErrorIfLocked(y,3);
 56:   STGetOperator_Private(st,&Op);
 57:   MatMult(Op,x,y);
 58:   PetscFunctionReturn(0);
 59: }

 61: PetscErrorCode STApplyMat_Generic(ST st,Mat B,Mat C)
 62: {
 63:   Mat            work;

 65:   if (st->M && st->P) {
 66:     MatMatMult(st->M,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&work);
 67:     STMatMatSolve(st,work,C);
 68:     MatDestroy(&work);
 69:   } else if (st->M) MatMatMult(st->M,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
 70:   else STMatMatSolve(st,B,C);
 71:   PetscFunctionReturn(0);
 72: }

 74: /*@
 75:    STApplyMat - Applies the spectral transformation operator to a matrix, for
 76:    instance (A - sB)^-1 B in the case of the shift-and-invert transformation
 77:    and generalized eigenproblem.

 79:    Collective on st

 81:    Input Parameters:
 82: +  st - the spectral transformation context
 83: -  X  - input matrix

 85:    Output Parameter:
 86: .  Y - output matrix

 88:    Level: developer

 90: .seealso: STApply()
 91: @*/
 92: PetscErrorCode STApplyMat(ST st,Mat X,Mat Y)
 93: {
 98:   STCheckMatrices(st,1);
101:   (*st->ops->applymat)(st,X,Y);
102:   PetscFunctionReturn(0);
103: }

105: PetscErrorCode STApplyTranspose_Generic(ST st,Vec x,Vec y)
106: {
107:   if (st->M && st->P) {
108:     STMatSolveTranspose(st,x,st->work[0]);
109:     MatMultTranspose(st->M,st->work[0],y);
110:   } else if (st->M) MatMultTranspose(st->M,x,y);
111:   else STMatSolveTranspose(st,x,y);
112:   PetscFunctionReturn(0);
113: }

115: /*@
116:    STApplyTranspose - Applies the transpose of the operator to a vector, for
117:    instance B^T(A - sB)^-T in the case of the shift-and-invert transformation
118:    and generalized eigenproblem.

120:    Collective on st

122:    Input Parameters:
123: +  st - the spectral transformation context
124: -  x  - input vector

126:    Output Parameter:
127: .  y - output vector

129:    Level: developer

131: .seealso: STApply(), STApplyHermitianTranspose()
132: @*/
133: PetscErrorCode STApplyTranspose(ST st,Vec x,Vec y)
134: {
135:   Mat            Op;

141:   STCheckMatrices(st,1);
143:   VecSetErrorIfLocked(y,3);
145:   STGetOperator_Private(st,&Op);
146:   MatMultTranspose(Op,x,y);
147:   PetscFunctionReturn(0);
148: }

150: /*@
151:    STApplyHermitianTranspose - Applies the hermitian-transpose of the operator
152:    to a vector, for instance B^H(A - sB)^-H in the case of the shift-and-invert
153:    transformation and generalized eigenproblem.

155:    Collective on st

157:    Input Parameters:
158: +  st - the spectral transformation context
159: -  x  - input vector

161:    Output Parameter:
162: .  y - output vector

164:    Note:
165:    Currently implemented via STApplyTranspose() with appropriate conjugation.

167:    Level: developer

169: .seealso: STApply(), STApplyTranspose()
170: @*/
171: PetscErrorCode STApplyHermitianTranspose(ST st,Vec x,Vec y)
172: {
173:   Mat            Op;

179:   STCheckMatrices(st,1);
181:   VecSetErrorIfLocked(y,3);
183:   STGetOperator_Private(st,&Op);
184:   MatMultHermitianTranspose(Op,x,y);
185:   PetscFunctionReturn(0);
186: }

188: /*@
189:    STGetBilinearForm - Returns the matrix used in the bilinear form with a
190:    generalized problem with semi-definite B.

192:    Not collective, though a parallel Mat may be returned

194:    Input Parameters:
195: .  st - the spectral transformation context

197:    Output Parameter:
198: .  B - output matrix

200:    Notes:
201:    The output matrix B must be destroyed after use. It will be NULL in
202:    case of standard eigenproblems.

204:    Level: developer

206: .seealso: BVSetMatrix()
207: @*/
208: PetscErrorCode STGetBilinearForm(ST st,Mat *B)
209: {
213:   STCheckMatrices(st,1);
214:   (*st->ops->getbilinearform)(st,B);
215:   PetscFunctionReturn(0);
216: }

218: PetscErrorCode STGetBilinearForm_Default(ST st,Mat *B)
219: {
220:   if (st->nmat==1) *B = NULL;
221:   else {
222:     *B = st->A[1];
223:     PetscObjectReference((PetscObject)*B);
224:   }
225:   PetscFunctionReturn(0);
226: }

228: static PetscErrorCode MatMult_STOperator(Mat Op,Vec x,Vec y)
229: {
230:   ST             st;

232:   MatShellGetContext(Op,&st);
233:   STSetUp(st);
234:   PetscLogEventBegin(ST_Apply,st,x,y,0);
235:   if (st->D) { /* with balancing */
236:     VecPointwiseDivide(st->wb,x,st->D);
237:     (*st->ops->apply)(st,st->wb,y);
238:     VecPointwiseMult(y,y,st->D);
239:   } else (*st->ops->apply)(st,x,y);
240:   PetscLogEventEnd(ST_Apply,st,x,y,0);
241:   PetscFunctionReturn(0);
242: }

244: static PetscErrorCode MatMultTranspose_STOperator(Mat Op,Vec x,Vec y)
245: {
246:   ST             st;

248:   MatShellGetContext(Op,&st);
249:   STSetUp(st);
250:   PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);
251:   if (st->D) { /* with balancing */
252:     VecPointwiseMult(st->wb,x,st->D);
253:     (*st->ops->applytrans)(st,st->wb,y);
254:     VecPointwiseDivide(y,y,st->D);
255:   } else (*st->ops->applytrans)(st,x,y);
256:   PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);
257:   PetscFunctionReturn(0);
258: }

260: #if defined(PETSC_USE_COMPLEX)
261: static PetscErrorCode MatMultHermitianTranspose_STOperator(Mat Op,Vec x,Vec y)
262: {
263:   ST             st;

265:   MatShellGetContext(Op,&st);
266:   STSetUp(st);
267:   PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);
268:   if (!st->wht) {
269:     MatCreateVecs(st->A[0],&st->wht,NULL);
270:     PetscLogObjectParent((PetscObject)st,(PetscObject)st->wht);
271:   }
272:   VecCopy(x,st->wht);
273:   VecConjugate(st->wht);
274:   if (st->D) { /* with balancing */
275:     VecPointwiseMult(st->wb,st->wht,st->D);
276:     (*st->ops->applytrans)(st,st->wb,y);
277:     VecPointwiseDivide(y,y,st->D);
278:   } else (*st->ops->applytrans)(st,st->wht,y);
279:   VecConjugate(y);
280:   PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);
281:   PetscFunctionReturn(0);
282: }
283: #endif

285: static PetscErrorCode MatMatMult_STOperator(Mat Op,Mat B,Mat C,void *ctx)
286: {
287:   ST             st;

289:   MatShellGetContext(Op,&st);
290:   STSetUp(st);
291:   PetscLogEventBegin(ST_Apply,st,B,C,0);
292:   STApplyMat_Generic(st,B,C);
293:   PetscLogEventEnd(ST_Apply,st,B,C,0);
294:   PetscFunctionReturn(0);
295: }

297: PetscErrorCode STGetOperator_Private(ST st,Mat *Op)
298: {
299:   PetscInt       m,n,M,N;
300:   Vec            v;
301:   VecType        vtype;

303:   if (!st->Op) {
304:     if (Op) *Op = NULL;
305:     /* create the shell matrix */
306:     MatGetLocalSize(st->A[0],&m,&n);
307:     MatGetSize(st->A[0],&M,&N);
308:     MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,st,&st->Op);
309:     MatShellSetOperation(st->Op,MATOP_MULT,(void(*)(void))MatMult_STOperator);
310:     MatShellSetOperation(st->Op,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator);
311: #if defined(PETSC_USE_COMPLEX)
312:     MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_STOperator);
313: #else
314:     MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator);
315: #endif
316:     if (!st->D && st->ops->apply==STApply_Generic) {
317:       MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSE,MATDENSE);
318:       MatShellSetMatProductOperation(st->Op,MATPRODUCT_AB,NULL,MatMatMult_STOperator,NULL,MATDENSECUDA,MATDENSECUDA);
319:     }
320:     /* make sure the shell matrix generates a vector of the same type as the problem matrices */
321:     MatCreateVecs(st->A[0],&v,NULL);
322:     VecGetType(v,&vtype);
323:     MatShellSetVecType(st->Op,vtype);
324:     VecDestroy(&v);
325:     /* build the operator matrices */
326:     STComputeOperator(st);
327:   }
328:   if (Op) *Op = st->Op;
329:   PetscFunctionReturn(0);
330: }

332: /*@
333:    STGetOperator - Returns a shell matrix that represents the operator of the
334:    spectral transformation.

336:    Collective on st

338:    Input Parameter:
339: .  st - the spectral transformation context

341:    Output Parameter:
342: .  Op - operator matrix

344:    Notes:
345:    The operator is defined in linear eigenproblems only, not in polynomial ones,
346:    so the call will fail if more than 2 matrices were passed in STSetMatrices().

348:    The returned shell matrix is essentially a wrapper to the STApply() and
349:    STApplyTranspose() operations. The operator can often be expressed as

351: $     Op = D*inv(K)*M*inv(D)

353:    where D is the balancing matrix, and M and K are two matrices corresponding
354:    to the numerator and denominator for spectral transformations that represent
355:    a rational matrix function. In the case of STSHELL, the inner part inv(K)*M
356:    is replaced by the user-provided operation from STShellSetApply().

358:    The preconditioner matrix K typically depends on the value of the shift, and
359:    its inverse is handled via an internal KSP object. Normal usage does not
360:    require explicitly calling STGetOperator(), but it can be used to force the
361:    creation of K and M, and then K is passed to the KSP. This is useful for
362:    setting options associated with the PCFactor (to set MUMPS options, for instance).

364:    The returned matrix must NOT be destroyed by the user. Instead, when no
365:    longer needed it must be returned with STRestoreOperator(). In particular,
366:    this is required before modifying the ST matrices or the shift.

368:    A NULL pointer can be passed in Op in case the matrix is not required but we
369:    want to force its creation. In this case, STRestoreOperator() should not be
370:    called.

372:    Level: advanced

374: .seealso: STApply(), STApplyTranspose(), STSetBalanceMatrix(), STShellSetApply(),
375:           STGetKSP(), STSetShift(), STRestoreOperator(), STSetMatrices()
376: @*/
377: PetscErrorCode STGetOperator(ST st,Mat *Op)
378: {
381:   STCheckMatrices(st,1);
382:   STCheckNotSeized(st,1);
384:   STGetOperator_Private(st,Op);
385:   if (Op) st->opseized = PETSC_TRUE;
386:   PetscFunctionReturn(0);
387: }

389: /*@
390:    STRestoreOperator - Restore the previously seized operator matrix.

392:    Collective on st

394:    Input Parameters:
395: +  st - the spectral transformation context
396: -  Op - operator matrix

398:    Notes:
399:    The arguments must match the corresponding call to STGetOperator().

401:    Level: advanced

403: .seealso: STGetOperator()
404: @*/
405: PetscErrorCode STRestoreOperator(ST st,Mat *Op)
406: {
411:   *Op = NULL;
412:   st->opseized = PETSC_FALSE;
413:   PetscFunctionReturn(0);
414: }

416: /*
417:    STComputeOperator - Computes the matrices that constitute the operator

419:       Op = D*inv(K)*M*inv(D).

421:    K and M are computed here (D is user-provided) from the system matrices
422:    and the shift sigma (whenever these are changed, this function recomputes
423:    K and M). This is used only in linear eigenproblems (nmat<3).

425:    K is the "preconditioner matrix": it is the denominator in rational operators,
426:    e.g. (A-sigma*B) in shift-and-invert. In non-rational transformations such
427:    as STFILTER, K=NULL which means identity. After computing K, it is passed to
428:    the internal KSP object via KSPSetOperators.

430:    M is the numerator in rational operators. If unused it is set to NULL (e.g.
431:    in STPRECOND).

433:    STSHELL does not compute anything here, but sets the flag as if it was ready.
434: */
435: PetscErrorCode STComputeOperator(ST st)
436: {
437:   PC             pc;

441:   if (!st->opready && st->ops->computeoperator) {
442:     PetscInfo(st,"Building the operator matrices\n");
443:     STCheckMatrices(st,1);
444:     if (!st->T) {
445:       PetscCalloc1(PetscMax(2,st->nmat),&st->T);
446:       PetscLogObjectMemory((PetscObject)st,PetscMax(2,st->nmat)*sizeof(Mat));
447:     }
448:     PetscLogEventBegin(ST_ComputeOperator,st,0,0,0);
449:     (*st->ops->computeoperator)(st);
450:     PetscLogEventEnd(ST_ComputeOperator,st,0,0,0);
451:     if (st->usesksp) {
452:       if (!st->ksp) STGetKSP(st,&st->ksp);
453:       if (st->P) {
454:         STSetDefaultKSP(st);
455:         ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P);
456:       } else {
457:         /* STPRECOND defaults to PCNONE if st->P is empty */
458:         KSPGetPC(st->ksp,&pc);
459:         PCSetType(pc,PCNONE);
460:       }
461:     }
462:   }
463:   st->opready = PETSC_TRUE;
464:   PetscFunctionReturn(0);
465: }

467: /*@
468:    STSetUp - Prepares for the use of a spectral transformation.

470:    Collective on st

472:    Input Parameter:
473: .  st - the spectral transformation context

475:    Level: advanced

477: .seealso: STCreate(), STApply(), STDestroy()
478: @*/
479: PetscErrorCode STSetUp(ST st)
480: {
481:   PetscInt       i,n,k;

485:   STCheckMatrices(st,1);
486:   switch (st->state) {
487:     case ST_STATE_INITIAL:
488:       PetscInfo(st,"Setting up new ST\n");
489:       if (!((PetscObject)st)->type_name) STSetType(st,STSHIFT);
490:       break;
491:     case ST_STATE_SETUP:
492:       PetscFunctionReturn(0);
493:     case ST_STATE_UPDATED:
494:       PetscInfo(st,"Setting up updated ST\n");
495:       break;
496:   }
497:   PetscLogEventBegin(ST_SetUp,st,0,0,0);
498:   if (st->state!=ST_STATE_UPDATED) {
499:     if (!(st->nmat<3 && st->opready)) {
500:       if (st->T) {
501:         for (i=0;i<PetscMax(2,st->nmat);i++) MatDestroy(&st->T[i]);
502:       }
503:       MatDestroy(&st->P);
504:     }
505:   }
506:   if (st->D) {
507:     MatGetLocalSize(st->A[0],NULL,&n);
508:     VecGetLocalSize(st->D,&k);
510:     if (!st->wb) {
511:       VecDuplicate(st->D,&st->wb);
512:       PetscLogObjectParent((PetscObject)st,(PetscObject)st->wb);
513:     }
514:   }
515:   if (st->nmat<3 && st->transform) STComputeOperator(st);
516:   else {
517:     if (!st->T) {
518:       PetscCalloc1(PetscMax(2,st->nmat),&st->T);
519:       PetscLogObjectMemory((PetscObject)st,PetscMax(2,st->nmat)*sizeof(Mat));
520:     }
521:   }
522:   if (st->ops->setup) (*st->ops->setup)(st);
523:   st->state = ST_STATE_SETUP;
524:   PetscLogEventEnd(ST_SetUp,st,0,0,0);
525:   PetscFunctionReturn(0);
526: }

528: /*
529:    Computes coefficients for the transformed polynomial,
530:    and stores the result in argument S.

532:    alpha - value of the parameter of the transformed polynomial
533:    beta - value of the previous shift (only used in inplace mode)
534:    k - index of first matrix included in the computation
535:    coeffs - coefficients of the expansion
536:    initial - true if this is the first time
537:    precond - whether the preconditioner matrix must be computed
538: */
539: PetscErrorCode STMatMAXPY_Private(ST st,PetscScalar alpha,PetscScalar beta,PetscInt k,PetscScalar *coeffs,PetscBool initial,PetscBool precond,Mat *S)
540: {
541:   PetscInt       *matIdx=NULL,nmat,i,ini=-1;
542:   PetscScalar    t=1.0,ta,gamma;
543:   PetscBool      nz=PETSC_FALSE;
544:   Mat            *A=precond?st->Psplit:st->A;
545:   MatStructure   str=precond?st->strp:st->str;

547:   nmat = st->nmat-k;
548:   switch (st->matmode) {
549:   case ST_MATMODE_INPLACE:
552:     if (initial) {
553:       PetscObjectReference((PetscObject)A[0]);
554:       *S = A[0];
555:       gamma = alpha;
556:     } else gamma = alpha-beta;
557:     if (gamma != 0.0) {
558:       if (st->nmat>1) MatAXPY(*S,gamma,A[1],str);
559:       else MatShift(*S,gamma);
560:     }
561:     break;
562:   case ST_MATMODE_SHELL:
564:     if (initial) {
565:       if (st->nmat>2) {
566:         PetscMalloc1(nmat,&matIdx);
567:         for (i=0;i<nmat;i++) matIdx[i] = k+i;
568:       }
569:       STMatShellCreate(st,alpha,nmat,matIdx,coeffs,S);
570:       PetscLogObjectParent((PetscObject)st,(PetscObject)*S);
571:       if (st->nmat>2) PetscFree(matIdx);
572:     } else STMatShellShift(*S,alpha);
573:     break;
574:   case ST_MATMODE_COPY:
575:     if (coeffs) {
576:       for (i=0;i<nmat && ini==-1;i++) {
577:         if (coeffs[i]!=0.0) ini = i;
578:         else t *= alpha;
579:       }
580:       if (coeffs[ini] != 1.0) nz = PETSC_TRUE;
581:       for (i=ini+1;i<nmat&&!nz;i++) if (coeffs[i]!=0.0) nz = PETSC_TRUE;
582:     } else { nz = PETSC_TRUE; ini = 0; }
583:     if ((alpha == 0.0 || !nz) && t==1.0) {
584:       PetscObjectReference((PetscObject)A[k+ini]);
585:       MatDestroy(S);
586:       *S = A[k+ini];
587:     } else {
588:       if (*S && *S!=A[k+ini]) {
589:         MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
590:         MatCopy(A[k+ini],*S,DIFFERENT_NONZERO_PATTERN);
591:       } else {
592:         MatDestroy(S);
593:         MatDuplicate(A[k+ini],MAT_COPY_VALUES,S);
594:         MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
595:         PetscLogObjectParent((PetscObject)st,(PetscObject)*S);
596:       }
597:       if (coeffs && coeffs[ini]!=1.0) MatScale(*S,coeffs[ini]);
598:       for (i=ini+k+1;i<PetscMax(2,st->nmat);i++) {
599:         t *= alpha;
600:         ta = t;
601:         if (coeffs) ta *= coeffs[i-k];
602:         if (ta!=0.0) {
603:           if (st->nmat>1) MatAXPY(*S,ta,A[i],str);
604:           else MatShift(*S,ta);
605:         }
606:       }
607:     }
608:   }
609:   MatSetOption(*S,MAT_SYMMETRIC,st->asymm);
610:   MatSetOption(*S,MAT_HERMITIAN,(PetscImaginaryPart(st->sigma)==0.0)?st->aherm:PETSC_FALSE);
611:   PetscFunctionReturn(0);
612: }

614: /*
615:    Computes the values of the coefficients required by STMatMAXPY_Private
616:    for the case of monomial basis.
617: */
618: PetscErrorCode STCoeffs_Monomial(ST st, PetscScalar *coeffs)
619: {
620:   PetscInt  k,i,ini,inip;

622:   /* Compute binomial coefficients */
623:   ini = (st->nmat*(st->nmat-1))/2;
624:   for (i=0;i<st->nmat;i++) coeffs[ini+i]=1.0;
625:   for (k=st->nmat-1;k>=1;k--) {
626:     inip = ini+1;
627:     ini = (k*(k-1))/2;
628:     coeffs[ini] = 1.0;
629:     for (i=1;i<k;i++) coeffs[ini+i] = coeffs[ini+i-1]+coeffs[inip+i-1];
630:   }
631:   PetscFunctionReturn(0);
632: }

634: /*@
635:    STPostSolve - Optional post-solve phase, intended for any actions that must
636:    be performed on the ST object after the eigensolver has finished.

638:    Collective on st

640:    Input Parameters:
641: .  st  - the spectral transformation context

643:    Level: developer

645: .seealso: EPSSolve()
646: @*/
647: PetscErrorCode STPostSolve(ST st)
648: {
651:   if (st->ops->postsolve) (*st->ops->postsolve)(st);
652:   PetscFunctionReturn(0);
653: }

655: /*@
656:    STBackTransform - Back-transformation phase, intended for
657:    spectral transformations which require to transform the computed
658:    eigenvalues back to the original eigenvalue problem.

660:    Not Collective

662:    Input Parameters:
663: +  st   - the spectral transformation context
664: .  n    - number of eigenvalues
665: .  eigr - real part of a computed eigenvalues
666: -  eigi - imaginary part of a computed eigenvalues

668:    Level: developer

670: .seealso: STIsInjective()
671: @*/
672: PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
673: {
676:   if (st->ops->backtransform) (*st->ops->backtransform)(st,n,eigr,eigi);
677:   PetscFunctionReturn(0);
678: }

680: /*@
681:    STIsInjective - Ask if this spectral transformation is injective or not
682:    (that is, if it corresponds to a one-to-one mapping). If not, then it
683:    does not make sense to call STBackTransform().

685:    Not collective

687:    Input Parameter:
688: .  st   - the spectral transformation context

690:    Output Parameter:
691: .  is - the answer

693:    Level: developer

695: .seealso: STBackTransform()
696: @*/
697: PetscErrorCode STIsInjective(ST st,PetscBool* is)
698: {
699:   PetscBool      shell;


705:   PetscObjectTypeCompare((PetscObject)st,STSHELL,&shell);
706:   if (shell) STIsInjective_Shell(st,is);
707:   else *is = st->ops->backtransform? PETSC_TRUE: PETSC_FALSE;
708:   PetscFunctionReturn(0);
709: }

711: /*@
712:    STMatSetUp - Build the preconditioner matrix used in STMatSolve().

714:    Collective on st

716:    Input Parameters:
717: +  st     - the spectral transformation context
718: .  sigma  - the shift
719: -  coeffs - the coefficients (may be NULL)

721:    Note:
722:    This function is not intended to be called by end users, but by SLEPc
723:    solvers that use ST. It builds matrix st->P as follows, then calls KSPSetUp().
724: .vb
725:     If (coeffs)  st->P = Sum_{i=0..nmat-1} coeffs[i]*sigma^i*A_i
726:     else         st->P = Sum_{i=0..nmat-1} sigma^i*A_i
727: .ve

729:    Level: developer

731: .seealso: STMatSolve()
732: @*/
733: PetscErrorCode STMatSetUp(ST st,PetscScalar sigma,PetscScalar *coeffs)
734: {
737:   STCheckMatrices(st,1);

739:   PetscLogEventBegin(ST_MatSetUp,st,0,0,0);
740:   STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_FALSE,&st->P);
741:   if (st->Psplit) STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_TRUE,&st->Pmat);
742:   ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P);
743:   KSPSetUp(st->ksp);
744:   PetscLogEventEnd(ST_MatSetUp,st,0,0,0);
745:   PetscFunctionReturn(0);
746: }

748: /*@
749:    STSetWorkVecs - Sets a number of work vectors into the ST object.

751:    Collective on st

753:    Input Parameters:
754: +  st - the spectral transformation context
755: -  nw - number of work vectors to allocate

757:    Developer Notes:
758:    This is SLEPC_EXTERN because it may be required by shell STs.

760:    Level: developer

762: .seealso: STMatCreateVecs()
763: @*/
764: PetscErrorCode STSetWorkVecs(ST st,PetscInt nw)
765: {
766:   PetscInt       i;

771:   if (st->nwork < nw) {
772:     VecDestroyVecs(st->nwork,&st->work);
773:     st->nwork = nw;
774:     PetscMalloc1(nw,&st->work);
775:     for (i=0;i<nw;i++) STMatCreateVecs(st,&st->work[i],NULL);
776:     PetscLogObjectParents(st,nw,st->work);
777:   }
778:   PetscFunctionReturn(0);
779: }