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      201856 : static PetscErrorCode STBackTransform_Sinvert(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
      17             : {
      18      201856 :   PetscInt    j;
      19             : #if !defined(PETSC_USE_COMPLEX)
      20             :   PetscScalar t;
      21             : #endif
      22             : 
      23      201856 :   PetscFunctionBegin;
      24             : #if !defined(PETSC_USE_COMPLEX)
      25             :   for (j=0;j<n;j++) {
      26             :     if (eigi[j] == 0) eigr[j] = 1.0 / eigr[j] + st->sigma;
      27             :     else {
      28             :       t = eigr[j] * eigr[j] + eigi[j] * eigi[j];
      29             :       eigr[j] = eigr[j] / t + st->sigma;
      30             :       eigi[j] = - eigi[j] / t;
      31             :     }
      32             :   }
      33             : #else
      34      616164 :   for (j=0;j<n;j++) {
      35      414308 :     eigr[j] = 1.0 / eigr[j] + st->sigma;
      36             :   }
      37             : #endif
      38      201856 :   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           3 :     if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
      46           2 :     else PetscCall(MatShift(st->A[0],st->sigma));
      47           3 :     st->Astate[0] = ((PetscObject)st->A[0])->state;
      48           3 :     st->state   = ST_STATE_INITIAL;
      49           3 :     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         140 : static PetscErrorCode STComputeOperator_Sinvert(ST st)
      61             : {
      62         140 :   PetscFunctionBegin;
      63             :   /* if the user did not set the shift, use the target value */
      64         140 :   if (!st->sigma_set) st->sigma = st->defsigma;
      65         140 :   PetscCall(PetscObjectReference((PetscObject)st->A[1]));
      66         140 :   PetscCall(MatDestroy(&st->T[0]));
      67         140 :   st->T[0] = st->A[1];
      68         140 :   PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[1]));
      69         140 :   PetscCall(PetscObjectReference((PetscObject)st->T[1]));
      70         140 :   PetscCall(MatDestroy(&st->P));
      71         140 :   st->P = st->T[1];
      72         140 :   st->M = (st->nmat>1)? st->T[0]: NULL;
      73         140 :   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         140 :   PetscFunctionReturn(PETSC_SUCCESS);
      77             : }
      78             : 
      79         255 : static PetscErrorCode STSetUp_Sinvert(ST st)
      80             : {
      81         255 :   PetscInt       k,nc,nmat=st->nmat,nsub;
      82         255 :   PetscScalar    *coeffs=NULL;
      83         255 :   PetscBool      islu;
      84         255 :   KSP            *subksp;
      85         255 :   PC             pc,pcsub;
      86         255 :   Mat            A00;
      87         255 :   MatType        type;
      88         255 :   PetscBool      flg;
      89         255 :   char           str[64];
      90             : 
      91         255 :   PetscFunctionBegin;
      92         255 :   if (nmat>1) PetscCall(STSetWorkVecs(st,1));
      93             :   /* if the user did not set the shift, use the target value */
      94         255 :   if (!st->sigma_set) st->sigma = st->defsigma;
      95         255 :   if (nmat>2) {  /* set-up matrices for polynomial eigenproblems */
      96         111 :     if (st->transform) {
      97          23 :       nc = (nmat*(nmat+1))/2;
      98          23 :       PetscCall(PetscMalloc1(nc,&coeffs));
      99             :       /* Compute coeffs */
     100          23 :       PetscCall(STCoeffs_Monomial(st,coeffs));
     101             :       /* T[0] = A_n */
     102          23 :       k = nmat-1;
     103          23 :       PetscCall(PetscObjectReference((PetscObject)st->A[k]));
     104          23 :       PetscCall(MatDestroy(&st->T[0]));
     105          23 :       st->T[0] = st->A[k];
     106          72 :       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]));
     107          23 :       PetscCall(PetscFree(coeffs));
     108          23 :       PetscCall(PetscObjectReference((PetscObject)st->T[nmat-1]));
     109          23 :       PetscCall(MatDestroy(&st->P));
     110          23 :       st->P = st->T[nmat-1];
     111          23 :       if (st->Psplit) {  /* build custom preconditioner from the split matrices */
     112           2 :         PetscCall(STMatMAXPY_Private(st,st->sigma,0.0,0,coeffs?coeffs+((nmat-1)*nmat)/2:NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
     113             :       }
     114          23 :       PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
     115             :     } else {
     116         485 :       for (k=0;k<nmat;k++) {
     117         397 :         PetscCall(PetscObjectReference((PetscObject)st->A[k]));
     118         397 :         PetscCall(MatDestroy(&st->T[k]));
     119         397 :         st->T[k] = st->A[k];
     120             :       }
     121             :     }
     122             :   }
     123         255 :   if (st->structured) {
     124             :     /* ./ex55 -st_type sinvert -eps_target 0 -eps_target_magnitude
     125             :               -st_ksp_type preonly -st_pc_type fieldsplit
     126             :               -st_fieldsplit_0_ksp_type preonly -st_fieldsplit_0_pc_type lu
     127             :               -st_fieldsplit_1_ksp_type preonly -st_fieldsplit_1_pc_type lu
     128             :               -st_pc_fieldsplit_type schur -st_pc_fieldsplit_schur_fact_type full
     129             :               -st_pc_fieldsplit_schur_precondition full */
     130           3 :     PetscCall(KSPGetPC(st->ksp,&pc));
     131           3 :     PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCLU,&islu));
     132           3 :     if (islu) {  /* assume PCLU means the user has not set any options */
     133           3 :       PetscCall(KSPSetType(st->ksp,KSPPREONLY));
     134           3 :       PetscCall(PCSetType(pc,PCFIELDSPLIT));
     135           3 :       PetscCall(PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR));
     136           3 :       PetscCall(PCFieldSplitSetSchurFactType(pc,PC_FIELDSPLIT_SCHUR_FACT_FULL));
     137           3 :       PetscCall(PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_FULL,NULL));
     138             :       /* hack to set Mat type of Schur complement equal to A00, assumes default prefixes */
     139           3 :       PetscCall(PetscOptionsHasName(((PetscObject)st)->options,((PetscObject)st)->prefix,"-st_fieldsplit_1_explicit_operator_mat_type",&flg));
     140           3 :       if (!flg) {
     141           3 :         PetscCall(MatNestGetSubMat(st->A[0],0,0,&A00));
     142           3 :         PetscCall(MatGetType(A00,&type));
     143           3 :         PetscCall(PetscSNPrintf(str,sizeof(str),"-%sst_fieldsplit_1_explicit_operator_mat_type %s",((PetscObject)st)->prefix?((PetscObject)st)->prefix:"",type));
     144           3 :         PetscCall(PetscOptionsInsertString(((PetscObject)st)->options,str));
     145             :       }
     146             :       /* set preonly+lu on block solvers */
     147           3 :       PetscCall(KSPSetUp(st->ksp));
     148           3 :       PetscCall(PCFieldSplitGetSubKSP(pc,&nsub,&subksp));
     149           3 :       PetscCall(KSPSetType(subksp[0],KSPPREONLY));
     150           3 :       PetscCall(KSPGetPC(subksp[0],&pcsub));
     151           3 :       PetscCall(PCSetType(pcsub,PCLU));
     152           3 :       PetscCall(KSPSetType(subksp[1],KSPPREONLY));
     153           3 :       PetscCall(KSPGetPC(subksp[1],&pcsub));
     154           3 :       PetscCall(PCSetType(pcsub,PCLU));
     155           3 :       PetscCall(PetscFree(subksp));
     156             :     }
     157             :   } else {
     158         252 :     if (st->P) PetscCall(KSPSetUp(st->ksp));
     159             :   }
     160         255 :   PetscFunctionReturn(PETSC_SUCCESS);
     161             : }
     162             : 
     163         623 : static PetscErrorCode STSetShift_Sinvert(ST st,PetscScalar newshift)
     164             : {
     165         623 :   PetscInt       nmat=PetscMax(st->nmat,2),k,nc;
     166         623 :   PetscScalar    *coeffs=NULL;
     167             : 
     168         623 :   PetscFunctionBegin;
     169         623 :   if (st->transform) {
     170         590 :     if (st->matmode == ST_MATMODE_COPY && nmat>2) {
     171           2 :       nc = (nmat*(nmat+1))/2;
     172           2 :       PetscCall(PetscMalloc1(nc,&coeffs));
     173             :       /* Compute coeffs */
     174           2 :       PetscCall(STCoeffs_Monomial(st,coeffs));
     175             :     }
     176        1186 :     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]));
     177         590 :     if (st->matmode == ST_MATMODE_COPY && nmat>2) PetscCall(PetscFree(coeffs));
     178         590 :     if (st->P!=st->T[nmat-1]) {
     179          18 :       PetscCall(PetscObjectReference((PetscObject)st->T[nmat-1]));
     180          18 :       PetscCall(MatDestroy(&st->P));
     181          18 :       st->P = st->T[nmat-1];
     182             :     }
     183         590 :     if (st->Psplit) {  /* build custom preconditioner from the split matrices */
     184          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));
     185             :     }
     186             :   }
     187         623 :   if (st->P) {
     188         621 :     PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
     189         621 :     PetscCall(KSPSetUp(st->ksp));
     190             :   }
     191         623 :   PetscFunctionReturn(PETSC_SUCCESS);
     192             : }
     193             : 
     194         240 : SLEPC_EXTERN PetscErrorCode STCreate_Sinvert(ST st)
     195             : {
     196         240 :   PetscFunctionBegin;
     197         240 :   st->usesksp = PETSC_TRUE;
     198             : 
     199         240 :   st->ops->apply           = STApply_Generic;
     200         240 :   st->ops->applytrans      = STApplyTranspose_Generic;
     201         240 :   st->ops->applyhermtrans  = STApplyHermitianTranspose_Generic;
     202         240 :   st->ops->backtransform   = STBackTransform_Sinvert;
     203         240 :   st->ops->setshift        = STSetShift_Sinvert;
     204         240 :   st->ops->getbilinearform = STGetBilinearForm_Default;
     205         240 :   st->ops->setup           = STSetUp_Sinvert;
     206         240 :   st->ops->computeoperator = STComputeOperator_Sinvert;
     207         240 :   st->ops->postsolve       = STPostSolve_Sinvert;
     208         240 :   st->ops->checknullspace  = STCheckNullSpace_Default;
     209         240 :   st->ops->setdefaultksp   = STSetDefaultKSP_Default;
     210         240 :   PetscFunctionReturn(PETSC_SUCCESS);
     211             : }