Actual source code: stsles.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 related to the KSP object associated with it
 12: */

 14: #include <slepc/private/stimpl.h>            /*I "slepcst.h" I*/

 16: /*
 17:    This is used to set a default type for the KSP and PC objects.
 18:    It is called at STSetFromOptions (before KSPSetFromOptions)
 19:    and also at STSetUp (in case STSetFromOptions was not called).
 20: */
 21: PetscErrorCode STSetDefaultKSP(ST st)
 22: {

 28:   if (!st->ksp) { STGetKSP(st,&st->ksp); }
 29:   if (st->ops->setdefaultksp) { (*st->ops->setdefaultksp)(st); }
 30:   return(0);
 31: }

 33: /*
 34:    This is done by all ST types except PRECOND.
 35:    The default is an LU direct solver, or GMRES+Jacobi if matmode=shell.
 36: */
 37: PetscErrorCode STSetDefaultKSP_Default(ST st)
 38: {
 40:   PC             pc;
 41:   PCType         pctype;
 42:   KSPType        ksptype;

 45:   KSPGetPC(st->ksp,&pc);
 46:   KSPGetType(st->ksp,&ksptype);
 47:   PCGetType(pc,&pctype);
 48:   if (!pctype && !ksptype) {
 49:     if (st->shift_matrix == ST_MATMODE_SHELL) {
 50:       KSPSetType(st->ksp,KSPGMRES);
 51:       PCSetType(pc,PCJACOBI);
 52:     } else {
 53:       KSPSetType(st->ksp,KSPPREONLY);
 54:       PCSetType(pc,PCLU);
 55:     }
 56:   }
 57:   KSPSetErrorIfNotConverged(st->ksp,PETSC_TRUE);
 58:   return(0);
 59: }

 61: /*@
 62:    STMatMult - Computes the matrix-vector product y = T[k] x, where T[k] is
 63:    the k-th matrix of the spectral transformation.

 65:    Collective on ST

 67:    Input Parameters:
 68: +  st - the spectral transformation context
 69: .  k  - index of matrix to use
 70: -  x  - the vector to be multiplied

 72:    Output Parameter:
 73: .  y - the result

 75:    Level: developer

 77: .seealso: STMatMultTranspose()
 78: @*/
 79: PetscErrorCode STMatMult(ST st,PetscInt k,Vec x,Vec y)
 80: {

 88:   STCheckMatrices(st,1);
 89:   if (k<0 || k>=PetscMax(2,st->nmat)) SETERRQ1(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",st->nmat);
 90:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
 91:   VecSetErrorIfLocked(y,3);

 93:   if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
 94:   VecLockReadPush(x);
 95:   PetscLogEventBegin(ST_MatMult,st,x,y,0);
 96:   if (!st->T[k]) {
 97:     /* T[k]=NULL means identity matrix */
 98:     VecCopy(x,y);
 99:   } else {
100:     MatMult(st->T[k],x,y);
101:   }
102:   PetscLogEventEnd(ST_MatMult,st,x,y,0);
103:   VecLockReadPop(x);
104:   return(0);
105: }

107: /*@
108:    STMatMultTranspose - Computes the matrix-vector product y = T[k]' x, where T[k] is
109:    the k-th matrix of the spectral transformation.

111:    Collective on ST

113:    Input Parameters:
114: +  st - the spectral transformation context
115: .  k  - index of matrix to use
116: -  x  - the vector to be multiplied

118:    Output Parameter:
119: .  y - the result

121:    Level: developer

123: .seealso: STMatMult()
124: @*/
125: PetscErrorCode STMatMultTranspose(ST st,PetscInt k,Vec x,Vec y)
126: {

134:   STCheckMatrices(st,1);
135:   if (k<0 || k>=PetscMax(2,st->nmat)) SETERRQ1(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",st->nmat);
136:   if (x == y) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
137:   VecSetErrorIfLocked(y,3);

139:   if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
140:   VecLockReadPush(x);
141:   PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0);
142:   if (!st->T[k]) {
143:     /* T[k]=NULL means identity matrix */
144:     VecCopy(x,y);
145:   } else {
146:     MatMultTranspose(st->T[k],x,y);
147:   }
148:   PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0);
149:   VecLockReadPop(x);
150:   return(0);
151: }

153: /*@
154:    STMatSolve - Solves P x = b, where P is the preconditioner matrix of
155:    the spectral transformation, using a KSP object stored internally.

157:    Collective on ST

159:    Input Parameters:
160: +  st - the spectral transformation context
161: -  b  - right hand side vector

163:    Output Parameter:
164: .  x - computed solution

166:    Level: developer

168: .seealso: STMatSolveTranspose()
169: @*/
170: PetscErrorCode STMatSolve(ST st,Vec b,Vec x)
171: {
173:   PetscInt       its;
174:   PetscBool      flg;

180:   STCheckMatrices(st,1);
181:   if (x == b) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
182:   VecSetErrorIfLocked(x,3);

184:   if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
185:   VecLockReadPush(b);
186:   PetscLogEventBegin(ST_MatSolve,st,b,x,0);
187:   PetscObjectTypeCompareAny((PetscObject)st,&flg,STPRECOND,STSHELL,"");
188:   if (!flg && !st->P) {
189:     /* P=NULL means identity matrix */
190:     VecCopy(b,x);
191:   } else {
192:     if (!st->ksp) { STGetKSP(st,&st->ksp); }
193:     KSPSolve(st->ksp,b,x);
194:     KSPGetIterationNumber(st->ksp,&its);
195:     PetscInfo1(st,"Linear solve iterations=%D\n",its);
196:   }
197:   PetscLogEventEnd(ST_MatSolve,st,b,x,0);
198:   VecLockReadPop(b);
199:   return(0);
200: }

202: /*@
203:    STMatSolveTranspose - Solves P' x = b, where P is the preconditioner matrix of
204:    the spectral transformation, using a KSP object stored internally.

206:    Collective on ST

208:    Input Parameters:
209: .  st - the spectral transformation context
210: .  b  - right hand side vector

212:    Output Parameter:
213: .  x - computed solution

215:    Level: developer

217: .seealso: STMatSolve()
218: @*/
219: PetscErrorCode STMatSolveTranspose(ST st,Vec b,Vec x)
220: {
222:   PetscInt       its;
223:   PetscBool      flg;

229:   STCheckMatrices(st,1);
230:   if (x == b) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
231:   VecSetErrorIfLocked(x,3);

233:   if (st->state!=ST_STATE_SETUP) { STSetUp(st); }
234:   VecLockReadPush(b);
235:   PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0);
236:   PetscObjectTypeCompareAny((PetscObject)st,&flg,STPRECOND,STSHELL,"");
237:   if (!flg && !st->P) {
238:     /* P=NULL means identity matrix */
239:     VecCopy(b,x);
240:   } else {
241:     if (!st->ksp) { STGetKSP(st,&st->ksp); }
242:     KSPSolveTranspose(st->ksp,b,x);
243:     KSPGetIterationNumber(st->ksp,&its);
244:     PetscInfo1(st,"Linear solve iterations=%D\n",its);
245:   }
246:   PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0);
247:   VecLockReadPop(b);
248:   return(0);
249: }

251: /*
252:    STMatSetHermitian - Sets the Hermitian flag to the ST matrix.

254:    Input Parameters:
255: .  st - the spectral transformation context
256: .  M  - matrix
257: */
258: PetscErrorCode STMatSetHermitian(ST st,Mat M)
259: {
261:   PetscBool      set,aherm,mherm;
262:   PetscInt       i;

265: #if defined(PETSC_USE_COMPLEX)
266:   if (PetscImaginaryPart(st->sigma)!=0.0) mherm = PETSC_FALSE;
267:   else {
268:     mherm = PETSC_TRUE;
269:     for (i=0;i<st->nmat;i++) {
270:       MatIsHermitianKnown(st->A[i],&set,&aherm);
271:       if (!set) aherm = PETSC_FALSE;
272:       if (!aherm) { mherm = PETSC_FALSE; break; }
273:       if (PetscRealPart(st->sigma)==0.0) break;
274:     }
275:   }
276:   MatSetOption(M,MAT_HERMITIAN,mherm);
277: #else
278:   mherm = PETSC_TRUE;
279:   for (i=0;i<st->nmat;i++) {
280:     MatIsSymmetricKnown(st->A[i],&set,&aherm);
281:     if (!set) aherm = PETSC_FALSE;
282:     if (!aherm) { mherm = PETSC_FALSE; break; }
283:   }
284:   MatSetOption(M,MAT_SYMMETRIC,mherm);
285: #endif
286:   return(0);
287: }

289: PetscErrorCode STCheckFactorPackage(ST st)
290: {
292:   PC             pc;
293:   PetscMPIInt    size;
294:   PetscBool      flg;
295:   MatSolverType  stype;

298:   MPI_Comm_size(PetscObjectComm((PetscObject)st),&size);
299:   if (size==1) return(0);
300:   KSPGetPC(st->ksp,&pc);
301:   PCFactorGetMatSolverType(pc,&stype);
302:   if (stype) {   /* currently selected PC is a factorization */
303:     PetscStrcmp(stype,MATSOLVERPETSC,&flg);
304:     if (flg) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"You chose to solve linear systems with a factorization, but in parallel runs you need to select an external package; see the users guide for details");
305:   }
306:   return(0);
307: }

309: /*@
310:    STSetKSP - Sets the KSP object associated with the spectral
311:    transformation.

313:    Collective on ST

315:    Input Parameters:
316: +  st   - the spectral transformation context
317: -  ksp  - the linear system context

319:    Level: advanced
320: @*/
321: PetscErrorCode STSetKSP(ST st,KSP ksp)
322: {

329:   PetscObjectReference((PetscObject)ksp);
330:   KSPDestroy(&st->ksp);
331:   st->ksp = ksp;
332:   PetscLogObjectParent((PetscObject)st,(PetscObject)st->ksp);
333:   return(0);
334: }

336: /*@
337:    STGetKSP - Gets the KSP object associated with the spectral
338:    transformation.

340:    Not Collective

342:    Input Parameter:
343: .  st - the spectral transformation context

345:    Output Parameter:
346: .  ksp  - the linear system context

348:    Level: intermediate
349: @*/
350: PetscErrorCode STGetKSP(ST st,KSP* ksp)
351: {

357:   if (!st->ksp) {
358:     KSPCreate(PetscObjectComm((PetscObject)st),&st->ksp);
359:     PetscObjectIncrementTabLevel((PetscObject)st->ksp,(PetscObject)st,1);
360:     KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix);
361:     KSPAppendOptionsPrefix(st->ksp,"st_");
362:     PetscLogObjectParent((PetscObject)st,(PetscObject)st->ksp);
363:     PetscObjectSetOptions((PetscObject)st->ksp,((PetscObject)st)->options);
364:     KSPSetTolerances(st->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
365:   }
366:   *ksp = st->ksp;
367:   return(0);
368: }

370: PetscErrorCode STCheckNullSpace_Default(ST st,BV V)
371: {
373:   PetscInt       nc,i,c;
374:   PetscReal      norm;
375:   Vec            *T,w,vi;
376:   Mat            A;
377:   PC             pc;
378:   MatNullSpace   nullsp;

381:   BVGetNumConstraints(V,&nc);
382:   PetscMalloc1(nc,&T);
383:   if (!st->ksp) { STGetKSP(st,&st->ksp); }
384:   KSPGetPC(st->ksp,&pc);
385:   PCGetOperators(pc,&A,NULL);
386:   MatCreateVecs(A,NULL,&w);
387:   c = 0;
388:   for (i=0;i<nc;i++) {
389:     BVGetColumn(V,-nc+i,&vi);
390:     MatMult(A,vi,w);
391:     VecNorm(w,NORM_2,&norm);
392:     if (norm < 1e-8) {
393:       PetscInfo2(st,"Vector %D norm=%g\n",i,(double)norm);
394:       BVCreateVec(V,T+c);
395:       VecCopy(vi,T[c]);
396:       c++;
397:     }
398:     BVRestoreColumn(V,-nc+i,&vi);
399:   }
400:   VecDestroy(&w);
401:   if (c>0) {
402:     MatNullSpaceCreate(PetscObjectComm((PetscObject)st),PETSC_FALSE,c,T,&nullsp);
403:     MatSetNullSpace(A,nullsp);
404:     MatNullSpaceDestroy(&nullsp);
405:     VecDestroyVecs(c,&T);
406:   } else {
407:     PetscFree(T);
408:   }
409:   return(0);
410: }

412: /*@
413:    STCheckNullSpace - Given a basis vectors object, this function tests each
414:    of its constraint vectors to be a nullspace vector of the coefficient
415:    matrix of the associated KSP object. All these nullspace vectors are passed
416:    to the KSP object.

418:    Collective on ST

420:    Input Parameters:
421: +  st - the spectral transformation context
422: -  V  - basis vectors to be checked

424:    Note:
425:    This function allows to handle singular pencils and to solve some problems
426:    in which the nullspace is important (see the users guide for details).

428:    Level: developer

430: .seealso: EPSSetDeflationSpace()
431: @*/
432: PetscErrorCode STCheckNullSpace(ST st,BV V)
433: {
435:   PetscInt       nc;

442:   if (!st->state) SETERRQ(PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must call STSolve() first");

444:   BVGetNumConstraints(V,&nc);
445:   if (nc && st->ops->checknullspace) {
446:     (*st->ops->checknullspace)(st,V);
447:   }
448:   return(0);
449: }