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: }