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             :    HZ iteration for generalized symmetric-indefinite eigenproblem.
      12             :    Based on Matlab code from David Watkins.
      13             : 
      14             :    References:
      15             : 
      16             :        [1] D.S. Watkins, The Matrix Eigenvalue Problem: GR and Krylov Subspace
      17             :            Methods, SIAM, 2007.
      18             : 
      19             :        [2] M.A. Brebner, J. Grad, "Eigenvalues of Ax = lambda Bx for real
      20             :            symmetric matrices A and B computed by reduction to pseudosymmetric
      21             :            form and the HR process", Linear Alg. Appl. 43:99-118, 1982.
      22             : */
      23             : 
      24             : #include <slepc/private/dsimpl.h>
      25             : #include <slepcblaslapack.h>
      26             : 
      27             : /*
      28             :    Sets up a 2-by-2 matrix to eliminate y in the vector [x y]'.
      29             :    Transformation is rotator if sygn = 1 and hyperbolic if sygn = -1.
      30             : */
      31         355 : static PetscErrorCode UnifiedRotation(PetscReal x,PetscReal y,PetscReal sygn,PetscReal *rot,PetscReal *rcond,PetscBool *swap)
      32             : {
      33         355 :   PetscReal nrm,c,s;
      34             : 
      35         355 :   PetscFunctionBegin;
      36         355 :   *swap = PETSC_FALSE;
      37         355 :   if (y == 0) {
      38           0 :     rot[0] = 1.0; rot[1] = 0.0; rot[2] = 0.0; rot[3] = 1.0;
      39           0 :     *rcond = 1.0;
      40             :   } else {
      41         355 :     nrm = PetscMax(PetscAbs(x),PetscAbs(y));
      42         355 :     c = x/nrm; s = y/nrm;
      43         355 :     PetscCheck(sygn==1.0 || sygn==-1,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Value of sygn sent to transetup must be 1 or -1");
      44         355 :     if (sygn == 1.0) {  /* set up a rotator */
      45         239 :       nrm = SlepcAbs(c,s);
      46         239 :       c = c/nrm; s = s/nrm;
      47             :       /* rot = [c s; -s c]; */
      48         239 :       rot[0] = c; rot[1] = -s; rot[2] = s; rot[3] = c;
      49         239 :       *rcond = 1.0;
      50             :     } else {  /* sygn == -1, set up a hyperbolic transformation */
      51         116 :       nrm = c*c-s*s;
      52         116 :       if (nrm > 0) nrm = PetscSqrtReal(nrm);
      53             :       else {
      54          59 :         PetscCheck(nrm<0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Breakdown in construction of hyperbolic transformation");
      55          59 :         nrm = PetscSqrtReal(-nrm);
      56          59 :         *swap = PETSC_TRUE;
      57             :       }
      58         116 :       c = c/nrm; s = s/nrm;
      59             :       /* rot = [c -s; -s c]; */
      60         116 :       rot[0] = c; rot[1] = -s; rot[2] = -s; rot[3] = c;
      61         121 :       *rcond = PetscAbs(PetscAbs(s)-PetscAbs(c))/(PetscAbs(s)+PetscAbs(c));
      62             :     }
      63             :   }
      64         355 :   PetscFunctionReturn(PETSC_SUCCESS);
      65             : }
      66             : 
      67          37 : static PetscErrorCode HZStep(PetscBLASInt ntop,PetscBLASInt nn,PetscReal tr,PetscReal dt,PetscReal *aa,PetscReal *bb,PetscReal *dd,PetscScalar *uu,PetscInt n,PetscInt ld,PetscBool *flag)
      68             : {
      69          37 :   PetscBLASInt   one=1;
      70          37 :   PetscInt       k,jj,ii;
      71          37 :   PetscBLASInt   n_;
      72          37 :   PetscReal      bulge10,bulge20,bulge30,bulge31,bulge41,bulge42;
      73          37 :   PetscReal      sygn,rcond=1.0,worstcond,rot[4],buf[2],t;
      74          37 :   PetscScalar    rtmp;
      75          37 :   PetscBool      swap;
      76             : 
      77          37 :   PetscFunctionBegin;
      78          37 :   *flag = PETSC_FALSE;
      79          37 :   worstcond = 1.0;
      80          37 :   PetscCall(PetscBLASIntCast(n,&n_));
      81             : 
      82             :   /* Build initial bulge that sets step in motion */
      83          37 :   bulge10 = dd[ntop+1]*(aa[ntop]*(aa[ntop] - dd[ntop]*tr) + dt*dd[ntop]*dd[ntop]) + dd[ntop]*bb[ntop]*bb[ntop];
      84          37 :   bulge20 = bb[ntop]*(dd[ntop+1]*aa[ntop] + dd[ntop]*aa[ntop+1] - dd[ntop]*dd[ntop+1]*tr);
      85          37 :   bulge30 = bb[ntop]*bb[ntop+1]*dd[ntop];
      86          37 :   bulge31 = 0.0;
      87          37 :   bulge41 = 0.0;
      88          37 :   bulge42 = 0.0;
      89             : 
      90             :   /* Chase the bulge */
      91         237 :   for (jj=ntop;jj<nn-1;jj++) {
      92             : 
      93             :     /* Check for trivial bulge */
      94         487 :     if (jj>ntop && PetscMax(PetscMax(PetscAbs(bulge10),PetscAbs(bulge20)),PetscAbs(bulge30))<PETSC_MACHINE_EPSILON*(PetscAbs(aa[jj]) + PetscAbs(aa[jj+1]))) {
      95           8 :       bb[jj-1] = 0.0;  /* deflate and move on */
      96             : 
      97             :     } else { /* carry out the step */
      98             : 
      99             :       /* Annihilate tip entry bulge30 */
     100         192 :       if (bulge30 != 0.0) {
     101             : 
     102             :         /* Make an interchange if necessary to ensure that the
     103             :            first transformation is othogonal, not hyperbolic.  */
     104         163 :         if (dd[jj+1] != dd[jj+2]) { /* make an interchange */
     105          45 :           if (dd[jj] != dd[jj+1]) {  /* interchange 1st and 2nd */
     106          20 :             buf[0] = bulge20; bulge20 = bulge10; bulge10 = buf[0];
     107          20 :             buf[0] = aa[jj]; aa[jj] = aa[jj+1]; aa[jj+1] = buf[0];
     108          20 :             buf[0] = bb[jj+1]; bb[jj+1] = bulge31; bulge31 = buf[0];
     109          20 :             buf[0] = dd[jj]; dd[jj] = dd[jj+1]; dd[jj+1] = buf[0];
     110         220 :             for (k=0;k<n;k++) {
     111         200 :               rtmp = uu[k+jj*ld]; uu[k+jj*ld] = uu[k+(jj+1)*ld]; uu[k+(jj+1)*ld] = rtmp;
     112             :             }
     113             :           } else {  /* interchange 1st and 3rd */
     114          25 :             buf[0] = bulge30; bulge30 = bulge10; bulge10 = buf[0];
     115          25 :             buf[0] = aa[jj]; aa[jj] = aa[jj+2]; aa[jj+2] = buf[0];
     116          25 :             buf[0] = bb[jj]; bb[jj] = bb[jj+1]; bb[jj+1] = buf[0];
     117          25 :             buf[0] = dd[jj]; dd[jj] = dd[jj+2]; dd[jj+2] = buf[0];
     118          25 :             if (jj + 2 < nn-1) {
     119          21 :               bulge41 = bb[jj+2];
     120          21 :               bb[jj+2] = 0;
     121             :             }
     122         275 :             for (k=0;k<n;k++) {
     123         250 :               rtmp = uu[k+jj*ld]; uu[k+jj*ld] = uu[k+(jj+2)*ld]; uu[k+(jj+2)*ld] = rtmp;
     124             :             }
     125             :           }
     126             :         }
     127             : 
     128             :         /* Set up transforming matrix rot. */
     129         163 :         PetscCall(UnifiedRotation(bulge20,bulge30,1,rot,&rcond,&swap));
     130             : 
     131             :         /* Apply transforming matrix rot to T. */
     132         163 :         bulge20 = rot[0]*bulge20 + rot[2]*bulge30;
     133         163 :         buf[0] = rot[0]*bb[jj] + rot[2]*bulge31;
     134         163 :         buf[1] = rot[1]*bb[jj] + rot[3]*bulge31;
     135         163 :         bb[jj] = buf[0];
     136         163 :         bulge31 = buf[1];
     137         163 :         buf[0] = rot[0]*rot[0]*aa[jj+1] + 2.0*rot[0]*rot[2]*bb[jj+1] + rot[2]*rot[2]*aa[jj+2];
     138         163 :         buf[1] = rot[1]*rot[1]*aa[jj+1] + 2.0*rot[3]*rot[1]*bb[jj+1] + rot[3]*rot[3]*aa[jj+2];
     139         163 :         bb[jj+1] = rot[1]*rot[0]*aa[jj+1] + rot[3]*rot[2]*aa[jj+2] + (rot[3]*rot[0] + rot[1]*rot[2])*bb[jj+1];
     140         163 :         aa[jj+1] = buf[0];
     141         163 :         aa[jj+2] = buf[1];
     142         163 :         if (jj + 2 < nn-1) {
     143         126 :           bulge42 = bb[jj+2]*rot[2];
     144         126 :           bb[jj+2] = bb[jj+2]*rot[3];
     145             :         }
     146             : 
     147             :         /* Accumulate transforming matrix */
     148         163 :         PetscCallBLAS("BLASrot",BLASMIXEDrot_(&n_,uu+(jj+1)*ld,&one,uu+(jj+2)*ld,&one,&rot[0],&rot[2]));
     149             :       }
     150             : 
     151             :       /* Annihilate inner entry bulge20 */
     152         192 :       if (bulge20 != 0.0) {
     153             : 
     154             :         /* Begin by setting up transforming matrix rot */
     155         192 :         sygn = dd[jj]*dd[jj+1];
     156         192 :         PetscCall(UnifiedRotation(bulge10,bulge20,sygn,rot,&rcond,&swap));
     157         192 :         if (rcond<PETSC_MACHINE_EPSILON) {
     158           0 :           *flag = PETSC_TRUE;
     159           0 :           PetscFunctionReturn(PETSC_SUCCESS);
     160             :         }
     161         192 :         if (rcond < worstcond) worstcond = rcond;
     162             : 
     163             :         /* Apply transforming matrix rot to T */
     164         192 :         if (jj > ntop) bb[jj-1] = rot[0]*bulge10 + rot[2]*bulge20;
     165         192 :         buf[0] = rot[0]*rot[0]*aa[jj] + 2*rot[0]*rot[2]*bb[jj] + rot[2]*rot[2]*aa[jj+1];
     166         192 :         buf[1] = rot[1]*rot[1]*aa[jj] + 2*rot[3]*rot[1]*bb[jj] + rot[3]*rot[3]*aa[jj+1];
     167         192 :         bb[jj] = rot[1]*rot[0]*aa[jj] + rot[3]*rot[2]*aa[jj+1] + (rot[3]*rot[0] + rot[1]*rot[2])*bb[jj];
     168         192 :         aa[jj] = buf[0];
     169         192 :         aa[jj+1] = buf[1];
     170         192 :         if (jj + 1 < nn-1) {
     171             :           /* buf = [ bulge31 bb(jj+1) ] * rot' */
     172         163 :           buf[0] = rot[0]*bulge31 + rot[2]*bb[jj+1];
     173         163 :           buf[1] = rot[1]*bulge31 + rot[3]*bb[jj+1];
     174         163 :           bulge31 = buf[0];
     175         163 :           bb[jj+1] = buf[1];
     176             :         }
     177         192 :         if (jj + 2 < nn-1) {
     178             :           /* buf = [bulge41 bulge42] * rot' */
     179         126 :           buf[0] = rot[0]*bulge41 + rot[2]*bulge42;
     180         126 :           buf[1] = rot[1]*bulge41 + rot[3]*bulge42;
     181         126 :           bulge41 = buf[0];
     182         126 :           bulge42 = buf[1];
     183             :         }
     184             : 
     185             :         /* Apply transforming matrix rot to D */
     186         192 :         if (swap == 1) {
     187          59 :           buf[0] = dd[jj]; dd[jj] = dd[jj+1]; dd[jj+1] = buf[0];
     188             :         }
     189             : 
     190             :         /* Accumulate transforming matrix, uu(jj:jj+1,:) = rot*uu(jj:jj+1,:) */
     191         192 :         if (sygn==1) {
     192          76 :           PetscCallBLAS("BLASrot",BLASMIXEDrot_(&n_,uu+jj*ld,&one,uu+(jj+1)*ld,&one,&rot[0],&rot[2]));
     193             :         } else {
     194         116 :           if (PetscAbsReal(rot[0])>PetscAbsReal(rot[1])) { /* Type I */
     195          57 :             t = rot[1]/rot[0];
     196         627 :             for (ii=0;ii<n;ii++) {
     197         570 :               uu[jj*ld+ii] = rot[0]*uu[jj*ld+ii] + rot[1]*uu[(jj+1)*ld+ii];
     198         570 :               uu[(jj+1)*ld+ii] = t*uu[jj*ld+ii] + uu[(jj+1)*ld+ii]/rot[0];
     199             :             }
     200             :           } else { /* Type II */
     201          59 :             t = rot[0]/rot[1];
     202         649 :             for (ii=0;ii<n;ii++) {
     203         590 :               rtmp = uu[jj*ld+ii];
     204         590 :               uu[jj*ld+ii] = rot[0]*uu[jj*ld+ii] + rot[1]*uu[(jj+1)*ld+ii];
     205         590 :               uu[(jj+1)*ld+ii] = t*uu[jj*ld+ii] + rtmp/rot[1];
     206             :             }
     207             :           }
     208             :         }
     209             :       }
     210             :     }
     211             : 
     212             :     /* Adjust bulge for next step */
     213         200 :     bulge10 = bb[jj];
     214         200 :     bulge20 = bulge31;
     215         200 :     bulge30 = bulge41;
     216         200 :     bulge31 = bulge42;
     217         200 :     bulge41 = 0.0;
     218         200 :     bulge42 = 0.0;
     219             :   }
     220          37 :   PetscFunctionReturn(PETSC_SUCCESS);
     221             : }
     222             : 
     223           3 : static PetscErrorCode HZIteration(PetscBLASInt nn,PetscBLASInt cgd,PetscReal *aa,PetscReal *bb,PetscReal *dd,PetscScalar *uu,PetscBLASInt ld)
     224             : {
     225           3 :   PetscBLASInt   j2,one=1,its,nits,nstop,jj,ntop,nbot,ntry;
     226           3 :   PetscReal      htr,det,dis,dif,tn,kt,c,s,tr,dt;
     227           3 :   PetscBool      flag=PETSC_FALSE;
     228             : 
     229           3 :   PetscFunctionBegin;
     230           3 :   its = 0;
     231           3 :   nbot = nn-1;
     232           3 :   nits = 0;
     233           3 :   nstop = 40*(nn - cgd);
     234             : 
     235          57 :   while (nbot >= cgd && nits < nstop) {
     236             : 
     237             :     /* Check for zeros on the subdiagonal */
     238          54 :     jj = nbot - 1;
     239         265 :     while (jj>=cgd && PetscAbs(bb[jj])>PETSC_MACHINE_EPSILON*(PetscAbs(aa[jj])+PetscAbs(aa[jj+1]))) jj = jj-1;
     240          54 :     if (jj>=cgd) bb[jj]=0;
     241          54 :     ntop = jj + 1;  /* starting point for step */
     242          54 :     if (ntop == nbot) {  /* isolate single eigenvalue */
     243             :       nbot = ntop - 1;
     244             :       its = 0;
     245          48 :     } else if (ntop+1 == nbot) {  /* isolate pair of eigenvalues */
     246          11 :       htr = 0.5*(aa[ntop]*dd[ntop] + aa[nbot]*dd[nbot]);
     247          11 :       det = dd[ntop]*dd[nbot]*(aa[ntop]*aa[nbot]-bb[ntop]*bb[ntop]);
     248          11 :       dis = htr*htr - det;
     249          11 :       if (dis > 0) {  /* distinct real eigenvalues */
     250           8 :         if (dd[ntop] == dd[nbot]) {  /* separate the eigenvalues by a Jacobi rotator */
     251           6 :           dif = aa[ntop]-aa[nbot];
     252           6 :           if (2.0*PetscAbs(bb[ntop])<=dif) {
     253           3 :             tn = 2*bb[ntop]/dif;
     254           3 :             tn = tn/(1.0 + PetscSqrtReal(1.0+tn*tn));
     255             :           } else {
     256           3 :             kt = dif/(2.0*bb[ntop]);
     257           5 :             tn = PetscSign(kt)/(PetscAbsReal(kt)+PetscSqrtReal(1.0+kt*kt));
     258             :           }
     259           6 :           c = 1.0/PetscSqrtReal(1.0 + tn*tn);
     260           6 :           s = c*tn;
     261           6 :           aa[ntop] = aa[ntop] + tn*bb[ntop];
     262           6 :           aa[nbot] = aa[nbot] - tn*bb[ntop];
     263           6 :           bb[ntop] = 0;
     264           6 :           j2 = nn-cgd;
     265           6 :           PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,uu+ntop*ld+cgd,&one,uu+nbot*ld+cgd,&one,&c,&s));
     266             :         }
     267             :       }
     268             :       nbot = ntop - 1;
     269             :     } else {  /* Do an HZ iteration */
     270          37 :       its = its + 1;
     271          37 :       nits = nits + 1;
     272          37 :       tr = aa[nbot-1]*dd[nbot-1] + aa[nbot]*dd[nbot];
     273          37 :       dt = dd[nbot-1]*dd[nbot]*(aa[nbot-1]*aa[nbot]-bb[nbot-1]*bb[nbot-1]);
     274          37 :       for (ntry=1;ntry<=6;ntry++) {
     275          37 :         PetscCall(HZStep(ntop,nbot+1,tr,dt,aa,bb,dd,uu,nn,ld,&flag));
     276          37 :         if (!flag) break;
     277           0 :         PetscCheck(ntry<6,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Unable to complete hz step after six tries");
     278           0 :         tr = 0.9*tr; dt = 0.81*dt;
     279             :       }
     280             :     }
     281             :   }
     282           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     283             : }
     284             : 
     285           4 : PetscErrorCode DSSolve_GHIEP_HZ(DS ds,PetscScalar *wr,PetscScalar *wi)
     286             : {
     287           4 :   PetscInt          i,off;
     288           4 :   PetscBLASInt      n,l,n1,ld = 0;
     289           4 :   const PetscScalar *A,*B;
     290           4 :   PetscScalar       *Q;
     291           4 :   PetscReal         *d,*e,*s;
     292             : 
     293           4 :   PetscFunctionBegin;
     294             : #if !defined(PETSC_USE_COMPLEX)
     295             :   PetscAssertPointer(wi,3);
     296             : #endif
     297           4 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     298           4 :   PetscCall(PetscBLASIntCast(ds->n-ds->l,&n1));
     299           4 :   off = ds->l + ds->l*ld;
     300           4 :   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
     301           4 :   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
     302           4 :   e = d + ld;
     303             : #if defined(PETSC_USE_DEBUG)
     304             :   /* Check signature */
     305           4 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
     306          39 :   for (i=0;i<ds->n;i++) {
     307          35 :     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
     308          35 :     PetscCheck(de==1.0 || de==-1.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diagonal elements of the signature matrix must be 1 or -1");
     309             :   }
     310           4 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
     311             : #endif
     312             :   /* Quick return */
     313           4 :   if (n1 == 1) {
     314           1 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
     315           6 :     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
     316           1 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
     317           1 :     PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
     318           1 :     if (ds->compact) {
     319           1 :       wr[ds->l] = d[ds->l]/s[ds->l]; wi[ds->l] = 0.0;
     320             :     } else {
     321           0 :       PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
     322           0 :       PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
     323           0 :       d[ds->l] = PetscRealPart(A[off]);
     324           0 :       s[ds->l] = PetscRealPart(B[off]);
     325           0 :       PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
     326           0 :       PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
     327           0 :       wr[ds->l] = d[ds->l]/s[ds->l]; wi[ds->l] = 0.0;
     328             :     }
     329           1 :     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
     330           1 :     PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
     331           1 :     PetscFunctionReturn(PETSC_SUCCESS);
     332             :   }
     333             :   /* Reduce to pseudotriadiagonal form */
     334           3 :   PetscCall(DSIntermediate_GHIEP(ds));
     335           3 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
     336           3 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     337           3 :   PetscCall(PetscBLASIntCast(ds->l,&l));
     338           3 :   PetscCall(HZIteration(n,l,d,e,s,Q,ld));
     339           3 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
     340           3 :   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
     341           3 :   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
     342           3 :   if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE));
     343             :   /* Undo from diagonal the blocks with real eigenvalues*/
     344           3 :   PetscCall(DSGHIEPRealBlocks(ds));
     345             : 
     346             :   /* Recover eigenvalues from diagonal */
     347           3 :   PetscCall(DSGHIEPComplexEigs(ds,0,ds->n,wr,wi));
     348             : #if defined(PETSC_USE_COMPLEX)
     349           3 :   if (wi) {
     350          31 :     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
     351             :   }
     352             : #endif
     353           3 :   ds->t = ds->n;
     354           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     355             : }