Actual source code: stsolve.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }