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             :    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          97 : static PetscErrorCode dvd_improvex_precond_d(dvdDashboard *d)
      45             : {
      46          97 :   dvdPCWrapper   *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;
      47             : 
      48          97 :   PetscFunctionBegin;
      49             :   /* Free local data */
      50          97 :   PetscCall(PCDestroy(&dvdpc->pc));
      51          97 :   PetscCall(PetscFree(d->improvex_precond_data));
      52          97 :   PetscFunctionReturn(PETSC_SUCCESS);
      53             : }
      54             : 
      55       58595 : static PetscErrorCode dvd_static_precond_PC_0(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
      56             : {
      57       58595 :   dvdPCWrapper   *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;
      58             : 
      59       58595 :   PetscFunctionBegin;
      60       58595 :   PetscCall(PCApply(dvdpc->pc,x,Px));
      61       58595 :   PetscFunctionReturn(PETSC_SUCCESS);
      62             : }
      63             : 
      64             : /*
      65             :   Create a trivial preconditioner
      66             : */
      67       10401 : static PetscErrorCode dvd_precond_none(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
      68             : {
      69       10401 :   PetscFunctionBegin;
      70       10401 :   PetscCall(VecCopy(x,Px));
      71       10401 :   PetscFunctionReturn(PETSC_SUCCESS);
      72             : }
      73             : 
      74             : /*
      75             :   Create a static preconditioner from a PC
      76             : */
      77         291 : PetscErrorCode dvd_static_precond_PC(dvdDashboard *d,dvdBlackboard *b,PC pc)
      78             : {
      79         291 :   dvdPCWrapper   *dvdpc;
      80         291 :   Mat            P;
      81         291 :   PetscBool      t0,t1,t2;
      82             : 
      83         291 :   PetscFunctionBegin;
      84             :   /* Setup the step */
      85         291 :   if (b->state >= DVD_STATE_CONF) {
      86             :     /* If the preconditioner is valid */
      87          97 :     if (pc) {
      88          97 :       PetscCall(PetscNew(&dvdpc));
      89          97 :       dvdpc->pc = pc;
      90          97 :       PetscCall(PetscObjectReference((PetscObject)pc));
      91          97 :       d->improvex_precond_data = dvdpc;
      92          97 :       d->improvex_precond = dvd_static_precond_PC_0;
      93             : 
      94             :       /* PC saves the matrix associated with the linear system, and it has to
      95             :          be initialize to a valid matrix */
      96          97 :       PetscCall(PCGetOperatorsSet(pc,NULL,&t0));
      97          97 :       PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCNONE,&t1));
      98          97 :       PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&t2));
      99          97 :       if (t0 && !t1) {
     100          41 :         PetscCall(PCGetOperators(pc,NULL,&P));
     101          41 :         PetscCall(PetscObjectReference((PetscObject)P));
     102          41 :         PetscCall(PCSetOperators(pc,P,P));
     103          41 :         PetscCall(PCSetReusePreconditioner(pc,PETSC_TRUE));
     104          41 :         PetscCall(MatDestroy(&P));
     105          56 :       } else if (t2) {
     106           0 :         PetscCall(PCSetOperators(pc,d->A,d->A));
     107           0 :         PetscCall(PCSetReusePreconditioner(pc,PETSC_TRUE));
     108             :       } else {
     109          56 :         d->improvex_precond = dvd_precond_none;
     110             :       }
     111             : 
     112          97 :       PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_precond_d));
     113             : 
     114             :     /* Else, use no preconditioner */
     115           0 :     } else d->improvex_precond = dvd_precond_none;
     116             :   }
     117         291 :   PetscFunctionReturn(PETSC_SUCCESS);
     118             : }
     119             : 
     120          23 : static PetscErrorCode dvd_harm_d(dvdDashboard *d)
     121             : {
     122          23 :   PetscFunctionBegin;
     123             :   /* Free local data */
     124          23 :   PetscCall(PetscFree(d->calcpairs_W_data));
     125          23 :   PetscFunctionReturn(PETSC_SUCCESS);
     126             : }
     127             : 
     128          23 : static PetscErrorCode dvd_harm_transf(dvdHarmonic *dvdh,PetscScalar t)
     129             : {
     130          23 :   PetscFunctionBegin;
     131          23 :   switch (dvdh->mode) {
     132          13 :   case DVD_HARM_RR:    /* harmonic RR */
     133          13 :     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 0.0; dvdh->Pb = -1.0;
     134          13 :     break;
     135           0 :   case DVD_HARM_RRR:   /* relative harmonic RR */
     136           0 :     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 1.0; dvdh->Pb = 0.0;
     137           0 :     break;
     138           0 :   case DVD_HARM_REIGS: /* rightmost eigenvalues */
     139           0 :     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 1.0; dvdh->Pb = -PetscConj(t);
     140           0 :     break;
     141          10 :   case DVD_HARM_LEIGS: /* largest eigenvalues */
     142          10 :     dvdh->Wa = 0.0; dvdh->Wb = 1.0; dvdh->Pa = 1.0; dvdh->Pb = 0.0;
     143          10 :     break;
     144           0 :   case DVD_HARM_NONE:
     145             :   default:
     146           0 :     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Harmonic type not supported");
     147             :   }
     148             : 
     149             :   /* Check the transformation does not change the sign of the imaginary part */
     150             : #if !defined(PETSC_USE_COMPLEX)
     151             :   if (dvdh->Pb*dvdh->Wa - dvdh->Wb*dvdh->Pa < 0.0) {
     152             :     dvdh->Pa *= -1.0;
     153             :     dvdh->Pb *= -1.0;
     154             :   }
     155             : #endif
     156          23 :   PetscFunctionReturn(PETSC_SUCCESS);
     157             : }
     158             : 
     159        1300 : static PetscErrorCode dvd_harm_updateW(dvdDashboard *d)
     160             : {
     161        1300 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     162        1300 :   PetscInt       l,k;
     163        1300 :   BV             BX = d->BX?d->BX:d->eps->V;
     164             : 
     165        1300 :   PetscFunctionBegin;
     166             :   /* Update the target if it is necessary */
     167        1300 :   if (!data->withTarget) PetscCall(dvd_harm_transf(data,d->eigr[0]));
     168             : 
     169             :   /* W(i) <- Wa*AV(i) - Wb*BV(i) */
     170        1300 :   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
     171        1300 :   PetscAssert(k==l+d->V_new_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
     172        1300 :   PetscCall(BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e));
     173        1300 :   PetscCall(BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e));
     174        1300 :   PetscCall(BVSetActiveColumns(BX,l+d->V_new_s,l+d->V_new_e));
     175        1300 :   PetscCall(BVCopy(d->AX,d->W));
     176        1300 :   PetscCall(BVScale(d->W,data->Wa));
     177        1300 :   PetscCall(BVMult(d->W,-data->Wb,1.0,BX,NULL));
     178        1300 :   PetscCall(BVSetActiveColumns(d->W,l,k));
     179        1300 :   PetscCall(BVSetActiveColumns(d->AX,l,k));
     180        1300 :   PetscCall(BVSetActiveColumns(BX,l,k));
     181        1300 :   PetscFunctionReturn(PETSC_SUCCESS);
     182             : }
     183             : 
     184        1300 : static PetscErrorCode dvd_harm_proj(dvdDashboard *d)
     185             : {
     186        1300 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     187        1300 :   PetscInt       i,j,l0,l,k,ld;
     188        1300 :   PetscScalar    h,g,*H,*G;
     189             : 
     190        1300 :   PetscFunctionBegin;
     191        1300 :   PetscCall(BVGetActiveColumns(d->eps->V,&l0,&k));
     192        1300 :   l = l0 + d->V_new_s;
     193        1300 :   k = l0 + d->V_new_e;
     194        1300 :   PetscCall(MatGetSize(d->H,&ld,NULL));
     195        1300 :   PetscCall(MatDenseGetArray(d->H,&H));
     196        1300 :   PetscCall(MatDenseGetArray(d->G,&G));
     197             :   /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
     198             :   /* Right part */
     199        3113 :   for (i=l;i<k;i++) {
     200       23137 :     for (j=l0;j<k;j++) {
     201       21324 :       h = H[ld*i+j];
     202       21324 :       g = G[ld*i+j];
     203       21324 :       H[ld*i+j] = data->Pa*h - data->Pb*g;
     204       21324 :       G[ld*i+j] = data->Wa*h - data->Wb*g;
     205             :     }
     206             :   }
     207             :   /* Left part */
     208       14966 :   for (i=l0;i<l;i++) {
     209       31691 :     for (j=l;j<k;j++) {
     210       18025 :       h = H[ld*i+j];
     211       18025 :       g = G[ld*i+j];
     212       18025 :       H[ld*i+j] = data->Pa*h - data->Pb*g;
     213       18025 :       G[ld*i+j] = data->Wa*h - data->Wb*g;
     214             :     }
     215             :   }
     216        1300 :   PetscCall(MatDenseRestoreArray(d->H,&H));
     217        1300 :   PetscCall(MatDenseRestoreArray(d->G,&G));
     218        1300 :   PetscFunctionReturn(PETSC_SUCCESS);
     219             : }
     220             : 
     221          57 : PetscErrorCode dvd_harm_updateproj(dvdDashboard *d)
     222             : {
     223          57 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     224          57 :   PetscInt       i,j,l,k,ld;
     225          57 :   PetscScalar    h,g,*H,*G;
     226             : 
     227          57 :   PetscFunctionBegin;
     228          57 :   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
     229          57 :   k = l + d->V_tra_s;
     230          57 :   PetscCall(MatGetSize(d->H,&ld,NULL));
     231          57 :   PetscCall(MatDenseGetArray(d->H,&H));
     232          57 :   PetscCall(MatDenseGetArray(d->G,&G));
     233             :   /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
     234             :   /* Right part */
     235         114 :   for (i=l;i<k;i++) {
     236         135 :     for (j=0;j<l;j++) {
     237          78 :       h = H[ld*i+j];
     238          78 :       g = G[ld*i+j];
     239          78 :       H[ld*i+j] = data->Pa*h - data->Pb*g;
     240          78 :       G[ld*i+j] = data->Wa*h - data->Wb*g;
     241             :     }
     242             :   }
     243             :   /* Lower triangular part */
     244         135 :   for (i=0;i<l;i++) {
     245         156 :     for (j=l;j<k;j++) {
     246          78 :       h = H[ld*i+j];
     247          78 :       g = G[ld*i+j];
     248          78 :       H[ld*i+j] = data->Pa*h - data->Pb*g;
     249          78 :       G[ld*i+j] = data->Wa*h - data->Wb*g;
     250             :     }
     251             :   }
     252          57 :   PetscCall(MatDenseRestoreArray(d->H,&H));
     253          57 :   PetscCall(MatDenseRestoreArray(d->G,&G));
     254          57 :   PetscFunctionReturn(PETSC_SUCCESS);
     255             : }
     256             : 
     257       51207 : static PetscErrorCode dvd_harm_backtrans(dvdHarmonic *data,PetscScalar *ar,PetscScalar *ai)
     258             : {
     259       51207 :   PetscScalar xr;
     260             : #if !defined(PETSC_USE_COMPLEX)
     261             :   PetscScalar xi, k;
     262             : #endif
     263             : 
     264       51207 :   PetscFunctionBegin;
     265       51207 :   xr = *ar;
     266             : #if !defined(PETSC_USE_COMPLEX)
     267             :   xi = *ai;
     268             :   if (PetscUnlikely(xi != 0.0)) {
     269             :     k = (data->Pa - data->Wa*xr)*(data->Pa - data->Wa*xr) + data->Wa*data->Wa*xi*xi;
     270             :     *ar = (data->Pb*data->Pa - (data->Pb*data->Wa + data->Wb*data->Pa)*xr + data->Wb*data->Wa*(xr*xr + xi*xi))/k;
     271             :     *ai = (data->Pb*data->Wa - data->Wb*data->Pa)*xi/k;
     272             :   } else
     273             : #endif
     274       51207 :     *ar = (data->Pb - data->Wb*xr) / (data->Pa - data->Wa*xr);
     275       51207 :   PetscFunctionReturn(PETSC_SUCCESS);
     276             : }
     277             : 
     278       34157 : static PetscErrorCode dvd_harm_eig_backtrans(dvdDashboard *d,PetscScalar ar,PetscScalar ai,PetscScalar *br,PetscScalar *bi)
     279             : {
     280       34157 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     281             : 
     282       34157 :   PetscFunctionBegin;
     283       34157 :   PetscCall(dvd_harm_backtrans(data,&ar,&ai));
     284       34157 :   *br = ar;
     285       34157 :   *bi = ai;
     286       34157 :   PetscFunctionReturn(PETSC_SUCCESS);
     287             : }
     288             : 
     289        1513 : static PetscErrorCode dvd_harm_eigs_trans(dvdDashboard *d)
     290             : {
     291        1513 :   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
     292        1513 :   PetscInt       i,l,k;
     293             : 
     294        1513 :   PetscFunctionBegin;
     295        1513 :   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
     296       18563 :   for (i=0;i<k-l;i++) PetscCall(dvd_harm_backtrans(data,&d->eigr[i],&d->eigi[i]));
     297        1513 :   PetscFunctionReturn(PETSC_SUCCESS);
     298             : }
     299             : 
     300          69 : PetscErrorCode dvd_harm_conf(dvdDashboard *d,dvdBlackboard *b,HarmType_t mode,PetscBool fixedTarget,PetscScalar t)
     301             : {
     302          69 :   dvdHarmonic    *dvdh;
     303             : 
     304          69 :   PetscFunctionBegin;
     305             :   /* Set the problem to GNHEP:
     306             :      d->G maybe is upper triangular due to biorthogonality of V and W */
     307          69 :   d->sEP = d->sA = d->sB = 0;
     308             : 
     309             :   /* Setup the step */
     310          69 :   if (b->state >= DVD_STATE_CONF) {
     311          23 :     PetscCall(PetscNew(&dvdh));
     312          23 :     dvdh->withTarget = fixedTarget;
     313          23 :     dvdh->mode = mode;
     314          23 :     if (fixedTarget) PetscCall(dvd_harm_transf(dvdh, t));
     315          23 :     d->calcpairs_W_data = dvdh;
     316          23 :     d->calcpairs_W = dvd_harm_updateW;
     317          23 :     d->calcpairs_proj_trans = dvd_harm_proj;
     318          23 :     d->calcpairs_eigs_trans = dvd_harm_eigs_trans;
     319          23 :     d->calcpairs_eig_backtrans = dvd_harm_eig_backtrans;
     320             : 
     321          23 :     PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_harm_d));
     322             :   }
     323          69 :   PetscFunctionReturn(PETSC_SUCCESS);
     324             : }