Actual source code: stsolve.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcst.h" I*/
16: /*@
17: STApply - Applies the spectral transformation operator to a vector, for
18: instance (A - sB)^-1 B in the case of the shift-and-invert transformation
19: and generalized eigenproblem.
21: Collective on ST and Vec
23: Input Parameters:
24: + st - the spectral transformation context
25: - x - input vector
27: Output Parameter:
28: . y - output vector
30: Level: developer
32: .seealso: STApplyTranspose()
33: @*/
34: PetscErrorCode STApply(ST st,Vec x,Vec y)
35: {
43: STCheckMatrices(st,1);
44: if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
45: VecSetErrorIfLocked(y,3);
47: if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
49: if (!st->ops->apply) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST does not have apply");
50: VecLockReadPush(x);
51: PetscLogEventBegin(ST_Apply,st,x,y,0);
52: if (st->D) { /* with balancing */
53: VecPointwiseDivide(st->wb,x,st->D);
54: (*st->ops->apply)(st,st->wb,y);
55: VecPointwiseMult(y,y,st->D);
56: } else {
57: (*st->ops->apply)(st,x,y);
58: }
59: PetscLogEventEnd(ST_Apply,st,x,y,0);
60: VecLockReadPop(x);
61: return(0);
62: }
64: /*@
65: STApplyTranspose - Applies the transpose of the operator to a vector, for
66: instance B^T(A - sB)^-T in the case of the shift-and-invert transformation
67: and generalized eigenproblem.
69: Collective on ST and Vec
71: Input Parameters:
72: + st - the spectral transformation context
73: - x - input vector
75: Output Parameter:
76: . y - output vector
78: Level: developer
80: .seealso: STApply()
81: @*/
82: PetscErrorCode STApplyTranspose(ST st,Vec x,Vec y)
83: {
91: STCheckMatrices(st,1);
92: if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
93: VecSetErrorIfLocked(y,3);
95: if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
97: if (!st->ops->applytrans) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST does not have applytrans");
98: VecLockReadPush(x);
99: PetscLogEventBegin(ST_ApplyTranspose,st,x,y,0);
100: if (st->D) { /* with balancing */
101: VecPointwiseMult(st->wb,x,st->D);
102: (*st->ops->applytrans)(st,st->wb,y);
103: VecPointwiseDivide(y,y,st->D);
104: } else {
105: (*st->ops->applytrans)(st,x,y);
106: }
107: PetscLogEventEnd(ST_ApplyTranspose,st,x,y,0);
108: VecLockReadPop(x);
109: return(0);
110: }
112: /*@
113: STGetBilinearForm - Returns the matrix used in the bilinear form with a
114: generalized problem with semi-definite B.
116: Not collective, though a parallel Mat may be returned
118: Input Parameters:
119: . st - the spectral transformation context
121: Output Parameter:
122: . B - output matrix
124: Notes:
125: The output matrix B must be destroyed after use. It will be NULL in
126: case of standard eigenproblems.
128: Level: developer
129: @*/
130: PetscErrorCode STGetBilinearForm(ST st,Mat *B)
131: {
138: STCheckMatrices(st,1);
139: (*st->ops->getbilinearform)(st,B);
140: return(0);
141: }
143: PetscErrorCode STGetBilinearForm_Default(ST st,Mat *B)
144: {
148: if (st->nmat==1) *B = NULL;
149: else {
150: *B = st->A[1];
151: PetscObjectReference((PetscObject)*B);
152: }
153: return(0);
154: }
156: static PetscErrorCode MatMult_STOperator(Mat Op,Vec x,Vec y)
157: {
159: ST st;
162: MatShellGetContext(Op,(void**)&st);
163: STApply(st,x,y);
164: return(0);
165: }
167: static PetscErrorCode MatMultTranspose_STOperator(Mat Op,Vec x,Vec y)
168: {
170: ST st;
173: MatShellGetContext(Op,(void**)&st);
174: STApplyTranspose(st,x,y);
175: return(0);
176: }
178: /*@
179: STGetOperator - Returns a shell matrix that represents the spectral transformation.
181: Collective on ST
183: Input Parameters:
184: . st - the spectral transformation context
186: Output Parameter:
187: . Op - output matrix
189: Notes:
190: The returned shell matrix is essentially a wrapper to the STApply() and
191: STApplyTranspose() operations. It must be destroyed after use.
193: Level: advanced
195: .seealso: STApply(), STApplyTranspose()
196: @*/
197: PetscErrorCode STGetOperator(ST st,Mat *Op)
198: {
200: PetscInt m,n,M,N;
206: STCheckMatrices(st,1);
208: MatGetLocalSize(st->A[0],&m,&n);
209: MatGetSize(st->A[0],&M,&N);
210: MatCreateShell(PetscObjectComm((PetscObject)st),m,n,M,N,st,Op);
211: MatShellSetOperation(*Op,MATOP_MULT,(void(*)(void))MatMult_STOperator);
212: MatShellSetOperation(*Op,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_STOperator);
213: return(0);
214: }
216: /*@
217: STSetUp - Prepares for the use of a spectral transformation.
219: Collective on ST
221: Input Parameter:
222: . st - the spectral transformation context
224: Level: advanced
226: .seealso: STCreate(), STApply(), STDestroy()
227: @*/
228: PetscErrorCode STSetUp(ST st)
229: {
230: PetscInt i,n,k;
236: STCheckMatrices(st,1);
237: switch (st->state) {
238: case ST_STATE_INITIAL:
239: PetscInfo(st,"Setting up new ST\n");
240: if (!((PetscObject)st)->type_name) {
241: STSetType(st,STSHIFT);
242: }
243: break;
244: case ST_STATE_SETUP:
245: return(0);
246: case ST_STATE_UPDATED:
247: PetscInfo(st,"Setting up updated ST\n");
248: break;
249: }
250: PetscLogEventBegin(ST_SetUp,st,0,0,0);
251: if (!st->T) {
252: PetscMalloc1(PetscMax(2,st->nmat),&st->T);
253: PetscLogObjectMemory((PetscObject)st,PetscMax(2,st->nmat)*sizeof(Mat));
254: for (i=0;i<PetscMax(2,st->nmat);i++) st->T[i] = NULL;
255: } else if (st->state!=ST_STATE_UPDATED) {
256: for (i=0;i<PetscMax(2,st->nmat);i++) {
257: MatDestroy(&st->T[i]);
258: }
259: }
260: if (st->state!=ST_STATE_UPDATED) { MatDestroy(&st->P); }
261: if (st->D) {
262: MatGetLocalSize(st->A[0],NULL,&n);
263: VecGetLocalSize(st->D,&k);
264: if (n != k) SETERRQ2(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_SIZ,"Balance matrix has wrong dimension %D (should be %D)",k,n);
265: if (!st->wb) {
266: VecDuplicate(st->D,&st->wb);
267: PetscLogObjectParent((PetscObject)st,(PetscObject)st->wb);
268: }
269: }
270: STSetDefaultKSP(st);
271: if (st->ops->setup) { (*st->ops->setup)(st); }
272: st->state = ST_STATE_SETUP;
273: PetscLogEventEnd(ST_SetUp,st,0,0,0);
274: return(0);
275: }
277: /*
278: Computes coefficients for the transformed polynomial,
279: and stores the result in argument S.
281: alpha - value of the parameter of the transformed polynomial
282: beta - value of the previous shift (only used in inplace mode)
283: k - number of A matrices involved in the computation
284: coeffs - coefficients of the expansion
285: initial - true if this is the first time (only relevant for shell mode)
286: */
287: PetscErrorCode STMatMAXPY_Private(ST st,PetscScalar alpha,PetscScalar beta,PetscInt k,PetscScalar *coeffs,PetscBool initial,Mat *S)
288: {
290: PetscInt *matIdx=NULL,nmat,i,ini=-1;
291: PetscScalar t=1.0,ta,gamma;
292: PetscBool nz=PETSC_FALSE;
295: nmat = st->nmat-k;
296: switch (st->shift_matrix) {
297: case ST_MATMODE_INPLACE:
298: if (st->nmat>2) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"ST_MATMODE_INPLACE not supported for polynomial eigenproblems");
299: if (initial) {
300: PetscObjectReference((PetscObject)st->A[0]);
301: *S = st->A[0];
302: gamma = alpha;
303: } else gamma = alpha-beta;
304: if (gamma != 0.0) {
305: if (st->nmat>1) {
306: MatAXPY(*S,gamma,st->A[1],st->str);
307: } else {
308: MatShift(*S,gamma);
309: }
310: }
311: break;
312: case ST_MATMODE_SHELL:
313: if (initial) {
314: if (st->nmat>2) {
315: PetscMalloc1(nmat,&matIdx);
316: for (i=0;i<nmat;i++) matIdx[i] = k+i;
317: }
318: STMatShellCreate(st,alpha,nmat,matIdx,coeffs,S);
319: PetscLogObjectParent((PetscObject)st,(PetscObject)*S);
320: if (st->nmat>2) { PetscFree(matIdx); }
321: } else {
322: STMatShellShift(*S,alpha);
323: }
324: break;
325: case ST_MATMODE_COPY:
326: if (coeffs) {
327: for (i=0;i<nmat && ini==-1;i++) {
328: if (coeffs[i]!=0.0) ini = i;
329: else t *= alpha;
330: }
331: if (coeffs[ini] != 1.0) nz = PETSC_TRUE;
332: for (i=ini+1;i<nmat&&!nz;i++) if (coeffs[i]!=0.0) nz = PETSC_TRUE;
333: } else { nz = PETSC_TRUE; ini = 0; }
334: if ((alpha == 0.0 || !nz) && t==1.0) {
335: PetscObjectReference((PetscObject)st->A[k+ini]);
336: MatDestroy(S);
337: *S = st->A[k+ini];
338: } else {
339: if (*S && *S!=st->A[k+ini]) {
340: MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
341: MatCopy(st->A[k+ini],*S,DIFFERENT_NONZERO_PATTERN);
342: } else {
343: MatDestroy(S);
344: MatDuplicate(st->A[k+ini],MAT_COPY_VALUES,S);
345: MatSetOption(*S,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
346: PetscLogObjectParent((PetscObject)st,(PetscObject)*S);
347: }
348: if (coeffs && coeffs[ini]!=1.0) {
349: MatScale(*S,coeffs[ini]);
350: }
351: for (i=ini+k+1;i<PetscMax(2,st->nmat);i++) {
352: t *= alpha;
353: ta = t;
354: if (coeffs) ta *= coeffs[i-k];
355: if (ta!=0.0) {
356: if (st->nmat>1) {
357: MatAXPY(*S,ta,st->A[i],st->str);
358: } else {
359: MatShift(*S,ta);
360: }
361: }
362: }
363: }
364: }
365: STMatSetHermitian(st,*S);
366: return(0);
367: }
369: /*
370: Computes the values of the coefficients required by STMatMAXPY_Private
371: for the case of monomial basis.
372: */
373: PetscErrorCode STCoeffs_Monomial(ST st, PetscScalar *coeffs)
374: {
375: PetscInt k,i,ini,inip;
378: /* Compute binomial coefficients */
379: ini = (st->nmat*(st->nmat-1))/2;
380: for (i=0;i<st->nmat;i++) coeffs[ini+i]=1.0;
381: for (k=st->nmat-1;k>=1;k--) {
382: inip = ini+1;
383: ini = (k*(k-1))/2;
384: coeffs[ini] = 1.0;
385: for (i=1;i<k;i++) coeffs[ini+i] = coeffs[ini+i-1]+coeffs[inip+i-1];
386: }
387: return(0);
388: }
390: /*@
391: STPostSolve - Optional post-solve phase, intended for any actions that must
392: be performed on the ST object after the eigensolver has finished.
394: Collective on ST
396: Input Parameters:
397: . st - the spectral transformation context
399: Level: developer
401: .seealso: EPSSolve()
402: @*/
403: PetscErrorCode STPostSolve(ST st)
404: {
410: if (st->ops->postsolve) {
411: (*st->ops->postsolve)(st);
412: }
413: return(0);
414: }
416: /*@
417: STBackTransform - Back-transformation phase, intended for
418: spectral transformations which require to transform the computed
419: eigenvalues back to the original eigenvalue problem.
421: Not Collective
423: Input Parameters:
424: + st - the spectral transformation context
425: eigr - real part of a computed eigenvalue
426: - eigi - imaginary part of a computed eigenvalue
428: Level: developer
430: .seealso: STIsInjective()
431: @*/
432: PetscErrorCode STBackTransform(ST st,PetscInt n,PetscScalar* eigr,PetscScalar* eigi)
433: {
439: if (st->ops->backtransform) {
440: (*st->ops->backtransform)(st,n,eigr,eigi);
441: }
442: return(0);
443: }
445: /*@
446: STIsInjective - Ask if this spectral transformation is injective or not
447: (that is, if it corresponds to a one-to-one mapping). If not, then it
448: does not make sense to call STBackTransform().
450: Not collective
452: Input Parameter:
453: . st - the spectral transformation context
455: Output Parameter:
456: . is - the answer
458: Level: developer
460: .seealso: STBackTransform()
461: @*/
462: PetscErrorCode STIsInjective(ST st,PetscBool* is)
463: {
465: PetscBool shell;
472: PetscObjectTypeCompare((PetscObject)st,STSHELL,&shell);
473: if (shell) {
474: STIsInjective_Shell(st,is);
475: } else *is = st->ops->backtransform? PETSC_TRUE: PETSC_FALSE;
476: return(0);
477: }
479: /*@
480: STMatSetUp - Build the preconditioner matrix used in STMatSolve().
482: Collective on ST
484: Input Parameters:
485: + st - the spectral transformation context
486: . sigma - the shift
487: - coeffs - the coefficients (may be NULL)
489: Note:
490: This function is not intended to be called by end users, but by SLEPc
491: solvers that use ST. It builds matrix st->P as follows, then calls KSPSetUp().
492: .vb
493: If (coeffs): st->P = Sum_{i=0:nmat-1} coeffs[i]*sigma^i*A_i.
494: else st->P = Sum_{i=0:nmat-1} sigma^i*A_i
495: .ve
497: Level: developer
499: .seealso: STMatSolve()
500: @*/
501: PetscErrorCode STMatSetUp(ST st,PetscScalar sigma,PetscScalar *coeffs)
502: {
508: STCheckMatrices(st,1);
510: PetscLogEventBegin(ST_MatSetUp,st,0,0,0);
511: STMatMAXPY_Private(st,sigma,0.0,0,coeffs,PETSC_TRUE,&st->P);
512: if (!st->ksp) { STGetKSP(st,&st->ksp); }
513: STCheckFactorPackage(st);
514: KSPSetOperators(st->ksp,st->P,st->P);
515: KSPSetUp(st->ksp);
516: PetscLogEventEnd(ST_MatSetUp,st,0,0,0);
517: return(0);
518: }
520: /*@
521: STSetWorkVecs - Sets a number of work vectors into the ST object.
523: Collective on ST
525: Input Parameters:
526: + st - the spectral transformation context
527: - nw - number of work vectors to allocate
529: Developers Note:
530: This is SLEPC_EXTERN because it may be required by shell STs.
532: Level: developer
533: @*/
534: PetscErrorCode STSetWorkVecs(ST st,PetscInt nw)
535: {
537: PetscInt i;
542: if (nw <= 0) SETERRQ1(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %D",nw);
543: if (st->nwork < nw) {
544: VecDestroyVecs(st->nwork,&st->work);
545: st->nwork = nw;
546: PetscMalloc1(nw,&st->work);
547: for (i=0;i<nw;i++) { STMatCreateVecs(st,&st->work[i],NULL); }
548: PetscLogObjectParents(st,nw,st->work);
549: }
550: return(0);
551: }