Actual source code: stsolve.c

slepc-3.18.2 2023-01-26
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:   return 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);
 55:   STGetOperator_Private(st,&Op);
 56:   MatMult(Op,x,y);
 57:   return 0;
 58: }

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

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

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

 78:    Collective on st

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

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

 87:    Level: developer

 89: .seealso: STApply()
 90: @*/
 91: PetscErrorCode STApplyMat(ST st,Mat X,Mat Y)
 92: {
 97:   STCheckMatrices(st,1);
 99:   PetscUseTypeMethod(st,applymat,X,Y);
100:   return 0;
101: }

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

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

118:    Collective on st

120:    Input Parameters:
121: +  st - the spectral transformation context
122: -  x  - input vector

124:    Output Parameter:
125: .  y - output vector

127:    Level: developer

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

139:   STCheckMatrices(st,1);
141:   VecSetErrorIfLocked(y,3);
142:   STGetOperator_Private(st,&Op);
143:   MatMultTranspose(Op,x,y);
144:   return 0;
145: }

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

152:    Collective on st

154:    Input Parameters:
155: +  st - the spectral transformation context
156: -  x  - input vector

158:    Output Parameter:
159: .  y - output vector

161:    Note:
162:    Currently implemented via STApplyTranspose() with appropriate conjugation.

164:    Level: developer

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

176:   STCheckMatrices(st,1);
178:   VecSetErrorIfLocked(y,3);
179:   STGetOperator_Private(st,&Op);
180:   MatMultHermitianTranspose(Op,x,y);
181:   return 0;
182: }

184: /*@
185:    STGetBilinearForm - Returns the matrix used in the bilinear form with a
186:    generalized problem with semi-definite B.

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

190:    Input Parameters:
191: .  st - the spectral transformation context

193:    Output Parameter:
194: .  B - output matrix

196:    Notes:
197:    The output matrix B must be destroyed after use. It will be NULL in
198:    case of standard eigenproblems.

200:    Level: developer

202: .seealso: BVSetMatrix()
203: @*/
204: PetscErrorCode STGetBilinearForm(ST st,Mat *B)
205: {
209:   STCheckMatrices(st,1);
210:   PetscUseTypeMethod(st,getbilinearform,B);
211:   return 0;
212: }

214: PetscErrorCode STGetBilinearForm_Default(ST st,Mat *B)
215: {
216:   if (st->nmat==1) *B = NULL;
217:   else {
218:     *B = st->A[1];
219:     PetscObjectReference((PetscObject)*B);
220:   }
221:   return 0;
222: }

224: static PetscErrorCode MatMult_STOperator(Mat Op,Vec x,Vec y)
225: {
226:   ST             st;

228:   MatShellGetContext(Op,&st);
229:   STSetUp(st);
230:   PetscLogEventBegin(ST_Apply,st,x,y,0);
231:   if (st->D) { /* with balancing */
232:     VecPointwiseDivide(st->wb,x,st->D);
233:     PetscUseTypeMethod(st,apply,st->wb,y);
234:     VecPointwiseMult(y,y,st->D);
235:   } else PetscUseTypeMethod(st,apply,x,y);
236:   PetscLogEventEnd(ST_Apply,st,x,y,0);
237:   return 0;
238: }

240: static PetscErrorCode MatMultTranspose_STOperator(Mat Op,Vec x,Vec y)
241: {
242:   ST             st;

244:   MatShellGetContext(Op,&st);
245:   STSetUp(st);
246:   PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);
247:   if (st->D) { /* with balancing */
248:     VecPointwiseMult(st->wb,x,st->D);
249:     PetscUseTypeMethod(st,applytrans,st->wb,y);
250:     VecPointwiseDivide(y,y,st->D);
251:   } else PetscUseTypeMethod(st,applytrans,x,y);
252:   PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);
253:   return 0;
254: }

256: #if defined(PETSC_USE_COMPLEX)
257: static PetscErrorCode MatMultHermitianTranspose_STOperator(Mat Op,Vec x,Vec y)
258: {
259:   ST             st;

261:   MatShellGetContext(Op,&st);
262:   STSetUp(st);
263:   PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);
264:   if (!st->wht) MatCreateVecs(st->A[0],&st->wht,NULL);
265:   VecCopy(x,st->wht);
266:   VecConjugate(st->wht);
267:   if (st->D) { /* with balancing */
268:     VecPointwiseMult(st->wb,st->wht,st->D);
269:     PetscUseTypeMethod(st,applytrans,st->wb,y);
270:     VecPointwiseDivide(y,y,st->D);
271:   } else PetscUseTypeMethod(st,applytrans,st->wht,y);
272:   VecConjugate(y);
273:   PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);
274:   return 0;
275: }
276: #endif

278: static PetscErrorCode MatMatMult_STOperator(Mat Op,Mat B,Mat C,void *ctx)
279: {
280:   ST             st;

282:   MatShellGetContext(Op,&st);
283:   STSetUp(st);
284:   PetscLogEventBegin(ST_Apply,st,B,C,0);
285:   STApplyMat_Generic(st,B,C);
286:   PetscLogEventEnd(ST_Apply,st,B,C,0);
287:   return 0;
288: }

290: PetscErrorCode STGetOperator_Private(ST st,Mat *Op)
291: {
292:   PetscInt       m,n,M,N;
293:   Vec            v;
294:   VecType        vtype;

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

325: /*@
326:    STGetOperator - Returns a shell matrix that represents the operator of the
327:    spectral transformation.

329:    Collective on st

331:    Input Parameter:
332: .  st - the spectral transformation context

334:    Output Parameter:
335: .  Op - operator matrix

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

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

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

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

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

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

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

365:    Level: advanced

367: .seealso: STApply(), STApplyTranspose(), STSetBalanceMatrix(), STShellSetApply(),
368:           STGetKSP(), STSetShift(), STRestoreOperator(), STSetMatrices()
369: @*/
370: PetscErrorCode STGetOperator(ST st,Mat *Op)
371: {
374:   STCheckMatrices(st,1);
375:   STCheckNotSeized(st,1);
377:   STGetOperator_Private(st,Op);
378:   if (Op) st->opseized = PETSC_TRUE;
379:   return 0;
380: }

382: /*@
383:    STRestoreOperator - Restore the previously seized operator matrix.

385:    Collective on st

387:    Input Parameters:
388: +  st - the spectral transformation context
389: -  Op - operator matrix

391:    Notes:
392:    The arguments must match the corresponding call to STGetOperator().

394:    Level: advanced

396: .seealso: STGetOperator()
397: @*/
398: PetscErrorCode STRestoreOperator(ST st,Mat *Op)
399: {
404:   *Op = NULL;
405:   st->opseized = PETSC_FALSE;
406:   return 0;
407: }

409: /*
410:    STComputeOperator - Computes the matrices that constitute the operator

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

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

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

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

426:    STSHELL does not compute anything here, but sets the flag as if it was ready.
427: */
428: PetscErrorCode STComputeOperator(ST st)
429: {
430:   PC             pc;

434:   if (!st->opready && st->ops->computeoperator) {
435:     PetscInfo(st,"Building the operator matrices\n");
436:     STCheckMatrices(st,1);
437:     if (!st->T) PetscCalloc1(PetscMax(2,st->nmat),&st->T);
438:     PetscLogEventBegin(ST_ComputeOperator,st,0,0,0);
439:     PetscUseTypeMethod(st,computeoperator);
440:     PetscLogEventEnd(ST_ComputeOperator,st,0,0,0);
441:     if (st->usesksp) {
442:       if (!st->ksp) STGetKSP(st,&st->ksp);
443:       if (st->P) {
444:         STSetDefaultKSP(st);
445:         ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P);
446:       } else {
447:         /* STPRECOND defaults to PCNONE if st->P is empty */
448:         KSPGetPC(st->ksp,&pc);
449:         PCSetType(pc,PCNONE);
450:       }
451:     }
452:   }
453:   st->opready = PETSC_TRUE;
454:   return 0;
455: }

457: /*@
458:    STSetUp - Prepares for the use of a spectral transformation.

460:    Collective on st

462:    Input Parameter:
463: .  st - the spectral transformation context

465:    Level: advanced

467: .seealso: STCreate(), STApply(), STDestroy()
468: @*/
469: PetscErrorCode STSetUp(ST st)
470: {
471:   PetscInt       i,n,k;

475:   STCheckMatrices(st,1);
476:   switch (st->state) {
477:     case ST_STATE_INITIAL:
478:       PetscInfo(st,"Setting up new ST\n");
479:       if (!((PetscObject)st)->type_name) STSetType(st,STSHIFT);
480:       break;
481:     case ST_STATE_SETUP:
482:       return 0;
483:     case ST_STATE_UPDATED:
484:       PetscInfo(st,"Setting up updated ST\n");
485:       break;
486:   }
487:   PetscLogEventBegin(ST_SetUp,st,0,0,0);
488:   if (st->state!=ST_STATE_UPDATED) {
489:     if (!(st->nmat<3 && st->opready)) {
490:       if (st->T) {
491:         for (i=0;i<PetscMax(2,st->nmat);i++) MatDestroy(&st->T[i]);
492:       }
493:       MatDestroy(&st->P);
494:     }
495:   }
496:   if (st->D) {
497:     MatGetLocalSize(st->A[0],NULL,&n);
498:     VecGetLocalSize(st->D,&k);
500:     if (!st->wb) VecDuplicate(st->D,&st->wb);
501:   }
502:   if (st->nmat<3 && st->transform) STComputeOperator(st);
503:   else {
504:     if (!st->T) PetscCalloc1(PetscMax(2,st->nmat),&st->T);
505:   }
506:   PetscTryTypeMethod(st,setup);
507:   st->state = ST_STATE_SETUP;
508:   PetscLogEventEnd(ST_SetUp,st,0,0,0);
509:   return 0;
510: }

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

516:    alpha - value of the parameter of the transformed polynomial
517:    beta - value of the previous shift (only used in inplace mode)
518:    k - index of first matrix included in the computation
519:    coeffs - coefficients of the expansion
520:    initial - true if this is the first time
521:    precond - whether the preconditioner matrix must be computed
522: */
523: PetscErrorCode STMatMAXPY_Private(ST st,PetscScalar alpha,PetscScalar beta,PetscInt k,PetscScalar *coeffs,PetscBool initial,PetscBool precond,Mat *S)
524: {
525:   PetscInt       *matIdx=NULL,nmat,i,ini=-1;
526:   PetscScalar    t=1.0,ta,gamma;
527:   PetscBool      nz=PETSC_FALSE;
528:   Mat            *A=precond?st->Psplit:st->A;
529:   MatStructure   str=precond?st->strp:st->str;

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

596: /*
597:    Computes the values of the coefficients required by STMatMAXPY_Private
598:    for the case of monomial basis.
599: */
600: PetscErrorCode STCoeffs_Monomial(ST st, PetscScalar *coeffs)
601: {
602:   PetscInt  k,i,ini,inip;

604:   /* Compute binomial coefficients */
605:   ini = (st->nmat*(st->nmat-1))/2;
606:   for (i=0;i<st->nmat;i++) coeffs[ini+i]=1.0;
607:   for (k=st->nmat-1;k>=1;k--) {
608:     inip = ini+1;
609:     ini = (k*(k-1))/2;
610:     coeffs[ini] = 1.0;
611:     for (i=1;i<k;i++) coeffs[ini+i] = coeffs[ini+i-1]+coeffs[inip+i-1];
612:   }
613:   return 0;
614: }

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

620:    Collective on st

622:    Input Parameters:
623: .  st  - the spectral transformation context

625:    Level: developer

627: .seealso: EPSSolve()
628: @*/
629: PetscErrorCode STPostSolve(ST st)
630: {
633:   PetscTryTypeMethod(st,postsolve);
634:   return 0;
635: }

637: /*@
638:    STBackTransform - Back-transformation phase, intended for
639:    spectral transformations which require to transform the computed
640:    eigenvalues back to the original eigenvalue problem.

642:    Not Collective

644:    Input Parameters:
645: +  st   - the spectral transformation context
646: .  n    - number of eigenvalues
647: .  eigr - real part of a computed eigenvalues
648: -  eigi - imaginary part of a computed eigenvalues

650:    Level: developer

652: .seealso: STIsInjective()
653: @*/
654: PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
655: {
658:   PetscTryTypeMethod(st,backtransform,n,eigr,eigi);
659:   return 0;
660: }

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

667:    Not collective

669:    Input Parameter:
670: .  st   - the spectral transformation context

672:    Output Parameter:
673: .  is - the answer

675:    Level: developer

677: .seealso: STBackTransform()
678: @*/
679: PetscErrorCode STIsInjective(ST st,PetscBool* is)
680: {
681:   PetscBool      shell;


687:   PetscObjectTypeCompare((PetscObject)st,STSHELL,&shell);
688:   if (shell) STIsInjective_Shell(st,is);
689:   else *is = st->ops->backtransform? PETSC_TRUE: PETSC_FALSE;
690:   return 0;
691: }

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

696:    Collective on st

698:    Input Parameters:
699: +  st     - the spectral transformation context
700: .  sigma  - the shift
701: -  coeffs - the coefficients (may be NULL)

703:    Note:
704:    This function is not intended to be called by end users, but by SLEPc
705:    solvers that use ST. It builds matrix st->P as follows, then calls KSPSetUp().
706: .vb
707:     If (coeffs)  st->P = Sum_{i=0..nmat-1} coeffs[i]*sigma^i*A_i
708:     else         st->P = Sum_{i=0..nmat-1} sigma^i*A_i
709: .ve

711:    Level: developer

713: .seealso: STMatSolve()
714: @*/
715: PetscErrorCode STMatSetUp(ST st,PetscScalar sigma,PetscScalar *coeffs)
716: {
719:   STCheckMatrices(st,1);

721:   PetscLogEventBegin(ST_MatSetUp,st,0,0,0);
722:   STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_FALSE,&st->P);
723:   if (st->Psplit) STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,PETSC_TRUE,&st->Pmat);
724:   ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P);
725:   KSPSetUp(st->ksp);
726:   PetscLogEventEnd(ST_MatSetUp,st,0,0,0);
727:   return 0;
728: }

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

733:    Collective on st

735:    Input Parameters:
736: +  st - the spectral transformation context
737: -  nw - number of work vectors to allocate

739:    Developer Notes:
740:    This is SLEPC_EXTERN because it may be required by shell STs.

742:    Level: developer

744: .seealso: STMatCreateVecs()
745: @*/
746: PetscErrorCode STSetWorkVecs(ST st,PetscInt nw)
747: {
748:   PetscInt       i;

753:   if (st->nwork < nw) {
754:     VecDestroyVecs(st->nwork,&st->work);
755:     st->nwork = nw;
756:     PetscMalloc1(nw,&st->work);
757:     for (i=0;i<nw;i++) STMatCreateVecs(st,&st->work[i],NULL);
758:   }
759:   return 0;
760: }