LCOV - code coverage report
Current view: top level - eps/impls/davidson - dvdutils.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 157 166 94.6 %
Date: 2019-07-19 01:26:51 Functions: 13 13 100.0 %

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2019, 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             :    SLEPc eigensolver: "davidson"
      12             : 
      13             :    Some utils
      14             : */
      15             : 
      16             : #include "davidson.h"
      17             : 
      18             : typedef struct {
      19             :   PC pc;
      20             : } dvdPCWrapper;
      21             : 
      22             : /*
      23             :   Configure the harmonics.
      24             :   switch (mode) {
      25             :   DVD_HARM_RR:    harmonic RR
      26             :   DVD_HARM_RRR:   relative harmonic RR
      27             :   DVD_HARM_REIGS: rightmost eigenvalues
      28             :   DVD_HARM_LEIGS: largest eigenvalues
      29             :   }
      30             :   fixedTarged, if true use the target instead of the best eigenvalue
      31             :   target, the fixed target to be used
      32             : */
      33             : typedef struct {
      34             :   PetscScalar Wa,Wb;       /* span{W} = span{Wa*AV - Wb*BV} */
      35             :   PetscScalar Pa,Pb;       /* H=W'*(Pa*AV - Pb*BV), G=W'*(Wa*AV - Wb*BV) */
      36             :   PetscBool   withTarget;
      37             :   HarmType_t  mode;
      38             : } dvdHarmonic;
      39             : 
      40             : typedef struct {
      41             :   Vec diagA, diagB;
      42             : } dvdJacobiPrecond;
      43             : 
      44          90 : static PetscErrorCode dvd_improvex_precond_d(dvdDashboard *d)
      45             : {
      46             :   PetscErrorCode ierr;
      47          90 :   dvdPCWrapper   *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;
      48             : 
      49          90 :   PetscFunctionBegin;
      50             :   /* Free local data */
      51          90 :   ierr = PCDestroy(&dvdpc->pc);CHKERRQ(ierr);
      52          90 :   ierr = PetscFree(d->improvex_precond_data);CHKERRQ(ierr);
      53          90 :   PetscFunctionReturn(0);
      54             : }
      55             : 
      56       27656 : static PetscErrorCode dvd_static_precond_PC_0(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
      57             : {
      58             :   PetscErrorCode ierr;
      59       27656 :   dvdPCWrapper   *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;
      60             : 
      61       27656 :   PetscFunctionBegin;
      62       27656 :   ierr = PCApply(dvdpc->pc,x,Px);CHKERRQ(ierr);
      63       27656 :   PetscFunctionReturn(0);
      64             : }
      65             : 
      66             : /*
      67             :   Create a trivial preconditioner
      68             : */
      69        7217 : static PetscErrorCode dvd_precond_none(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
      70             : {
      71             :   PetscErrorCode ierr;
      72             : 
      73        7217 :   PetscFunctionBegin;
      74        7217 :   ierr = VecCopy(x,Px);CHKERRQ(ierr);
      75        7217 :   PetscFunctionReturn(0);
      76             : }
      77             : 
      78             : /*
      79             :   Create a static preconditioner from a PC
      80             : */
      81         270 : PetscErrorCode dvd_static_precond_PC(dvdDashboard *d,dvdBlackboard *b,PC pc)
      82             : {
      83             :   PetscErrorCode ierr;
      84             :   dvdPCWrapper   *dvdpc;
      85             :   Mat            P;
      86             :   PetscBool      t0,t1,t2;
      87             : 
      88         270 :   PetscFunctionBegin;
      89             :   /* Setup the step */
      90         270 :   if (b->state >= DVD_STATE_CONF) {
      91             :     /* If the preconditioner is valid */
      92          90 :     if (pc) {
      93          90 :       ierr = PetscNewLog(d->eps,&dvdpc);CHKERRQ(ierr);
      94          90 :       dvdpc->pc = pc;
      95          90 :       ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
      96          90 :       d->improvex_precond_data = dvdpc;
      97          90 :       d->improvex_precond = dvd_static_precond_PC_0;
      98             : 
      99             :       /* PC saves the matrix associated with the linear system, and it has to
     100             :          be initialize to a valid matrix */
     101          90 :       ierr = PCGetOperatorsSet(pc,NULL,&t0);CHKERRQ(ierr);
     102          90 :       ierr = PetscObjectTypeCompare((PetscObject)pc,PCNONE,&t1);CHKERRQ(ierr);
     103          90 :       ierr = PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&t2);CHKERRQ(ierr);
     104          90 :       if (t0 && !t1) {
     105          48 :         ierr = PCGetOperators(pc,NULL,&P);CHKERRQ(ierr);
     106          48 :         ierr = PetscObjectReference((PetscObject)P);CHKERRQ(ierr);
     107          48 :         ierr = PCSetOperators(pc,P,P);CHKERRQ(ierr);
     108          48 :         ierr = PCSetReusePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr);
     109          48 :         ierr = MatDestroy(&P);CHKERRQ(ierr);
     110          42 :       } else if (t2) {
     111           0 :         ierr = PCSetOperators(pc,d->A,d->A);CHKERRQ(ierr);
     112           0 :         ierr = PCSetReusePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr);
     113             :       } else {
     114          42 :         d->improvex_precond = dvd_precond_none;
     115             :       }
     116             : 
     117          90 :       ierr = EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_precond_d);CHKERRQ(ierr);
     118             : 
     119             :     /* Else, use no preconditioner */
     120           0 :     } else d->improvex_precond = dvd_precond_none;
     121             :   }
     122         270 :   PetscFunctionReturn(0);
     123             : }
     124             : 
     125          21 : static PetscErrorCode dvd_harm_d(dvdDashboard *d)
     126             : {
     127             :   PetscErrorCode ierr;
     128             : 
     129          21 :   PetscFunctionBegin;
     130             :   /* Free local data */
     131          21 :   ierr = PetscFree(d->calcpairs_W_data);CHKERRQ(ierr);
     132          21 :   PetscFunctionReturn(0);
     133             : }
     134             : 
     135          21 : static PetscErrorCode dvd_harm_transf(dvdHarmonic *dvdh,PetscScalar t)
     136             : {
     137          21 :   PetscFunctionBegin;
     138          21 :   switch (dvdh->mode) {
     139             :   case DVD_HARM_RR:    /* harmonic RR */
     140          13 :     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 0.0; dvdh->Pb = -1.0;
     141          13 :     break;
     142             :   case DVD_HARM_RRR:   /* relative harmonic RR */
     143           0 :     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 1.0; dvdh->Pb = 0.0;
     144           0 :     break;
     145             :   case DVD_HARM_REIGS: /* rightmost eigenvalues */
     146           0 :     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 1.0; dvdh->Pb = -PetscConj(t);
     147           0 :     break;
     148             :   case DVD_HARM_LEIGS: /* largest eigenvalues */
     149           8 :     dvdh->Wa = 0.0; dvdh->Wb = 1.0; dvdh->Pa = 1.0; dvdh->Pb = 0.0;
     150           8 :     break;
     151             :   case DVD_HARM_NONE:
     152             :   default:
     153           0 :     SETERRQ(PETSC_COMM_SELF,1,"Harmonic type not supported");
     154             :   }
     155             : 
     156             :   /* Check the transformation does not change the sign of the imaginary part */
     157             : #if !defined(PETSC_USE_COMPLEX)
     158          21 :   if (dvdh->Pb*dvdh->Wa - dvdh->Wb*dvdh->Pa < 0.0) {
     159          21 :     dvdh->Pa *= -1.0;
     160          21 :     dvdh->Pb *= -1.0;
     161             :   }
     162             : #endif
     163          21 :   PetscFunctionReturn(0);
     164             : }
     165             : 
     166        1081 : static PetscErrorCode dvd_harm_updateW(dvdDashboard *d)
     167             : {
     168        1081 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     169             :   PetscErrorCode ierr;
     170             :   PetscInt       l,k;
     171        1081 :   BV             BX = d->BX?d->BX:d->eps->V;
     172             : 
     173        1081 :   PetscFunctionBegin;
     174             :   /* Update the target if it is necessary */
     175        1081 :   if (!data->withTarget) {
     176           0 :     ierr = dvd_harm_transf(data,d->eigr[0]);CHKERRQ(ierr);
     177             :   }
     178             : 
     179             :   /* W(i) <- Wa*AV(i) - Wb*BV(i) */
     180        1081 :   ierr = BVGetActiveColumns(d->eps->V,&l,&k);CHKERRQ(ierr);
     181        1081 :   if (k != l+d->V_new_s) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
     182        1081 :   ierr = BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr);
     183        1081 :   ierr = BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr);
     184        1081 :   ierr = BVSetActiveColumns(BX,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr);
     185        1081 :   ierr = BVCopy(d->AX,d->W);CHKERRQ(ierr);
     186        1081 :   ierr = BVScale(d->W,data->Wa);CHKERRQ(ierr);
     187        1081 :   ierr = BVMult(d->W,-data->Wb,1.0,BX,NULL);CHKERRQ(ierr);
     188        1081 :   ierr = BVSetActiveColumns(d->W,l,k);CHKERRQ(ierr);
     189        1081 :   ierr = BVSetActiveColumns(d->AX,l,k);CHKERRQ(ierr);
     190        1081 :   ierr = BVSetActiveColumns(BX,l,k);CHKERRQ(ierr);
     191        1081 :   PetscFunctionReturn(0);
     192             : }
     193             : 
     194        1081 : static PetscErrorCode dvd_harm_proj(dvdDashboard *d)
     195             : {
     196             :   PetscErrorCode ierr;
     197        1081 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     198             :   PetscInt       i,j,l0,l,k,ld;
     199             :   PetscScalar    h,g,*H,*G;
     200             : 
     201        1081 :   PetscFunctionBegin;
     202        1081 :   ierr = BVGetActiveColumns(d->eps->V,&l0,&k);CHKERRQ(ierr);
     203        1081 :   l = l0 + d->V_new_s;
     204        1081 :   k = l0 + d->V_new_e;
     205        1081 :   ierr = MatGetSize(d->H,&ld,NULL);CHKERRQ(ierr);
     206        1081 :   ierr = MatDenseGetArray(d->H,&H);CHKERRQ(ierr);
     207        1081 :   ierr = MatDenseGetArray(d->G,&G);CHKERRQ(ierr);
     208             :   /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
     209             :   /* Right part */
     210        1582 :   for (i=l;i<k;i++) {
     211       19829 :     for (j=l0;j<k;j++) {
     212       18247 :       h = H[ld*i+j];
     213       18247 :       g = G[ld*i+j];
     214       18247 :       H[ld*i+j] = data->Pa*h - data->Pb*g;
     215       18247 :       G[ld*i+j] = data->Wa*h - data->Wb*g;
     216             :     }
     217             :   }
     218             :   /* Left part */
     219       11977 :   for (i=l0;i<l;i++) {
     220       15243 :     for (j=l;j<k;j++) {
     221       15243 :       h = H[ld*i+j];
     222       15243 :       g = G[ld*i+j];
     223       15243 :       H[ld*i+j] = data->Pa*h - data->Pb*g;
     224       15243 :       G[ld*i+j] = data->Wa*h - data->Wb*g;
     225             :     }
     226             :   }
     227        1081 :   ierr = MatDenseRestoreArray(d->H,&H);CHKERRQ(ierr);
     228        1081 :   ierr = MatDenseRestoreArray(d->G,&G);CHKERRQ(ierr);
     229        1081 :   PetscFunctionReturn(0);
     230             : }
     231             : 
     232          57 : PetscErrorCode dvd_harm_updateproj(dvdDashboard *d)
     233             : {
     234             :   PetscErrorCode ierr;
     235          57 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     236             :   PetscInt       i,j,l,k,ld;
     237             :   PetscScalar    h,g,*H,*G;
     238             : 
     239          57 :   PetscFunctionBegin;
     240          57 :   ierr = BVGetActiveColumns(d->eps->V,&l,&k);CHKERRQ(ierr);
     241          57 :   k = l + d->V_tra_s;
     242          57 :   ierr = MatGetSize(d->H,&ld,NULL);CHKERRQ(ierr);
     243          57 :   ierr = MatDenseGetArray(d->H,&H);CHKERRQ(ierr);
     244          57 :   ierr = MatDenseGetArray(d->G,&G);CHKERRQ(ierr);
     245             :   /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
     246             :   /* Right part */
     247         116 :   for (i=l;i<k;i++) {
     248          86 :     for (j=0;j<l;j++) {
     249          86 :       h = H[ld*i+j];
     250          86 :       g = G[ld*i+j];
     251          86 :       H[ld*i+j] = data->Pa*h - data->Pb*g;
     252          86 :       G[ld*i+j] = data->Wa*h - data->Wb*g;
     253             :     }
     254             :   }
     255             :   /* Lower triangular part */
     256          84 :   for (i=0;i<l;i++) {
     257          86 :     for (j=l;j<k;j++) {
     258          86 :       h = H[ld*i+j];
     259          86 :       g = G[ld*i+j];
     260          86 :       H[ld*i+j] = data->Pa*h - data->Pb*g;
     261          86 :       G[ld*i+j] = data->Wa*h - data->Wb*g;
     262             :     }
     263             :   }
     264          57 :   ierr = MatDenseRestoreArray(d->H,&H);CHKERRQ(ierr);
     265          57 :   ierr = MatDenseRestoreArray(d->G,&G);CHKERRQ(ierr);
     266          57 :   PetscFunctionReturn(0);
     267             : }
     268             : 
     269       42224 : static PetscErrorCode dvd_harm_backtrans(dvdHarmonic *data,PetscScalar *ar,PetscScalar *ai)
     270             : {
     271             :   PetscScalar xr;
     272             : #if !defined(PETSC_USE_COMPLEX)
     273             :   PetscScalar xi, k;
     274             : #endif
     275             : 
     276       42224 :   PetscFunctionBegin;
     277       42224 :   PetscValidPointer(ar,2);
     278       42224 :   xr = *ar;
     279             : #if !defined(PETSC_USE_COMPLEX)
     280       42224 :   PetscValidPointer(ai,3);
     281       42224 :   xi = *ai;
     282             : 
     283       42224 :   if (xi != 0.0) {
     284        2500 :     k = (data->Pa - data->Wa*xr)*(data->Pa - data->Wa*xr) + data->Wa*data->Wa*xi*xi;
     285        2500 :     *ar = (data->Pb*data->Pa - (data->Pb*data->Wa + data->Wb*data->Pa)*xr + data->Wb*data->Wa*(xr*xr + xi*xi))/k;
     286        2500 :     *ai = (data->Pb*data->Wa - data->Wb*data->Pa)*xi/k;
     287             :   } else
     288             : #endif
     289       39724 :     *ar = (data->Pb - data->Wb*xr) / (data->Pa - data->Wa*xr);
     290       42224 :   PetscFunctionReturn(0);
     291             : }
     292             : 
     293       28169 : static PetscErrorCode dvd_harm_eig_backtrans(dvdDashboard *d,PetscScalar ar,PetscScalar ai,PetscScalar *br,PetscScalar *bi)
     294             : {
     295       28169 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     296             :   PetscErrorCode ierr;
     297             : 
     298       28169 :   PetscFunctionBegin;
     299       28169 :   ierr = dvd_harm_backtrans(data,&ar,&ai);CHKERRQ(ierr);
     300       28169 :   *br = ar;
     301       28169 :   *bi = ai;
     302       28169 :   PetscFunctionReturn(0);
     303             : }
     304             : 
     305        1295 : static PetscErrorCode dvd_harm_eigs_trans(dvdDashboard *d)
     306             : {
     307        1295 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     308             :   PetscInt       i,l,k;
     309             :   PetscErrorCode ierr;
     310             : 
     311        1295 :   PetscFunctionBegin;
     312        1295 :   ierr = BVGetActiveColumns(d->eps->V,&l,&k);CHKERRQ(ierr);
     313       14055 :   for (i=0;i<k-l;i++) {
     314       14055 :     ierr = dvd_harm_backtrans(data,&d->eigr[i],&d->eigi[i]);CHKERRQ(ierr);
     315             :   }
     316        1295 :   PetscFunctionReturn(0);
     317             : }
     318             : 
     319          63 : PetscErrorCode dvd_harm_conf(dvdDashboard *d,dvdBlackboard *b,HarmType_t mode,PetscBool fixedTarget,PetscScalar t)
     320             : {
     321             :   PetscErrorCode ierr;
     322             :   dvdHarmonic    *dvdh;
     323             : 
     324          63 :   PetscFunctionBegin;
     325             :   /* Set the problem to GNHEP:
     326             :      d->G maybe is upper triangular due to biorthogonality of V and W */
     327          63 :   d->sEP = d->sA = d->sB = 0;
     328             : 
     329             :   /* Setup the step */
     330          63 :   if (b->state >= DVD_STATE_CONF) {
     331          21 :     ierr = PetscNewLog(d->eps,&dvdh);CHKERRQ(ierr);
     332          21 :     dvdh->withTarget = fixedTarget;
     333          21 :     dvdh->mode = mode;
     334          21 :     if (fixedTarget) dvd_harm_transf(dvdh, t);
     335          21 :     d->calcpairs_W_data = dvdh;
     336          21 :     d->calcpairs_W = dvd_harm_updateW;
     337          21 :     d->calcpairs_proj_trans = dvd_harm_proj;
     338          21 :     d->calcpairs_eigs_trans = dvd_harm_eigs_trans;
     339          21 :     d->calcpairs_eig_backtrans = dvd_harm_eig_backtrans;
     340             : 
     341          21 :     ierr = EPSDavidsonFLAdd(&d->destroyList,dvd_harm_d);CHKERRQ(ierr);
     342             :   }
     343          63 :   PetscFunctionReturn(0);
     344             : }
     345             : 

Generated by: LCOV version 1.13