Actual source code: stsolve.c
slepc-3.18.2 2023-01-26
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: }