Actual source code: stsolve.c

slepc-3.13.3 2020-06-14
Report Typos and Errors
  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 STApplyTranspose_Generic(ST st,Vec x,Vec y)
 70: {

 74:   if (st->M && st->P) {
 75:     STMatSolveTranspose(st,x,st->work[0]);
 76:     MatMultTranspose(st->M,st->work[0],y);
 77:   } else if (st->M) {
 78:     MatMultTranspose(st->M,x,y);
 79:   } else {
 80:     STMatSolveTranspose(st,x,y);
 81:   }
 82:   return(0);
 83: }

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

 90:    Collective on st

 92:    Input Parameters:
 93: +  st - the spectral transformation context
 94: -  x  - input vector

 96:    Output Parameter:
 97: .  y - output vector

 99:    Level: developer

101: .seealso: STApply(), STApplyHermitianTranspose()
102: @*/
103: PetscErrorCode STApplyTranspose(ST st,Vec x,Vec y)
104: {
106:   Mat            Op;

113:   STCheckMatrices(st,1);
114:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
115:   VecSetErrorIfLocked(y,3);
116:   if (!st->ops->applytrans) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST does not have applytrans");
117:   STGetOperator_Private(st,&Op);
118:   MatMultTranspose(Op,x,y);
119:   return(0);
120: }

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

127:    Collective on st

129:    Input Parameters:
130: +  st - the spectral transformation context
131: -  x  - input vector

133:    Output Parameter:
134: .  y - output vector

136:    Note:
137:    Currently implemented via STApplyTranspose() with appropriate conjugation.

139:    Level: developer

141: .seealso: STApply(), STApplyTranspose()
142: @*/
143: PetscErrorCode STApplyHermitianTranspose(ST st,Vec x,Vec y)
144: {
146:   Mat            Op;

153:   STCheckMatrices(st,1);
154:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
155:   VecSetErrorIfLocked(y,3);
156:   if (!st->ops->applytrans) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST does not have applytrans");
157:   STGetOperator_Private(st,&Op);
158:   MatMultHermitianTranspose(Op,x,y);
159:   return(0);
160: }

162: /*@
163:    STGetBilinearForm - Returns the matrix used in the bilinear form with a
164:    generalized problem with semi-definite B.

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

168:    Input Parameters:
169: .  st - the spectral transformation context

171:    Output Parameter:
172: .  B - output matrix

174:    Notes:
175:    The output matrix B must be destroyed after use. It will be NULL in
176:    case of standard eigenproblems.

178:    Level: developer
179: @*/
180: PetscErrorCode STGetBilinearForm(ST st,Mat *B)
181: {

188:   STCheckMatrices(st,1);
189:   (*st->ops->getbilinearform)(st,B);
190:   return(0);
191: }

193: PetscErrorCode STGetBilinearForm_Default(ST st,Mat *B)
194: {

198:   if (st->nmat==1) *B = NULL;
199:   else {
200:     *B = st->A[1];
201:     PetscObjectReference((PetscObject)*B);
202:   }
203:   return(0);
204: }

206: static PetscErrorCode MatMult_STOperator(Mat Op,Vec x,Vec y)
207: {
209:   ST             st;

212:   MatShellGetContext(Op,(void**)&st);
213:   STSetUp(st);
214:   PetscLogEventBegin(ST_Apply,st,x,y,0);
215:   if (st->D) { /* with balancing */
216:     VecPointwiseDivide(st->wb,x,st->D);
217:     (*st->ops->apply)(st,st->wb,y);
218:     VecPointwiseMult(y,y,st->D);
219:   } else {
220:     (*st->ops->apply)(st,x,y);
221:   }
222:   PetscLogEventEnd(ST_Apply,st,x,y,0);
223:   return(0);
224: }

226: static PetscErrorCode MatMultTranspose_STOperator(Mat Op,Vec x,Vec y)
227: {
229:   ST             st;

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

246: #if defined(PETSC_USE_COMPLEX)
247: static PetscErrorCode MatMultHermitianTranspose_STOperator(Mat Op,Vec x,Vec y)
248: {
250:   ST             st;

253:   MatShellGetContext(Op,(void**)&st);
254:   STSetUp(st);
255:   PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);
256:   if (!st->wht) {
257:     MatCreateVecs(st->A[0],&st->wht,NULL);
258:     PetscLogObjectParent((PetscObject)st,(PetscObject)st->wht);
259:   }
260:   VecCopy(x,st->wht);
261:   VecConjugate(st->wht);
262:   if (st->D) { /* with balancing */
263:     VecPointwiseMult(st->wb,st->wht,st->D);
264:     (*st->ops->applytrans)(st,st->wb,y);
265:     VecPointwiseDivide(y,y,st->D);
266:   } else {
267:     (*st->ops->applytrans)(st,st->wht,y);
268:   }
269:   VecConjugate(y);
270:   PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);
271:   return(0);
272: }
273: #endif

275: PetscErrorCode STGetOperator_Private(ST st,Mat *Op)
276: {
278:   PetscInt       m,n,M,N;

281:   if (!st->Op) {
282:     if (Op) *Op = NULL;
283:     MatGetLocalSize(st->A[0],&m,&n);
284:     MatGetSize(st->A[0],&M,&N);
285:     MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,st,&st->Op);
286:     MatShellSetOperation(st->Op,MATOP_MULT,(void(*)(void))MatMult_STOperator);
287:     MatShellSetOperation(st->Op,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator);
288: #if defined(PETSC_USE_COMPLEX)
289:     MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_STOperator);
290: #else
291:     MatShellSetOperation(st->Op,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator);
292: #endif
293:     STComputeOperator(st);
294:   }
295:   if (Op) *Op = st->Op;
296:   return(0);
297: }

299: /*@
300:    STGetOperator - Returns a shell matrix that represents the operator of the
301:    spectral transformation.

303:    Collective on st

305:    Input Parameter:
306: .  st - the spectral transformation context

308:    Output Parameter:
309: .  Op - operator matrix

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

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

318: .vb
319:       Op = D*inv(K)*M*inv(D)
320: .ve

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

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

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

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

341:    Level: advanced

343: .seealso: STApply(), STApplyTranspose(), STSetBalanceMatrix(), STShellSetApply(),
344:           STGetKSP(), STSetShift(), STRestoreOperator(), STSetMatrices()
345: @*/
346: PetscErrorCode STGetOperator(ST st,Mat *Op)
347: {

353:   STCheckMatrices(st,1);
354:   STCheckNotSeized(st,1);
355:   if (st->nmat>2) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"The operator is not defined in polynomial eigenproblems");
356:   STGetOperator_Private(st,Op);
357:   if (Op) st->opseized = PETSC_TRUE;
358:   return(0);
359: }

361: /*@
362:    STRestoreOperator - Restore the previously seized operator matrix.

364:    Collective on st

366:    Input Parameters:
367: +  st - the spectral transformation context
368: -  Op - operator matrix

370:    Notes:
371:    The arguments must match the corresponding call to STGetOperator().

373:    Level: advanced

375: .seealso: STGetOperator()
376: @*/
377: PetscErrorCode STRestoreOperator(ST st,Mat *Op)
378: {
383:   if (!st->opseized) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must be called after STGetOperator()");
384:   *Op = NULL;
385:   st->opseized = PETSC_FALSE;
386:   return(0);
387: }

389: /*
390:    STComputeOperator - Computes the matrices that constitute the operator

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

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

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

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

406:    STSHELL does not compute anything here, but sets the flag as if it was ready.
407: */
408: PetscErrorCode STComputeOperator(ST st)
409: {
411:   PC             pc;

416:   if (!st->opready && st->ops->computeoperator) {
417:     STCheckMatrices(st,1);
418:     if (!st->T) {
419:       PetscCalloc1(PetscMax(2,st->nmat),&st->T);
420:       PetscLogObjectMemory((PetscObject)st,PetscMax(2,st->nmat)*sizeof(Mat));
421:     }
422:     PetscLogEventBegin(ST_ComputeOperator,st,0,0,0);
423:     (*st->ops->computeoperator)(st);
424:     PetscLogEventEnd(ST_ComputeOperator,st,0,0,0);
425:     if (st->usesksp) {
426:       if (!st->ksp) { STGetKSP(st,&st->ksp); }
427:       if (st->P) {
428:         STSetDefaultKSP(st);
429:         STCheckFactorPackage(st);
430:         KSPSetOperators(st->ksp,st->P,st->P);
431:       } else {
432:         /* STPRECOND defaults to PCNONE if st->P is empty */
433:         KSPGetPC(st->ksp,&pc);
434:         PCSetType(pc,PCNONE);
435:       }
436:     }
437:   }
438:   st->opready = PETSC_TRUE;
439:   return(0);
440: }

442: /*@
443:    STSetUp - Prepares for the use of a spectral transformation.

445:    Collective on st

447:    Input Parameter:
448: .  st - the spectral transformation context

450:    Level: advanced

452: .seealso: STCreate(), STApply(), STDestroy()
453: @*/
454: PetscErrorCode STSetUp(ST st)
455: {
456:   PetscInt       i,n,k;

462:   STCheckMatrices(st,1);
463:   switch (st->state) {
464:     case ST_STATE_INITIAL:
465:       PetscInfo(st,"Setting up new ST\n");
466:       if (!((PetscObject)st)->type_name) {
467:         STSetType(st,STSHIFT);
468:       }
469:       break;
470:     case ST_STATE_SETUP:
471:       return(0);
472:     case ST_STATE_UPDATED:
473:       PetscInfo(st,"Setting up updated ST\n");
474:       break;
475:   }
476:   PetscLogEventBegin(ST_SetUp,st,0,0,0);
477:   if (st->state!=ST_STATE_UPDATED) {
478:     if (!(st->nmat<3 && st->opready)) {
479:       if (st->T) {
480:         for (i=0;i<PetscMax(2,st->nmat);i++) {
481:           MatDestroy(&st->T[i]);
482:         }
483:       }
484:       MatDestroy(&st->P);
485:     }
486:   }
487:   if (st->D) {
488:     MatGetLocalSize(st->A[0],NULL,&n);
489:     VecGetLocalSize(st->D,&k);
490:     if (n != k) SETERRQ2(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Balance matrix has wrong dimension %D (should be %D)",k,n);
491:     if (!st->wb) {
492:       VecDuplicate(st->D,&st->wb);
493:       PetscLogObjectParent((PetscObject)st,(PetscObject)st->wb);
494:     }
495:   }
496:   if (st->nmat<3 && st->transform) {
497:     STComputeOperator(st);
498:   } else {
499:     if (!st->T) {
500:       PetscCalloc1(PetscMax(2,st->nmat),&st->T);
501:       PetscLogObjectMemory((PetscObject)st,PetscMax(2,st->nmat)*sizeof(Mat));
502:     }
503:   }
504:   if (st->ops->setup) { (*st->ops->setup)(st); }
505:   st->state = ST_STATE_SETUP;
506:   PetscLogEventEnd(ST_SetUp,st,0,0,0);
507:   return(0);
508: }

510: /*
511:    Computes coefficients for the transformed polynomial,
512:    and stores the result in argument S.

514:    alpha - value of the parameter of the transformed polynomial
515:    beta - value of the previous shift (only used in inplace mode)
516:    k - index of first matrix included in the computation
517:    coeffs - coefficients of the expansion
518:    initial - true if this is the first time (only relevant for shell mode)
519: */
520: PetscErrorCode STMatMAXPY_Private(ST st,PetscScalar alpha,PetscScalar beta,PetscInt k,PetscScalar *coeffs,PetscBool initial,Mat *S)
521: {
523:   PetscInt       *matIdx=NULL,nmat,i,ini=-1;
524:   PetscScalar    t=1.0,ta,gamma;
525:   PetscBool      nz=PETSC_FALSE;

528:   nmat = st->nmat-k;
529:   switch (st->matmode) {
530:   case ST_MATMODE_INPLACE:
531:     if (st->nmat>2) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for polynomial eigenproblems");
532:     if (initial) {
533:       PetscObjectReference((PetscObject)st->A[0]);
534:       *S = st->A[0];
535:       gamma = alpha;
536:     } else gamma = alpha-beta;
537:     if (gamma != 0.0) {
538:       if (st->nmat>1) {
539:         MatAXPY(*S,gamma,st->A[1],st->str);
540:       } else {
541:         MatShift(*S,gamma);
542:       }
543:     }
544:     break;
545:   case ST_MATMODE_SHELL:
546:     if (initial) {
547:       if (st->nmat>2) {
548:         PetscMalloc1(nmat,&matIdx);
549:         for (i=0;i<nmat;i++) matIdx[i] = k+i;
550:       }
551:       STMatShellCreate(st,alpha,nmat,matIdx,coeffs,S);
552:       PetscLogObjectParent((PetscObject)st,(PetscObject)*S);
553:       if (st->nmat>2) { PetscFree(matIdx); }
554:     } else {
555:       STMatShellShift(*S,alpha);
556:     }
557:     break;
558:   case ST_MATMODE_COPY:
559:     if (coeffs) {
560:       for (i=0;i<nmat && ini==-1;i++) {
561:         if (coeffs[i]!=0.0) ini = i;
562:         else t *= alpha;
563:       }
564:       if (coeffs[ini] != 1.0) nz = PETSC_TRUE;
565:       for (i=ini+1;i<nmat&&!nz;i++) if (coeffs[i]!=0.0) nz = PETSC_TRUE;
566:     } else { nz = PETSC_TRUE; ini = 0; }
567:     if ((alpha == 0.0 || !nz) && t==1.0) {
568:       PetscObjectReference((PetscObject)st->A[k+ini]);
569:       MatDestroy(S);
570:       *S = st->A[k+ini];
571:     } else {
572:       if (*S && *S!=st->A[k+ini]) {
573:         MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
574:         MatCopy(st->A[k+ini],*S,DIFFERENT_NONZERO_PATTERN);
575:       } else {
576:         MatDestroy(S);
577:         MatDuplicate(st->A[k+ini],MAT_COPY_VALUES,S);
578:         MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
579:         PetscLogObjectParent((PetscObject)st,(PetscObject)*S);
580:       }
581:       if (coeffs && coeffs[ini]!=1.0) {
582:         MatScale(*S,coeffs[ini]);
583:       }
584:       for (i=ini+k+1;i<PetscMax(2,st->nmat);i++) {
585:         t *= alpha;
586:         ta = t;
587:         if (coeffs) ta *= coeffs[i-k];
588:         if (ta!=0.0) {
589:           if (st->nmat>1) {
590:             MatAXPY(*S,ta,st->A[i],st->str);
591:           } else {
592:             MatShift(*S,ta);
593:           }
594:         }
595:       }
596:     }
597:   }
598:   MatSetOption(*S,MAT_SYMMETRIC,st->asymm);
599:   MatSetOption(*S,MAT_HERMITIAN,(PetscImaginaryPart(st->sigma)==0.0)?st->aherm:PETSC_FALSE);
600:   return(0);
601: }

603: /*
604:    Computes the values of the coefficients required by STMatMAXPY_Private
605:    for the case of monomial basis.
606: */
607: PetscErrorCode STCoeffs_Monomial(ST st, PetscScalar *coeffs)
608: {
609:   PetscInt  k,i,ini,inip;

612:   /* Compute binomial coefficients */
613:   ini = (st->nmat*(st->nmat-1))/2;
614:   for (i=0;i<st->nmat;i++) coeffs[ini+i]=1.0;
615:   for (k=st->nmat-1;k>=1;k--) {
616:     inip = ini+1;
617:     ini = (k*(k-1))/2;
618:     coeffs[ini] = 1.0;
619:     for (i=1;i<k;i++) coeffs[ini+i] = coeffs[ini+i-1]+coeffs[inip+i-1];
620:   }
621:   return(0);
622: }

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

628:    Collective on st

630:    Input Parameters:
631: .  st  - the spectral transformation context

633:    Level: developer

635: .seealso: EPSSolve()
636: @*/
637: PetscErrorCode STPostSolve(ST st)
638: {

644:   if (st->ops->postsolve) {
645:     (*st->ops->postsolve)(st);
646:   }
647:   return(0);
648: }

650: /*@
651:    STBackTransform - Back-transformation phase, intended for
652:    spectral transformations which require to transform the computed
653:    eigenvalues back to the original eigenvalue problem.

655:    Not Collective

657:    Input Parameters:
658: +  st   - the spectral transformation context
659:    eigr - real part of a computed eigenvalue
660: -  eigi - imaginary part of a computed eigenvalue

662:    Level: developer

664: .seealso: STIsInjective()
665: @*/
666: PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
667: {

673:   if (st->ops->backtransform) {
674:     (*st->ops->backtransform)(st,n,eigr,eigi);
675:   }
676:   return(0);
677: }

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

684:    Not collective

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

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

692:    Level: developer

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


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

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

716:    Collective on st

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

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

731:    Level: developer

733: .seealso: STMatSolve()
734: @*/
735: PetscErrorCode STMatSetUp(ST st,PetscScalar sigma,PetscScalar *coeffs)
736: {

742:   STCheckMatrices(st,1);

744:   PetscLogEventBegin(ST_MatSetUp,st,0,0,0);
745:   STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,&st->P);
746:   if (!st->ksp) { STGetKSP(st,&st->ksp); }
747:   STCheckFactorPackage(st);
748:   KSPSetOperators(st->ksp,st->P,st->P);
749:   KSPSetUp(st->ksp);
750:   PetscLogEventEnd(ST_MatSetUp,st,0,0,0);
751:   return(0);
752: }

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

757:    Collective on st

759:    Input Parameters:
760: +  st - the spectral transformation context
761: -  nw - number of work vectors to allocate

763:    Developers Note:
764:    This is SLEPC_EXTERN because it may be required by shell STs.

766:    Level: developer
767: @*/
768: PetscErrorCode STSetWorkVecs(ST st,PetscInt nw)
769: {
771:   PetscInt       i;

776:   if (nw <= 0) SETERRQ1(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %D",nw);
777:   if (st->nwork < nw) {
778:     VecDestroyVecs(st->nwork,&st->work);
779:     st->nwork = nw;
780:     PetscMalloc1(nw,&st->work);
781:     for (i=0;i<nw;i++) { STMatCreateVecs(st,&st->work[i],NULL); }
782:     PetscLogObjectParents(st,nw,st->work);
783:   }
784:   return(0);
785: }