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             :    Shift spectral transformation, applies (A + sigma I) as operator, or
      12             :    inv(B)(A + sigma B) for generalized problems
      13             : */
      14             : 
      15             : #include <slepc/private/stimpl.h>
      16             : 
      17             : /*
      18             :    Special STApply() for the BSE structured matrix
      19             : 
      20             :        H = [ R  C; -C^H -R^T ].
      21             : 
      22             :    Assumes that H is a MATNEST and x,y are VECNEST.
      23             : 
      24             :    Only the upper part of the product y=H*x is computed.
      25             : 
      26             :        y1 = R*x1+C*x2
      27             : 
      28             :    The bottom part of the input vector x2 should normally contain
      29             :    either conj(x1) or -conj(x1).
      30             :    The bottom part of the output vector y2 is not referenced.
      31             : */
      32        1083 : static PetscErrorCode STApply_Shift_BSE(ST st,Vec x,Vec y)
      33             : {
      34        1083 :   Mat   H,R,C;
      35        1083 :   Vec   x1,x2,y1;
      36             : 
      37        1083 :   PetscFunctionBegin;
      38        1083 :   H = st->T[0];
      39        1083 :   PetscCall(MatNestGetSubMat(H,0,0,&R));
      40        1083 :   PetscCall(MatNestGetSubMat(H,0,1,&C));
      41        1083 :   PetscCall(VecNestGetSubVec(x,0,&x1));
      42        1083 :   PetscCall(VecNestGetSubVec(x,1,&x2));
      43        1083 :   PetscCall(VecNestGetSubVec(y,0,&y1));
      44        1083 :   PetscCall(MatMult(C,x2,y1));
      45        1083 :   PetscCall(MatMultAdd(R,x1,y1,y1));
      46        1083 :   PetscFunctionReturn(PETSC_SUCCESS);
      47             : }
      48             : 
      49      831891 : static PetscErrorCode STBackTransform_Shift(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
      50             : {
      51      831891 :   PetscInt j;
      52             : 
      53      831891 :   PetscFunctionBegin;
      54     2486090 :   for (j=0;j<n;j++) {
      55     1654199 :     eigr[j] += st->sigma;
      56             :   }
      57      831891 :   PetscFunctionReturn(PETSC_SUCCESS);
      58             : }
      59             : 
      60         561 : static PetscErrorCode STPostSolve_Shift(ST st)
      61             : {
      62         561 :   PetscFunctionBegin;
      63         561 :   if (st->matmode == ST_MATMODE_INPLACE) {
      64           4 :     if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
      65           3 :     else PetscCall(MatShift(st->A[0],st->sigma));
      66           4 :     st->Astate[0] = ((PetscObject)st->A[0])->state;
      67           4 :     st->state   = ST_STATE_INITIAL;
      68           4 :     st->opready = PETSC_FALSE;
      69             :   }
      70         561 :   PetscFunctionReturn(PETSC_SUCCESS);
      71             : }
      72             : 
      73             : /*
      74             :    Operator (shift):
      75             :                Op               P         M
      76             :    if nmat=1:  A-sI             NULL      A-sI
      77             :    if nmat=2:  B^-1 (A-sB)      B         A-sB
      78             : */
      79         412 : static PetscErrorCode STComputeOperator_Shift(ST st)
      80             : {
      81         412 :   PetscFunctionBegin;
      82         412 :   st->usesksp = (st->nmat>1)? PETSC_TRUE: PETSC_FALSE;
      83         412 :   PetscCall(PetscObjectReference((PetscObject)st->A[1]));
      84         412 :   PetscCall(MatDestroy(&st->T[1]));
      85         412 :   st->T[1] = st->A[1];
      86         412 :   PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[0]));
      87         412 :   if (st->nmat>1) PetscCall(PetscObjectReference((PetscObject)st->T[1]));
      88         412 :   PetscCall(MatDestroy(&st->P));
      89         412 :   st->P = (st->nmat>1)? st->T[1]: NULL;
      90         412 :   st->M = st->T[0];
      91         412 :   if (st->nmat>1 && st->Psplit) {  /* build custom preconditioner from the split matrices */
      92           1 :     PetscCall(MatDestroy(&st->Pmat));
      93           1 :     PetscCall(PetscObjectReference((PetscObject)st->Psplit[1]));
      94           1 :     st->Pmat = st->Psplit[1];
      95             :   }
      96         412 :   PetscFunctionReturn(PETSC_SUCCESS);
      97             : }
      98             : 
      99         489 : static PetscErrorCode STSetUp_Shift(ST st)
     100             : {
     101         489 :   PetscInt       k,nc,nmat=st->nmat;
     102         489 :   PetscScalar    *coeffs=NULL;
     103             : 
     104         489 :   PetscFunctionBegin;
     105         489 :   if (nmat>1) PetscCall(STSetWorkVecs(st,1));
     106         489 :   if (nmat>2) {  /* set-up matrices for polynomial eigenproblems */
     107          77 :     if (st->transform) {
     108          14 :       nc = (nmat*(nmat+1))/2;
     109          14 :       PetscCall(PetscMalloc1(nc,&coeffs));
     110             :       /* Compute coeffs */
     111          14 :       PetscCall(STCoeffs_Monomial(st,coeffs));
     112             :       /* T[n] = A_n */
     113          14 :       k = nmat-1;
     114          14 :       PetscCall(PetscObjectReference((PetscObject)st->A[k]));
     115          14 :       PetscCall(MatDestroy(&st->T[k]));
     116          14 :       st->T[k] = st->A[k];
     117          45 :       for (k=0;k<nmat-1;k++) PetscCall(STMatMAXPY_Private(st,nmat>2?st->sigma:-st->sigma,0.0,k,coeffs?coeffs+((nmat-k)*(nmat-k-1))/2:NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[k]));
     118          14 :       PetscCall(PetscFree(coeffs));
     119          14 :       PetscCall(PetscObjectReference((PetscObject)st->T[nmat-1]));
     120          14 :       PetscCall(MatDestroy(&st->P));
     121          14 :       st->P = st->T[nmat-1];
     122          14 :       if (st->Psplit) {  /* build custom preconditioner from the split matrices */
     123           1 :         PetscCall(STMatMAXPY_Private(st,st->sigma,0.0,nmat-1,coeffs?coeffs:NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
     124             :       }
     125          14 :       PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
     126             :     } else {
     127         276 :       for (k=0;k<nmat;k++) {
     128         213 :         PetscCall(PetscObjectReference((PetscObject)st->A[k]));
     129         213 :         PetscCall(MatDestroy(&st->T[k]));
     130         213 :         st->T[k] = st->A[k];
     131             :       }
     132             :     }
     133             :   }
     134         489 :   if (st->P) PetscCall(KSPSetUp(st->ksp));
     135         489 :   if (st->structured) st->ops->apply = STApply_Shift_BSE;
     136         489 :   PetscFunctionReturn(PETSC_SUCCESS);
     137             : }
     138             : 
     139          18 : static PetscErrorCode STSetShift_Shift(ST st,PetscScalar newshift)
     140             : {
     141          18 :   PetscInt       k,nc,nmat=PetscMax(st->nmat,2);
     142          18 :   PetscScalar    *coeffs=NULL;
     143             : 
     144          18 :   PetscFunctionBegin;
     145          18 :   if (st->transform) {
     146          16 :     if (st->matmode == ST_MATMODE_COPY && nmat>2) {
     147           2 :       nc = (nmat*(nmat+1))/2;
     148           2 :       PetscCall(PetscMalloc1(nc,&coeffs));
     149             :       /* Compute coeffs */
     150           2 :       PetscCall(STCoeffs_Monomial(st,coeffs));
     151             :     }
     152          38 :     for (k=0;k<nmat-1;k++) PetscCall(STMatMAXPY_Private(st,nmat>2?newshift:-newshift,nmat>2?st->sigma:-st->sigma,k,coeffs?coeffs+((nmat-k)*(nmat-k-1))/2:NULL,PETSC_FALSE,PETSC_FALSE,&st->T[k]));
     153          16 :     if (st->matmode == ST_MATMODE_COPY && nmat>2) PetscCall(PetscFree(coeffs));
     154          16 :     if (st->nmat<=2) st->M = st->T[0];
     155             :   }
     156          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     157             : }
     158             : 
     159         561 : SLEPC_EXTERN PetscErrorCode STCreate_Shift(ST st)
     160             : {
     161         561 :   PetscFunctionBegin;
     162         561 :   st->usesksp = PETSC_TRUE;
     163             : 
     164         561 :   st->ops->apply           = STApply_Generic;
     165         561 :   st->ops->applytrans      = STApplyTranspose_Generic;
     166         561 :   st->ops->applyhermtrans  = STApplyHermitianTranspose_Generic;
     167         561 :   st->ops->backtransform   = STBackTransform_Shift;
     168         561 :   st->ops->setshift        = STSetShift_Shift;
     169         561 :   st->ops->getbilinearform = STGetBilinearForm_Default;
     170         561 :   st->ops->setup           = STSetUp_Shift;
     171         561 :   st->ops->computeoperator = STComputeOperator_Shift;
     172         561 :   st->ops->postsolve       = STPostSolve_Shift;
     173         561 :   st->ops->setdefaultksp   = STSetDefaultKSP_Default;
     174         561 :   PetscFunctionReturn(PETSC_SUCCESS);
     175             : }