LCOV - code coverage report
Current view: top level - sys/classes/st/impls/sinvert - sinvert.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 98 98 100.0 %
Date: 2024-03-29 00:55:19 Functions: 6 6 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             :    Implements the shift-and-invert technique for eigenvalue problems
      12             : */
      13             : 
      14             : #include <slepc/private/stimpl.h>
      15             : 
      16      131172 : static PetscErrorCode STBackTransform_Sinvert(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
      17             : {
      18      131172 :   PetscInt    j;
      19             : #if !defined(PETSC_USE_COMPLEX)
      20      131172 :   PetscScalar t;
      21             : #endif
      22             : 
      23      131172 :   PetscFunctionBegin;
      24             : #if !defined(PETSC_USE_COMPLEX)
      25      401951 :   for (j=0;j<n;j++) {
      26      270779 :     if (eigi[j] == 0) eigr[j] = 1.0 / eigr[j] + st->sigma;
      27             :     else {
      28       51933 :       t = eigr[j] * eigr[j] + eigi[j] * eigi[j];
      29       51933 :       eigr[j] = eigr[j] / t + st->sigma;
      30       51933 :       eigi[j] = - eigi[j] / t;
      31             :     }
      32             :   }
      33             : #else
      34             :   for (j=0;j<n;j++) {
      35             :     eigr[j] = 1.0 / eigr[j] + st->sigma;
      36             :   }
      37             : #endif
      38      131172 :   PetscFunctionReturn(PETSC_SUCCESS);
      39             : }
      40             : 
      41         202 : static PetscErrorCode STPostSolve_Sinvert(ST st)
      42             : {
      43         202 :   PetscFunctionBegin;
      44         202 :   if (st->matmode == ST_MATMODE_INPLACE) {
      45           2 :     if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
      46           1 :     else PetscCall(MatShift(st->A[0],st->sigma));
      47           2 :     st->Astate[0] = ((PetscObject)st->A[0])->state;
      48           2 :     st->state   = ST_STATE_INITIAL;
      49           2 :     st->opready = PETSC_FALSE;
      50             :   }
      51         202 :   PetscFunctionReturn(PETSC_SUCCESS);
      52             : }
      53             : 
      54             : /*
      55             :    Operator (sinvert):
      56             :                Op               P         M
      57             :    if nmat=1:  (A-sI)^-1        A-sI      NULL
      58             :    if nmat=2:  (A-sB)^-1 B      A-sB      B
      59             : */
      60         142 : static PetscErrorCode STComputeOperator_Sinvert(ST st)
      61             : {
      62         142 :   PetscFunctionBegin;
      63             :   /* if the user did not set the shift, use the target value */
      64         142 :   if (!st->sigma_set) st->sigma = st->defsigma;
      65         142 :   PetscCall(PetscObjectReference((PetscObject)st->A[1]));
      66         142 :   PetscCall(MatDestroy(&st->T[0]));
      67         142 :   st->T[0] = st->A[1];
      68         142 :   PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[1]));
      69         142 :   PetscCall(PetscObjectReference((PetscObject)st->T[1]));
      70         142 :   PetscCall(MatDestroy(&st->P));
      71         142 :   st->P = st->T[1];
      72         142 :   st->M = (st->nmat>1)? st->T[0]: NULL;
      73         142 :   if (st->Psplit) {  /* build custom preconditioner from the split matrices */
      74           7 :     PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
      75             :   }
      76         142 :   PetscFunctionReturn(PETSC_SUCCESS);
      77             : }
      78             : 
      79         251 : static PetscErrorCode STSetUp_Sinvert(ST st)
      80             : {
      81         251 :   PetscInt       k,nc,nmat=st->nmat;
      82         251 :   PetscScalar    *coeffs=NULL;
      83             : 
      84         251 :   PetscFunctionBegin;
      85         251 :   if (nmat>1) PetscCall(STSetWorkVecs(st,1));
      86             :   /* if the user did not set the shift, use the target value */
      87         251 :   if (!st->sigma_set) st->sigma = st->defsigma;
      88         251 :   if (nmat>2) {  /* set-up matrices for polynomial eigenproblems */
      89         105 :     if (st->transform) {
      90          22 :       nc = (nmat*(nmat+1))/2;
      91          22 :       PetscCall(PetscMalloc1(nc,&coeffs));
      92             :       /* Compute coeffs */
      93          22 :       PetscCall(STCoeffs_Monomial(st,coeffs));
      94             :       /* T[0] = A_n */
      95          22 :       k = nmat-1;
      96          22 :       PetscCall(PetscObjectReference((PetscObject)st->A[k]));
      97          22 :       PetscCall(MatDestroy(&st->T[0]));
      98          22 :       st->T[0] = st->A[k];
      99          69 :       for (k=1;k<nmat;k++) PetscCall(STMatMAXPY_Private(st,nmat>2?st->sigma:-st->sigma,0.0,nmat-k-1,coeffs?coeffs+(k*(k+1))/2:NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[k]));
     100          22 :       PetscCall(PetscFree(coeffs));
     101          22 :       PetscCall(PetscObjectReference((PetscObject)st->T[nmat-1]));
     102          22 :       PetscCall(MatDestroy(&st->P));
     103          22 :       st->P = st->T[nmat-1];
     104          22 :       if (st->Psplit) {  /* build custom preconditioner from the split matrices */
     105           2 :         PetscCall(STMatMAXPY_Private(st,st->sigma,0.0,0,coeffs?coeffs+((nmat-1)*nmat)/2:NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
     106             :       }
     107          22 :       PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
     108             :     } else {
     109         465 :       for (k=0;k<nmat;k++) {
     110         382 :         PetscCall(PetscObjectReference((PetscObject)st->A[k]));
     111         382 :         PetscCall(MatDestroy(&st->T[k]));
     112         382 :         st->T[k] = st->A[k];
     113             :       }
     114             :     }
     115             :   }
     116         251 :   if (st->P) PetscCall(KSPSetUp(st->ksp));
     117         251 :   PetscFunctionReturn(PETSC_SUCCESS);
     118             : }
     119             : 
     120         409 : static PetscErrorCode STSetShift_Sinvert(ST st,PetscScalar newshift)
     121             : {
     122         409 :   PetscInt       nmat=PetscMax(st->nmat,2),k,nc;
     123         409 :   PetscScalar    *coeffs=NULL;
     124             : 
     125         409 :   PetscFunctionBegin;
     126         409 :   if (st->transform) {
     127         377 :     if (st->matmode == ST_MATMODE_COPY && nmat>2) {
     128           2 :       nc = (nmat*(nmat+1))/2;
     129           2 :       PetscCall(PetscMalloc1(nc,&coeffs));
     130             :       /* Compute coeffs */
     131           2 :       PetscCall(STCoeffs_Monomial(st,coeffs));
     132             :     }
     133         760 :     for (k=1;k<nmat;k++) PetscCall(STMatMAXPY_Private(st,nmat>2?newshift:-newshift,nmat>2?st->sigma:-st->sigma,nmat-k-1,coeffs?coeffs+(k*(k+1))/2:NULL,PETSC_FALSE,PETSC_FALSE,&st->T[k]));
     134         377 :     if (st->matmode == ST_MATMODE_COPY && nmat>2) PetscCall(PetscFree(coeffs));
     135         377 :     if (st->P!=st->T[nmat-1]) {
     136          12 :       PetscCall(PetscObjectReference((PetscObject)st->T[nmat-1]));
     137          12 :       PetscCall(MatDestroy(&st->P));
     138          12 :       st->P = st->T[nmat-1];
     139             :     }
     140         377 :     if (st->Psplit) {  /* build custom preconditioner from the split matrices */
     141          35 :       PetscCall(STMatMAXPY_Private(st,nmat>2?newshift:-newshift,nmat>2?st->sigma:-st->sigma,0,coeffs?coeffs+((nmat-1)*nmat)/2:NULL,PETSC_FALSE,PETSC_TRUE,&st->Pmat));
     142             :     }
     143             :   }
     144         409 :   if (st->P) {
     145         407 :     PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
     146         407 :     PetscCall(KSPSetUp(st->ksp));
     147             :   }
     148         409 :   PetscFunctionReturn(PETSC_SUCCESS);
     149             : }
     150             : 
     151         236 : SLEPC_EXTERN PetscErrorCode STCreate_Sinvert(ST st)
     152             : {
     153         236 :   PetscFunctionBegin;
     154         236 :   st->usesksp = PETSC_TRUE;
     155             : 
     156         236 :   st->ops->apply           = STApply_Generic;
     157         236 :   st->ops->applytrans      = STApplyTranspose_Generic;
     158         236 :   st->ops->applyhermtrans  = STApplyHermitianTranspose_Generic;
     159         236 :   st->ops->backtransform   = STBackTransform_Sinvert;
     160         236 :   st->ops->setshift        = STSetShift_Sinvert;
     161         236 :   st->ops->getbilinearform = STGetBilinearForm_Default;
     162         236 :   st->ops->setup           = STSetUp_Sinvert;
     163         236 :   st->ops->computeoperator = STComputeOperator_Sinvert;
     164         236 :   st->ops->postsolve       = STPostSolve_Sinvert;
     165         236 :   st->ops->checknullspace  = STCheckNullSpace_Default;
     166         236 :   st->ops->setdefaultksp   = STSetDefaultKSP_Default;
     167         236 :   PetscFunctionReturn(PETSC_SUCCESS);
     168             : }

Generated by: LCOV version 1.14