LCOV - code coverage report
Current view: top level - sys/classes/st/impls/cayley - cayley.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 155 170 91.2 %
Date: 2017-12-10 00:15:52 Functions: 16 17 94.1 %

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2017, 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 Cayley spectral transform
      12             : */
      13             : 
      14             : #include <slepc/private/stimpl.h>          /*I "slepcst.h" I*/
      15             : 
      16             : typedef struct {
      17             :   PetscScalar nu;
      18             :   PetscBool   nu_set;
      19             : } ST_CAYLEY;
      20             : 
      21         110 : PetscErrorCode STApply_Cayley(ST st,Vec x,Vec y)
      22             : {
      23             :   PetscErrorCode ierr;
      24             : 
      25         110 :   PetscFunctionBegin;
      26             :   /* standard eigenproblem: y = (A - sI)^-1 (A + tI)x */
      27             :   /* generalized eigenproblem: y = (A - sB)^-1 (A + tB)x */
      28         110 :   ierr = MatMult(st->T[0],x,st->work[0]);CHKERRQ(ierr);
      29         110 :   ierr = STMatSolve(st,st->work[0],y);CHKERRQ(ierr);
      30         110 :   PetscFunctionReturn(0);
      31             : }
      32             : 
      33           6 : PetscErrorCode STApplyTranspose_Cayley(ST st,Vec x,Vec y)
      34             : {
      35             :   PetscErrorCode ierr;
      36             : 
      37           6 :   PetscFunctionBegin;
      38             :   /* standard eigenproblem: y =  (A + tI)^T (A - sI)^-T x */
      39             :   /* generalized eigenproblem: y = (A + tB)^T (A - sB)^-T x */
      40           6 :   ierr = STMatSolveTranspose(st,x,st->work[0]);CHKERRQ(ierr);
      41           6 :   ierr = MatMultTranspose(st->T[0],st->work[0],y);CHKERRQ(ierr);
      42           6 :   PetscFunctionReturn(0);
      43             : }
      44             : 
      45          10 : static PetscErrorCode MatMult_Cayley(Mat B,Vec x,Vec y)
      46             : {
      47             :   PetscErrorCode ierr;
      48             :   ST             st;
      49             :   ST_CAYLEY      *ctx;
      50             :   PetscScalar    nu;
      51             : 
      52          10 :   PetscFunctionBegin;
      53          10 :   ierr = MatShellGetContext(B,(void**)&st);CHKERRQ(ierr);
      54          10 :   ctx = (ST_CAYLEY*)st->data;
      55          10 :   nu = ctx->nu;
      56             : 
      57          10 :   if (st->shift_matrix == ST_MATMODE_INPLACE) { nu = nu + st->sigma; };
      58             : 
      59          10 :   if (st->nmat>1) {
      60             :     /* generalized eigenproblem: y = (A + tB)x */
      61           4 :     ierr = MatMult(st->A[0],x,y);CHKERRQ(ierr);
      62           4 :     ierr = MatMult(st->A[1],x,st->work[1]);CHKERRQ(ierr);
      63           4 :     ierr = VecAXPY(y,nu,st->work[1]);CHKERRQ(ierr);
      64             :   } else {
      65             :     /* standard eigenproblem: y = (A + tI)x */
      66           6 :     ierr = MatMult(st->A[0],x,y);CHKERRQ(ierr);
      67           6 :     ierr = VecAXPY(y,nu,x);CHKERRQ(ierr);
      68             :   }
      69          10 :   PetscFunctionReturn(0);
      70             : }
      71             : 
      72           2 : static PetscErrorCode MatMultTranspose_Cayley(Mat B,Vec x,Vec y)
      73             : {
      74             :   PetscErrorCode ierr;
      75             :   ST             st;
      76             :   ST_CAYLEY      *ctx;
      77             :   PetscScalar    nu;
      78             : 
      79           2 :   PetscFunctionBegin;
      80           2 :   ierr = MatShellGetContext(B,(void**)&st);CHKERRQ(ierr);
      81           2 :   ctx = (ST_CAYLEY*)st->data;
      82           2 :   nu = ctx->nu;
      83             : 
      84           2 :   if (st->shift_matrix == ST_MATMODE_INPLACE) { nu = nu + st->sigma; };
      85           2 :   nu = PetscConj(nu);
      86             : 
      87           2 :   if (st->nmat>1) {
      88             :     /* generalized eigenproblem: y = (A + tB)x */
      89           0 :     ierr = MatMultTranspose(st->A[0],x,y);CHKERRQ(ierr);
      90           0 :     ierr = MatMultTranspose(st->A[1],x,st->work[1]);CHKERRQ(ierr);
      91           0 :     ierr = VecAXPY(y,nu,st->work[1]);CHKERRQ(ierr);
      92             :   } else {
      93             :     /* standard eigenproblem: y = (A + tI)x */
      94           2 :     ierr = MatMultTranspose(st->A[0],x,y);CHKERRQ(ierr);
      95           2 :     ierr = VecAXPY(y,nu,x);CHKERRQ(ierr);
      96             :   }
      97           2 :   PetscFunctionReturn(0);
      98             : }
      99             : 
     100           7 : PetscErrorCode STGetBilinearForm_Cayley(ST st,Mat *B)
     101             : {
     102             :   PetscErrorCode ierr;
     103             : 
     104           7 :   PetscFunctionBegin;
     105           7 :   ierr = STSetUp(st);CHKERRQ(ierr);
     106           7 :   *B = st->T[0];
     107           7 :   ierr = PetscObjectReference((PetscObject)*B);CHKERRQ(ierr);
     108           7 :   PetscFunctionReturn(0);
     109             : }
     110             : 
     111         757 : PetscErrorCode STBackTransform_Cayley(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
     112             : {
     113         757 :   ST_CAYLEY   *ctx = (ST_CAYLEY*)st->data;
     114             :   PetscInt    j;
     115             : #if !defined(PETSC_USE_COMPLEX)
     116             :   PetscScalar t,i,r;
     117             : #endif
     118             : 
     119         757 :   PetscFunctionBegin;
     120             : #if !defined(PETSC_USE_COMPLEX)
     121        2269 :   for (j=0;j<n;j++) {
     122        1512 :     if (eigi[j] == 0.0) eigr[j] = (ctx->nu + eigr[j] * st->sigma) / (eigr[j] - 1.0);
     123             :     else {
     124          98 :       r = eigr[j];
     125          98 :       i = eigi[j];
     126          98 :       r = st->sigma * (r * r + i * i - r) + ctx->nu * (r - 1);
     127          98 :       i = - st->sigma * i - ctx->nu * i;
     128          98 :       t = i * i + r * (r - 2.0) + 1.0;
     129          98 :       eigr[j] = r / t;
     130          98 :       eigi[j] = i / t;
     131             :     }
     132             :   }
     133             : #else
     134             :   for (j=0;j<n;j++) {
     135             :     eigr[j] = (ctx->nu + eigr[j] * st->sigma) / (eigr[j] - 1.0);
     136             :   }
     137             : #endif
     138         757 :   PetscFunctionReturn(0);
     139             : }
     140             : 
     141           5 : PetscErrorCode STPostSolve_Cayley(ST st)
     142             : {
     143             :   PetscErrorCode ierr;
     144             : 
     145           5 :   PetscFunctionBegin;
     146           5 :   if (st->shift_matrix == ST_MATMODE_INPLACE) {
     147           1 :     if (st->nmat>1) {
     148           0 :       ierr = MatAXPY(st->A[0],st->sigma,st->A[1],st->str);CHKERRQ(ierr);
     149             :     } else {
     150           1 :       ierr = MatShift(st->A[0],st->sigma);CHKERRQ(ierr);
     151             :     }
     152           1 :     st->Astate[0] = ((PetscObject)st->A[0])->state;
     153           1 :     st->state = ST_STATE_INITIAL;
     154             :   }
     155           5 :   PetscFunctionReturn(0);
     156             : }
     157             : 
     158          11 : PetscErrorCode STSetUp_Cayley(ST st)
     159             : {
     160             :   PetscErrorCode ierr;
     161             :   PetscInt       n,m;
     162          11 :   ST_CAYLEY      *ctx = (ST_CAYLEY*)st->data;
     163             : 
     164          11 :   PetscFunctionBegin;
     165          11 :   ierr = STSetWorkVecs(st,2);CHKERRQ(ierr);
     166             : 
     167             :   /* if the user did not set the shift, use the target value */
     168          11 :   if (!st->sigma_set) st->sigma = st->defsigma;
     169             : 
     170          11 :   if (!ctx->nu_set) ctx->nu = st->sigma;
     171          11 :   if (ctx->nu == 0.0 && st->sigma == 0.0) SETERRQ(PetscObjectComm((PetscObject)st),1,"Values of shift and antishift cannot be zero simultaneously");
     172             : 
     173             :   /* T[0] = A+nu*B */
     174          11 :   if (st->shift_matrix==ST_MATMODE_INPLACE) {
     175           3 :     ierr = MatGetLocalSize(st->A[0],&n,&m);CHKERRQ(ierr);
     176           3 :     ierr = MatCreateShell(PetscObjectComm((PetscObject)st),n,m,PETSC_DETERMINE,PETSC_DETERMINE,st,&st->T[0]);CHKERRQ(ierr);
     177           3 :     ierr = MatShellSetOperation(st->T[0],MATOP_MULT,(void(*)(void))MatMult_Cayley);CHKERRQ(ierr);
     178           3 :     ierr = MatShellSetOperation(st->T[0],MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Cayley);CHKERRQ(ierr);
     179           3 :     ierr = PetscLogObjectParent((PetscObject)st,(PetscObject)st->T[0]);CHKERRQ(ierr);
     180             :   } else {
     181           8 :     ierr = STMatMAXPY_Private(st,ctx->nu,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),&st->T[0]);CHKERRQ(ierr);
     182             :   }
     183             : 
     184             :   /* T[1] = A-sigma*B */
     185          11 :   ierr = STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),&st->T[1]);CHKERRQ(ierr);
     186          11 :   ierr = PetscObjectReference((PetscObject)st->T[1]);CHKERRQ(ierr);
     187          11 :   ierr = MatDestroy(&st->P);CHKERRQ(ierr);
     188          11 :   st->P = st->T[1];
     189          11 :   if (!st->ksp) { ierr = STGetKSP(st,&st->ksp);CHKERRQ(ierr); }
     190          11 :   ierr = STCheckFactorPackage(st);CHKERRQ(ierr);
     191          11 :   ierr = KSPSetOperators(st->ksp,st->P,st->P);CHKERRQ(ierr);
     192          11 :   ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
     193          11 :   PetscFunctionReturn(0);
     194             : }
     195             : 
     196           6 : PetscErrorCode STSetShift_Cayley(ST st,PetscScalar newshift)
     197             : {
     198             :   PetscErrorCode ierr;
     199           6 :   ST_CAYLEY      *ctx = (ST_CAYLEY*)st->data;
     200             : 
     201           6 :   PetscFunctionBegin;
     202           6 :   if (newshift==0.0 && (!ctx->nu_set || (ctx->nu_set && ctx->nu==0.0))) SETERRQ(PetscObjectComm((PetscObject)st),1,"Values of shift and antishift cannot be zero simultaneously");
     203             : 
     204           6 :   if (!ctx->nu_set) {
     205           6 :     if (st->shift_matrix!=ST_MATMODE_INPLACE) {
     206           4 :       ierr = STMatMAXPY_Private(st,newshift,ctx->nu,0,NULL,PETSC_FALSE,&st->T[0]);CHKERRQ(ierr);
     207             :     }
     208           6 :     ctx->nu = newshift;
     209             :   }
     210           6 :   ierr = STMatMAXPY_Private(st,-newshift,-st->sigma,0,NULL,PETSC_FALSE,&st->T[1]);CHKERRQ(ierr);
     211           6 :   if (st->P!=st->T[1]) {
     212           0 :     ierr = MatDestroy(&st->P);CHKERRQ(ierr);
     213           0 :     st->P = st->T[1];
     214           0 :     ierr = PetscObjectReference((PetscObject)st->P);CHKERRQ(ierr);
     215             :   }
     216           6 :   ierr = KSPSetOperators(st->ksp,st->P,st->P);CHKERRQ(ierr);
     217           6 :   ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
     218           6 :   PetscFunctionReturn(0);
     219             : }
     220             : 
     221           5 : PetscErrorCode STSetFromOptions_Cayley(PetscOptionItems *PetscOptionsObject,ST st)
     222             : {
     223             :   PetscErrorCode ierr;
     224             :   PetscScalar    nu;
     225             :   PetscBool      flg;
     226           5 :   ST_CAYLEY      *ctx = (ST_CAYLEY*)st->data;
     227             : 
     228           5 :   PetscFunctionBegin;
     229           5 :   ierr = PetscOptionsHead(PetscOptionsObject,"ST Cayley Options");CHKERRQ(ierr);
     230             : 
     231           5 :     ierr = PetscOptionsScalar("-st_cayley_antishift","Value of the antishift","STCayleySetAntishift",ctx->nu,&nu,&flg);CHKERRQ(ierr);
     232           5 :     if (flg) { ierr = STCayleySetAntishift(st,nu);CHKERRQ(ierr); }
     233             : 
     234           5 :   ierr = PetscOptionsTail();CHKERRQ(ierr);
     235           5 :   PetscFunctionReturn(0);
     236             : }
     237             : 
     238          10 : static PetscErrorCode STCayleySetAntishift_Cayley(ST st,PetscScalar newshift)
     239             : {
     240             :   PetscErrorCode ierr;
     241          10 :   ST_CAYLEY *ctx = (ST_CAYLEY*)st->data;
     242             : 
     243          10 :   PetscFunctionBegin;
     244          10 :   if (st->state && st->shift_matrix!=ST_MATMODE_INPLACE) {
     245           4 :     ierr = STMatMAXPY_Private(st,newshift,ctx->nu,0,NULL,PETSC_FALSE,&st->T[0]);CHKERRQ(ierr);
     246             :   }
     247          10 :   ctx->nu     = newshift;
     248          10 :   ctx->nu_set = PETSC_TRUE;
     249          10 :   PetscFunctionReturn(0);
     250             : }
     251             : 
     252             : /*@
     253             :    STCayleySetAntishift - Sets the value of the anti-shift for the Cayley
     254             :    spectral transformation.
     255             : 
     256             :    Logically Collective on ST
     257             : 
     258             :    Input Parameters:
     259             : +  st  - the spectral transformation context
     260             : -  nu  - the anti-shift
     261             : 
     262             :    Options Database Key:
     263             : .  -st_cayley_antishift - Sets the value of the anti-shift
     264             : 
     265             :    Level: intermediate
     266             : 
     267             :    Note:
     268             :    In the generalized Cayley transform, the operator can be expressed as
     269             :    OP = inv(A - sigma B)*(A + nu B). This function sets the value of nu.
     270             :    Use STSetShift() for setting sigma.
     271             : 
     272             : .seealso: STSetShift(), STCayleyGetAntishift()
     273             : @*/
     274          10 : PetscErrorCode STCayleySetAntishift(ST st,PetscScalar nu)
     275             : {
     276             :   PetscErrorCode ierr;
     277             : 
     278          10 :   PetscFunctionBegin;
     279          10 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     280          10 :   PetscValidLogicalCollectiveScalar(st,nu,2);
     281          10 :   ierr = PetscTryMethod(st,"STCayleySetAntishift_C",(ST,PetscScalar),(st,nu));CHKERRQ(ierr);
     282          10 :   PetscFunctionReturn(0);
     283             : }
     284             : 
     285          18 : static PetscErrorCode STCayleyGetAntishift_Cayley(ST st,PetscScalar *nu)
     286             : {
     287          18 :   ST_CAYLEY *ctx = (ST_CAYLEY*)st->data;
     288             : 
     289          18 :   PetscFunctionBegin;
     290          18 :   *nu = ctx->nu;
     291          18 :   PetscFunctionReturn(0);
     292             : }
     293             : 
     294             : /*@
     295             :    STCayleyGetAntishift - Gets the value of the anti-shift used in the Cayley
     296             :    spectral transformation.
     297             : 
     298             :    Not Collective
     299             : 
     300             :    Input Parameter:
     301             : .  st  - the spectral transformation context
     302             : 
     303             :    Output Parameter:
     304             : .  nu  - the anti-shift
     305             : 
     306             :    Level: intermediate
     307             : 
     308             : .seealso: STGetShift(), STCayleySetAntishift()
     309             : @*/
     310          18 : PetscErrorCode STCayleyGetAntishift(ST st,PetscScalar *nu)
     311             : {
     312             :   PetscErrorCode ierr;
     313             : 
     314          18 :   PetscFunctionBegin;
     315          18 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     316          18 :   PetscValidScalarPointer(nu,2);
     317          18 :   ierr = PetscUseMethod(st,"STCayleyGetAntishift_C",(ST,PetscScalar*),(st,nu));CHKERRQ(ierr);
     318          18 :   PetscFunctionReturn(0);
     319             : }
     320             : 
     321           0 : PetscErrorCode STView_Cayley(ST st,PetscViewer viewer)
     322             : {
     323             :   PetscErrorCode ierr;
     324             :   char           str[50];
     325           0 :   ST_CAYLEY      *ctx = (ST_CAYLEY*)st->data;
     326             :   PetscBool      isascii;
     327             : 
     328           0 :   PetscFunctionBegin;
     329           0 :   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
     330           0 :   if (isascii) {
     331           0 :     ierr = SlepcSNPrintfScalar(str,50,ctx->nu,PETSC_FALSE);CHKERRQ(ierr);
     332           0 :     ierr = PetscViewerASCIIPrintf(viewer,"  antishift: %s\n",str);CHKERRQ(ierr);
     333             :   }
     334           0 :   PetscFunctionReturn(0);
     335             : }
     336             : 
     337          11 : PetscErrorCode STDestroy_Cayley(ST st)
     338             : {
     339             :   PetscErrorCode ierr;
     340             : 
     341          11 :   PetscFunctionBegin;
     342          11 :   ierr = PetscFree(st->data);CHKERRQ(ierr);
     343          11 :   ierr = PetscObjectComposeFunction((PetscObject)st,"STCayleySetAntishift_C",NULL);CHKERRQ(ierr);
     344          11 :   ierr = PetscObjectComposeFunction((PetscObject)st,"STCayleyGetAntishift_C",NULL);CHKERRQ(ierr);
     345          11 :   PetscFunctionReturn(0);
     346             : }
     347             : 
     348          11 : PETSC_EXTERN PetscErrorCode STCreate_Cayley(ST st)
     349             : {
     350             :   PetscErrorCode ierr;
     351             :   ST_CAYLEY      *ctx;
     352             : 
     353          11 :   PetscFunctionBegin;
     354          11 :   ierr = PetscNewLog(st,&ctx);CHKERRQ(ierr);
     355          11 :   st->data = (void*)ctx;
     356             : 
     357          11 :   st->ops->apply           = STApply_Cayley;
     358          11 :   st->ops->getbilinearform = STGetBilinearForm_Cayley;
     359          11 :   st->ops->applytrans      = STApplyTranspose_Cayley;
     360          11 :   st->ops->postsolve       = STPostSolve_Cayley;
     361          11 :   st->ops->backtransform   = STBackTransform_Cayley;
     362          11 :   st->ops->setfromoptions  = STSetFromOptions_Cayley;
     363          11 :   st->ops->setup           = STSetUp_Cayley;
     364          11 :   st->ops->setshift        = STSetShift_Cayley;
     365          11 :   st->ops->destroy         = STDestroy_Cayley;
     366          11 :   st->ops->view            = STView_Cayley;
     367          11 :   st->ops->checknullspace  = STCheckNullSpace_Default;
     368          11 :   st->ops->setdefaultksp   = STSetDefaultKSP_Default;
     369          11 :   ierr = PetscObjectComposeFunction((PetscObject)st,"STCayleySetAntishift_C",STCayleySetAntishift_Cayley);CHKERRQ(ierr);
     370          11 :   ierr = PetscObjectComposeFunction((PetscObject)st,"STCayleyGetAntishift_C",STCayleyGetAntishift_Cayley);CHKERRQ(ierr);
     371          11 :   PetscFunctionReturn(0);
     372             : }
     373             : 

Generated by: LCOV version 1.13