LCOV - code coverage report
Current view: top level - sys/classes/st/interface - stsles.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 224 228 98.2 %
Date: 2024-03-29 00:34:52 Functions: 14 14 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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             : */
      13             : 
      14             : #include <slepc/private/stimpl.h>            /*I "slepcst.h" I*/
      15             : 
      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        1134 : PetscErrorCode STSetDefaultKSP(ST st)
      22             : {
      23        1134 :   PetscFunctionBegin;
      24        1134 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
      25        1134 :   PetscValidType(st,1);
      26        1134 :   if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
      27        1134 :   PetscTryTypeMethod(st,setdefaultksp);
      28        1134 :   PetscFunctionReturn(PETSC_SUCCESS);
      29             : }
      30             : 
      31             : /*
      32             :    This is done by all ST types except PRECOND.
      33             :    The default is an LU direct solver, or GMRES+Jacobi if matmode=shell.
      34             : */
      35         843 : PetscErrorCode STSetDefaultKSP_Default(ST st)
      36             : {
      37         843 :   PC             pc;
      38         843 :   PCType         pctype;
      39         843 :   KSPType        ksptype;
      40             : 
      41         843 :   PetscFunctionBegin;
      42         843 :   PetscCall(KSPGetPC(st->ksp,&pc));
      43         843 :   PetscCall(KSPGetType(st->ksp,&ksptype));
      44         843 :   PetscCall(PCGetType(pc,&pctype));
      45         843 :   if (!pctype && !ksptype) {
      46         587 :     if (st->Pmat || st->Psplit) {
      47          21 :       PetscCall(KSPSetType(st->ksp,KSPBCGS));
      48          21 :       PetscCall(PCSetType(pc,PCBJACOBI));
      49         566 :     } else if (st->matmode == ST_MATMODE_SHELL) {
      50          51 :       PetscCall(KSPSetType(st->ksp,KSPGMRES));
      51          51 :       PetscCall(PCSetType(pc,PCJACOBI));
      52             :     } else {
      53         515 :       PetscCall(KSPSetType(st->ksp,KSPPREONLY));
      54        1000 :       PetscCall(PCSetType(pc,st->asymm?PCCHOLESKY:PCLU));
      55             :     }
      56             :   }
      57         843 :   PetscCall(KSPSetErrorIfNotConverged(st->ksp,PETSC_TRUE));
      58         843 :   PetscFunctionReturn(PETSC_SUCCESS);
      59             : }
      60             : 
      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.
      64             : 
      65             :    Neighbor-wise Collective
      66             : 
      67             :    Input Parameters:
      68             : +  st - the spectral transformation context
      69             : .  k  - index of matrix to use
      70             : -  x  - the vector to be multiplied
      71             : 
      72             :    Output Parameter:
      73             : .  y - the result
      74             : 
      75             :    Level: developer
      76             : 
      77             : .seealso: STMatMultTranspose(), STMatMultHermitianTranspose()
      78             : @*/
      79       67062 : PetscErrorCode STMatMult(ST st,PetscInt k,Vec x,Vec y)
      80             : {
      81       67062 :   PetscFunctionBegin;
      82       67062 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
      83      268248 :   PetscValidLogicalCollectiveInt(st,k,2);
      84       67062 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
      85       67062 :   PetscValidHeaderSpecific(y,VEC_CLASSID,4);
      86       67062 :   STCheckMatrices(st,1);
      87       67062 :   PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
      88       67062 :   PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
      89       67062 :   PetscCall(VecSetErrorIfLocked(y,3));
      90             : 
      91       67062 :   if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
      92       67062 :   PetscCall(VecLockReadPush(x));
      93       67062 :   PetscCall(PetscLogEventBegin(ST_MatMult,st,x,y,0));
      94       67062 :   if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
      95       67062 :   else PetscCall(MatMult(st->T[k],x,y));
      96       67062 :   PetscCall(PetscLogEventEnd(ST_MatMult,st,x,y,0));
      97       67062 :   PetscCall(VecLockReadPop(x));
      98       67062 :   PetscFunctionReturn(PETSC_SUCCESS);
      99             : }
     100             : 
     101             : /*@
     102             :    STMatMultTranspose - Computes the matrix-vector product y = T[k]^T x, where T[k] is
     103             :    the k-th matrix of the spectral transformation.
     104             : 
     105             :    Neighbor-wise Collective
     106             : 
     107             :    Input Parameters:
     108             : +  st - the spectral transformation context
     109             : .  k  - index of matrix to use
     110             : -  x  - the vector to be multiplied
     111             : 
     112             :    Output Parameter:
     113             : .  y - the result
     114             : 
     115             :    Level: developer
     116             : 
     117             : .seealso: STMatMult(), STMatMultHermitianTranspose()
     118             : @*/
     119           6 : PetscErrorCode STMatMultTranspose(ST st,PetscInt k,Vec x,Vec y)
     120             : {
     121           6 :   PetscFunctionBegin;
     122           6 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     123          24 :   PetscValidLogicalCollectiveInt(st,k,2);
     124           6 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
     125           6 :   PetscValidHeaderSpecific(y,VEC_CLASSID,4);
     126           6 :   STCheckMatrices(st,1);
     127           6 :   PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
     128           6 :   PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
     129           6 :   PetscCall(VecSetErrorIfLocked(y,3));
     130             : 
     131           6 :   if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
     132           6 :   PetscCall(VecLockReadPush(x));
     133           6 :   PetscCall(PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0));
     134           6 :   if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
     135           6 :   else PetscCall(MatMultTranspose(st->T[k],x,y));
     136           6 :   PetscCall(PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0));
     137           6 :   PetscCall(VecLockReadPop(x));
     138           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     139             : }
     140             : 
     141             : /*@
     142             :    STMatMultHermitianTranspose - Computes the matrix-vector product y = T[k]^H x, where T[k] is
     143             :    the k-th matrix of the spectral transformation.
     144             : 
     145             :    Neighbor-wise Collective
     146             : 
     147             :    Input Parameters:
     148             : +  st - the spectral transformation context
     149             : .  k  - index of matrix to use
     150             : -  x  - the vector to be multiplied
     151             : 
     152             :    Output Parameter:
     153             : .  y - the result
     154             : 
     155             :    Level: developer
     156             : 
     157             : .seealso: STMatMult(), STMatMultTranspose()
     158             : @*/
     159           3 : PetscErrorCode STMatMultHermitianTranspose(ST st,PetscInt k,Vec x,Vec y)
     160             : {
     161           3 :   PetscFunctionBegin;
     162           3 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     163          12 :   PetscValidLogicalCollectiveInt(st,k,2);
     164           3 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
     165           3 :   PetscValidHeaderSpecific(y,VEC_CLASSID,4);
     166           3 :   STCheckMatrices(st,1);
     167           3 :   PetscCheck(k>=0 && k<PetscMax(2,st->nmat),PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,st->nmat);
     168           3 :   PetscCheck(x!=y,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
     169           3 :   PetscCall(VecSetErrorIfLocked(y,3));
     170             : 
     171           3 :   if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
     172           3 :   PetscCall(VecLockReadPush(x));
     173           3 :   PetscCall(PetscLogEventBegin(ST_MatMultTranspose,st,x,y,0));
     174           3 :   if (!st->T[k]) PetscCall(VecCopy(x,y)); /* T[k]=NULL means identity matrix */
     175           3 :   else PetscCall(MatMultHermitianTranspose(st->T[k],x,y));
     176           3 :   PetscCall(PetscLogEventEnd(ST_MatMultTranspose,st,x,y,0));
     177           3 :   PetscCall(VecLockReadPop(x));
     178           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     179             : }
     180             : 
     181             : /*@
     182             :    STMatSolve - Solves P x = b, where P is the preconditioner matrix of
     183             :    the spectral transformation, using a KSP object stored internally.
     184             : 
     185             :    Collective
     186             : 
     187             :    Input Parameters:
     188             : +  st - the spectral transformation context
     189             : -  b  - right hand side vector
     190             : 
     191             :    Output Parameter:
     192             : .  x - computed solution
     193             : 
     194             :    Level: developer
     195             : 
     196             : .seealso: STMatSolveTranspose(), STMatSolveHermitianTranspose(), STMatMatSolve()
     197             : @*/
     198       53349 : PetscErrorCode STMatSolve(ST st,Vec b,Vec x)
     199             : {
     200       53349 :   PetscFunctionBegin;
     201       53349 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     202       53349 :   PetscValidHeaderSpecific(b,VEC_CLASSID,2);
     203       53349 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
     204       53349 :   STCheckMatrices(st,1);
     205       53349 :   PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
     206       53349 :   PetscCall(VecSetErrorIfLocked(x,3));
     207             : 
     208       53349 :   if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
     209       53349 :   PetscCall(VecLockReadPush(b));
     210       53349 :   PetscCall(PetscLogEventBegin(ST_MatSolve,st,b,x,0));
     211       53349 :   if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
     212       53341 :   else PetscCall(KSPSolve(st->ksp,b,x));
     213       53349 :   PetscCall(PetscLogEventEnd(ST_MatSolve,st,b,x,0));
     214       53349 :   PetscCall(VecLockReadPop(b));
     215       53349 :   PetscFunctionReturn(PETSC_SUCCESS);
     216             : }
     217             : 
     218             : /*@
     219             :    STMatMatSolve - Solves P X = B, where P is the preconditioner matrix of
     220             :    the spectral transformation, using a KSP object stored internally.
     221             : 
     222             :    Collective
     223             : 
     224             :    Input Parameters:
     225             : +  st - the spectral transformation context
     226             : -  B  - right hand side vectors
     227             : 
     228             :    Output Parameter:
     229             : .  X - computed solutions
     230             : 
     231             :    Level: developer
     232             : 
     233             : .seealso: STMatSolve()
     234             : @*/
     235        3454 : PetscErrorCode STMatMatSolve(ST st,Mat B,Mat X)
     236             : {
     237        3454 :   PetscFunctionBegin;
     238        3454 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     239        3454 :   PetscValidHeaderSpecific(B,MAT_CLASSID,2);
     240        3454 :   PetscValidHeaderSpecific(X,MAT_CLASSID,3);
     241        3454 :   STCheckMatrices(st,1);
     242             : 
     243        3454 :   if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
     244        3454 :   PetscCall(PetscLogEventBegin(ST_MatSolve,st,B,X,0));
     245        3454 :   if (!st->P) PetscCall(MatCopy(B,X,SAME_NONZERO_PATTERN)); /* P=NULL means identity matrix */
     246        3454 :   else PetscCall(KSPMatSolve(st->ksp,B,X));
     247        3454 :   PetscCall(PetscLogEventEnd(ST_MatSolve,st,B,X,0));
     248        3454 :   PetscFunctionReturn(PETSC_SUCCESS);
     249             : }
     250             : 
     251             : /*@
     252             :    STMatSolveTranspose - Solves P^T x = b, where P is the preconditioner matrix of
     253             :    the spectral transformation, using a KSP object stored internally.
     254             : 
     255             :    Collective
     256             : 
     257             :    Input Parameters:
     258             : +  st - the spectral transformation context
     259             : -  b  - right hand side vector
     260             : 
     261             :    Output Parameter:
     262             : .  x - computed solution
     263             : 
     264             :    Level: developer
     265             : 
     266             : .seealso: STMatSolve()
     267             : @*/
     268          18 : PetscErrorCode STMatSolveTranspose(ST st,Vec b,Vec x)
     269             : {
     270          18 :   PetscFunctionBegin;
     271          18 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     272          18 :   PetscValidHeaderSpecific(b,VEC_CLASSID,2);
     273          18 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
     274          18 :   STCheckMatrices(st,1);
     275          18 :   PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
     276          18 :   PetscCall(VecSetErrorIfLocked(x,3));
     277             : 
     278          18 :   if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
     279          18 :   PetscCall(VecLockReadPush(b));
     280          18 :   PetscCall(PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0));
     281          18 :   if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
     282          18 :   else PetscCall(KSPSolveTranspose(st->ksp,b,x));
     283          18 :   PetscCall(PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0));
     284          18 :   PetscCall(VecLockReadPop(b));
     285          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     286             : }
     287             : 
     288             : /*@
     289             :    STMatSolveHermitianTranspose - Solves P^H x = b, where P is the preconditioner matrix of
     290             :    the spectral transformation, using a KSP object stored internally.
     291             : 
     292             :    Collective
     293             : 
     294             :    Input Parameters:
     295             : +  st - the spectral transformation context
     296             : -  b  - right hand side vector
     297             : 
     298             :    Output Parameter:
     299             : .  x - computed solution
     300             : 
     301             :    Level: developer
     302             : 
     303             : .seealso: STMatSolve()
     304             : @*/
     305         593 : PetscErrorCode STMatSolveHermitianTranspose(ST st,Vec b,Vec x)
     306             : {
     307         593 :   Vec  w;
     308             : 
     309         593 :   PetscFunctionBegin;
     310         593 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     311         593 :   PetscValidHeaderSpecific(b,VEC_CLASSID,2);
     312         593 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
     313         593 :   STCheckMatrices(st,1);
     314         593 :   PetscCheck(x!=b,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_IDN,"x and b must be different vectors");
     315         593 :   PetscCall(VecSetErrorIfLocked(x,3));
     316             : 
     317         593 :   if (st->state!=ST_STATE_SETUP) PetscCall(STSetUp(st));
     318         593 :   PetscCall(VecLockReadPush(b));
     319         593 :   PetscCall(PetscLogEventBegin(ST_MatSolveTranspose,st,b,x,0));
     320         593 :   if (!st->P) PetscCall(VecCopy(b,x)); /* P=NULL means identity matrix */
     321             :   else {
     322         593 :     PetscCall(VecDuplicate(b,&w));
     323         593 :     PetscCall(VecCopy(b,w));
     324         593 :     PetscCall(VecConjugate(w));
     325         593 :     PetscCall(KSPSolveTranspose(st->ksp,w,x));
     326         593 :     PetscCall(VecDestroy(&w));
     327         593 :     PetscCall(VecConjugate(x));
     328             :   }
     329         593 :   PetscCall(PetscLogEventEnd(ST_MatSolveTranspose,st,b,x,0));
     330         593 :   PetscCall(VecLockReadPop(b));
     331         593 :   PetscFunctionReturn(PETSC_SUCCESS);
     332             : }
     333             : 
     334        1335 : PetscErrorCode STCheckFactorPackage(ST st)
     335             : {
     336        1335 :   PC             pc;
     337        1335 :   PetscMPIInt    size;
     338        1335 :   PetscBool      flg;
     339        1335 :   MatSolverType  stype;
     340             : 
     341        1335 :   PetscFunctionBegin;
     342        1335 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)st),&size));
     343        1335 :   if (size==1) PetscFunctionReturn(PETSC_SUCCESS);
     344          82 :   PetscCall(KSPGetPC(st->ksp,&pc));
     345          82 :   PetscCall(PCFactorGetMatSolverType(pc,&stype));
     346          82 :   if (stype) {   /* currently selected PC is a factorization */
     347           0 :     PetscCall(PetscStrcmp(stype,MATSOLVERPETSC,&flg));
     348           0 :     PetscCheck(!flg,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");
     349             :   }
     350          82 :   PetscFunctionReturn(PETSC_SUCCESS);
     351             : }
     352             : 
     353             : /*@
     354             :    STSetKSP - Sets the KSP object associated with the spectral
     355             :    transformation.
     356             : 
     357             :    Collective
     358             : 
     359             :    Input Parameters:
     360             : +  st   - the spectral transformation context
     361             : -  ksp  - the linear system context
     362             : 
     363             :    Level: advanced
     364             : 
     365             : .seealso: STGetKSP()
     366             : @*/
     367           6 : PetscErrorCode STSetKSP(ST st,KSP ksp)
     368             : {
     369           6 :   PetscFunctionBegin;
     370           6 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     371           6 :   PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
     372           6 :   PetscCheckSameComm(st,1,ksp,2);
     373           6 :   STCheckNotSeized(st,1);
     374           6 :   PetscCall(PetscObjectReference((PetscObject)ksp));
     375           6 :   PetscCall(KSPDestroy(&st->ksp));
     376           6 :   st->ksp = ksp;
     377           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     378             : }
     379             : 
     380             : /*@
     381             :    STGetKSP - Gets the KSP object associated with the spectral
     382             :    transformation.
     383             : 
     384             :    Collective
     385             : 
     386             :    Input Parameter:
     387             : .  st - the spectral transformation context
     388             : 
     389             :    Output Parameter:
     390             : .  ksp  - the linear system context
     391             : 
     392             :    Level: intermediate
     393             : 
     394             : .seealso: STSetKSP()
     395             : @*/
     396        2638 : PetscErrorCode STGetKSP(ST st,KSP* ksp)
     397             : {
     398        2638 :   PetscFunctionBegin;
     399        2638 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     400        2638 :   PetscAssertPointer(ksp,2);
     401        2638 :   if (!st->ksp) {
     402         778 :     PetscCall(KSPCreate(PetscObjectComm((PetscObject)st),&st->ksp));
     403         778 :     PetscCall(PetscObjectIncrementTabLevel((PetscObject)st->ksp,(PetscObject)st,1));
     404         778 :     PetscCall(KSPSetOptionsPrefix(st->ksp,((PetscObject)st)->prefix));
     405         778 :     PetscCall(KSPAppendOptionsPrefix(st->ksp,"st_"));
     406         778 :     PetscCall(PetscObjectSetOptions((PetscObject)st->ksp,((PetscObject)st)->options));
     407         778 :     PetscCall(KSPSetTolerances(st->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
     408             :   }
     409        2638 :   *ksp = st->ksp;
     410        2638 :   PetscFunctionReturn(PETSC_SUCCESS);
     411             : }
     412             : 
     413           1 : PetscErrorCode STCheckNullSpace_Default(ST st,BV V)
     414             : {
     415           1 :   PetscInt       nc,i,c;
     416           1 :   PetscReal      norm;
     417           1 :   Vec            *T,w,vi;
     418           1 :   Mat            A;
     419           1 :   PC             pc;
     420           1 :   MatNullSpace   nullsp;
     421             : 
     422           1 :   PetscFunctionBegin;
     423           1 :   PetscCall(BVGetNumConstraints(V,&nc));
     424           1 :   PetscCall(PetscMalloc1(nc,&T));
     425           1 :   if (!st->ksp) PetscCall(STGetKSP(st,&st->ksp));
     426           1 :   PetscCall(KSPGetPC(st->ksp,&pc));
     427           1 :   PetscCall(PCGetOperators(pc,&A,NULL));
     428           1 :   PetscCall(MatCreateVecs(A,NULL,&w));
     429             :   c = 0;
     430           2 :   for (i=0;i<nc;i++) {
     431           1 :     PetscCall(BVGetColumn(V,-nc+i,&vi));
     432           1 :     PetscCall(MatMult(A,vi,w));
     433           1 :     PetscCall(VecNorm(w,NORM_2,&norm));
     434           1 :     if (norm < 10.0*PETSC_SQRT_MACHINE_EPSILON) {
     435           1 :       PetscCall(PetscInfo(st,"Vector %" PetscInt_FMT " included in the nullspace of OP, norm=%g\n",i,(double)norm));
     436           1 :       PetscCall(BVCreateVec(V,T+c));
     437           1 :       PetscCall(VecCopy(vi,T[c]));
     438           1 :       PetscCall(VecNormalize(T[c],NULL));
     439           1 :       c++;
     440           0 :     } else PetscCall(PetscInfo(st,"Vector %" PetscInt_FMT " discarded as possible nullspace of OP, norm=%g\n",i,(double)norm));
     441           1 :     PetscCall(BVRestoreColumn(V,-nc+i,&vi));
     442             :   }
     443           1 :   PetscCall(VecDestroy(&w));
     444           1 :   if (c>0) {
     445           1 :     PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)st),PETSC_FALSE,c,T,&nullsp));
     446           1 :     PetscCall(MatSetNullSpace(A,nullsp));
     447           1 :     PetscCall(MatNullSpaceDestroy(&nullsp));
     448           1 :     PetscCall(VecDestroyVecs(c,&T));
     449           0 :   } else PetscCall(PetscFree(T));
     450           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     451             : }
     452             : 
     453             : /*@
     454             :    STCheckNullSpace - Tests if constraint vectors are nullspace vectors.
     455             : 
     456             :    Collective
     457             : 
     458             :    Input Parameters:
     459             : +  st - the spectral transformation context
     460             : -  V  - basis vectors to be checked
     461             : 
     462             :    Notes:
     463             :    Given a basis vectors object, this function tests each of its constraint
     464             :    vectors to be a nullspace vector of the coefficient matrix of the
     465             :    associated KSP object. All these nullspace vectors are passed to the KSP
     466             :    object.
     467             : 
     468             :    This function allows handling singular pencils and solving some problems
     469             :    in which the nullspace is important (see the users guide for details).
     470             : 
     471             :    Level: developer
     472             : 
     473             : .seealso: EPSSetDeflationSpace()
     474             : @*/
     475          14 : PetscErrorCode STCheckNullSpace(ST st,BV V)
     476             : {
     477          14 :   PetscInt       nc;
     478             : 
     479          14 :   PetscFunctionBegin;
     480          14 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     481          14 :   PetscValidHeaderSpecific(V,BV_CLASSID,2);
     482          14 :   PetscValidType(st,1);
     483          14 :   PetscCheckSameComm(st,1,V,2);
     484          14 :   PetscCheck(st->state,PetscObjectComm((PetscObject)st),PETSC_ERR_ARG_WRONGSTATE,"Must call STSetUp() first");
     485             : 
     486          14 :   PetscCall(BVGetNumConstraints(V,&nc));
     487          14 :   if (nc) PetscTryTypeMethod(st,checknullspace,V);
     488          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     489             : }

Generated by: LCOV version 1.14