anchor.diff:     ! note: SGGEV seems to have memory issues in certain cases. The eigenvalues seem to be okay, but the eigenvectors vary wildly with different compiling options.
      !       DGGEV seems to work better, so I'm making these variables LAKi (which is set to R8Ki for now)   - bjj 4/25/2014
      ! bjj: This comes from the LAPACK documentation:
  -   !   Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) may easily over- or underflow, and BETA(j) may even be zero.
  -   !   Thus, the user should avoid naively computing the ratio alpha/beta.  However, ALPHAR and ALPHAI will be always less
  +   !   Note: the quotients AlphaR(j)/BETA(j) and AlphaI(j)/BETA(j) may easily over- or underflow, and BETA(j) may even be zero.
  +   !   Thus, the user should avoid naively computing the ratio Alpha/beta.  However, AlphaR and AlphaI will be always less
      !   than and usually comparable with norm(A) in magnitude, and BETA always less than and usually comparable with norm(B).    
  -   ! Omega2=ALPHAR/BETA  !Note this may not be correct if ALPHAI<>0 and/or BETA=0 TO INCLUDE ERROR CHECK, also they need to be sorted
  -   CALL  LAPACK_ggev('N','V',N ,K, M, ALPHAR, ALPHAI, BETA, VL, EigVect, WORK, LWORK, ErrStat2, ErrMsg2)
  +   ! Omega2=AlphaR/BETA  !Note this may not be correct if AlphaI<>0 and/or BETA=0 TO INCLUDE ERROR CHECK, also they need to be sorted
  +   CALL  LAPACK_ggev('N','V',N ,K, M, AlphaR, AlphaI, Beta, VL, EigVect, WORK, LWORK, ErrStat2, ErrMsg2)
      if(Failed()) return
   
      ! --- Determinign and sorting eigen frequencies 
  +   Omega2(:) =0.0_LaKi
      DO I=1,N !Initialize the key and calculate Omega
         KEY(I)=I
  -      IF ( AlphaR(I)<0.0_LAKi ) THEN
  -         !print*,'I', I, AlphaR(I), AlphaI(I), Beta(I), real(AlphaR(I)/Beta(I))
  -         if (bCheckSingularity) then
  -            ErrStat2=ErrID_Fatal
  -            ErrMsg2= 'Negative eigenvalue found, system may be singular (may contain rigid body modes)'
  -            if(Failed()) return
  +      if ( EqualRealNos(real(Beta(I),ReKi),0.0_ReKi) ) then
  +         ! --- Beta =0 
  +         print*,'Large eigenvalue found, system may be ill-conditioned'
  +         Omega2(I) = HUGE(1.0_LaKi) 
  +      elseif ( EqualRealNos(real(AlphaI(I),ReKi),0.0_ReKi) ) THEN
  +         ! --- Real Eigenvalues
  +         IF ( AlphaR(I)<0.0_LaKi ) THEN
  +            if ( (AlphaR(I)/Beta(I))<1e-6_LaKi ) then
  +               ! Tolerating very small negative eigenvalues
  +               print*,'Negative eigenvalue found with small norm, system may be singular (may contain rigid body modes)'
  +               Omega2(I)=0.0_LaKi
  +            else
  +               print*,'Negative eigenvalue found, system may be singular (may contain rigid body modes)'
  +               Omega2(I)=AlphaR(I)/Beta(I)
  +            endif
  +         else
  +            Omega2(I) = AlphaR(I)/Beta(I)
            endif
  -      !ELSE IF ( .not.EqualRealNos(AlphaI(I),0.0_LAKi )) THEN
  -      !   print*,'I', I, AlphaR(I), AlphaI(I), Beta(I), real(AlphaR(I)/Beta(I))
  -      !   ErrStat2=ErrID_Fatal
  -      !   ErrMsg2= 'Complex eigenvalue found, system may be singular (may contain rigid body modes)'
  -      !   if(Failed()) return
  -      ELSE IF ( EqualRealNos(real(Beta(I),ReKi),0.0_ReKi) ) THEN
  -         Omega2(I) = HUGE(Omega2)  ! bjj: should this be an error?
  -         if (bCheckSingularity) then
  -            ErrStat2=ErrID_Fatal
  -            ErrMsg2= 'Large eigenvalue found, system may be singular (may contain rigid body modes)'
  -            if(Failed()) return
  +      else
  +         ! --- Complex Eigenvalues
  +         normA = sqrt(AlphaR(I)**2 + AlphaI(I)**2)
  +         if ( (normA/Beta(I))<1e-6_LaKi ) then
  +            ! Tolerating very small eigenvalues with imaginary part
  +            print*,'Complex eigenvalue found with small norm, approximating as 0'
  +            Omega2(I) = 0.0_LaKi
  +         elseif ( abs(AlphaR(I))>1e3_LaKi*abs(AlphaI(I)) ) then
  +            ! Tolerating very small imaginary part compared to real part... (not pretty)
  +            print*,'Complex eigenvalue found with small Im compare to Re'
  +            !print*'AlphaR,AlphaI,beta=',AlphaR(I), AlphaI(I), beta(I), AlphaR(I)/beta(I), AlphaI(I)/beta(I)
  +            Omega2(I) = AlphaR(I)/Beta(I)
  +         else
  +            print*,'Complex eigenvalue found, AlphaR,AlphaI,beta=',AlphaR(I), AlphaI(I), beta(I), AlphaR(I)/beta(I), AlphaI(I)/beta(I)
  +            Omega2(I) = HUGE(1.0_LaKi)
            endif
  -      ELSE IF ( EqualRealNos(real(ALPHAR(I),ReKi),0.0_ReKi) ) THEN
  -         if (bCheckSingularity) then
  -            ErrStat2=ErrID_Fatal
  -            ErrMsg2= 'Eigenvalue zero found, system is singular (may contain rigid body modes)'
  -            if(Failed()) return
  -         endif
  -      ELSE
  -         Omega2(I) = REAL( ALPHAR(I)/BETA(I), ReKi )
  -      END IF           
  -   ENDDO  
  +      endif
  +   enddo  
      ! Sorting
      CALL ScaLAPACK_LASRT('I',N,Omega2,KEY,ErrStat2,ErrMsg2); if(Failed()) return 
       
  
positive.code:    LWORK=8*N + 16  !this is what the eigensolver wants  >> bjj: +16 because of MKL ?ggev documenation ( "lwork >= max(1, 8n+16) for real flavors"), though LAPACK documenation says 8n is fine
     !bjj: there seems to be a memory problem in *GGEV, so I'm making the WORK array larger to see if I can figure it out
     CALL AllocAry( Work,    LWORK, 'Work',   ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') 
     CALL AllocAry( ALPHAR,  N,     'ALPHAR', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve')
     CALL AllocAry( ALPHAI,  N,     'ALPHAI', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve')
     CALL AllocAry( BETA,    N,     'BETA',   ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve')
     CALL AllocAry( VL,      N,  N, 'VL',     ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve')
     CALL AllocAry( KEY,     N,     'KEY',    ErrStat2, ErrMsg2 ); if(Failed()) return
      
  
negative.code: 
     SUBROUTINE CleanupEigen()
        IF (ALLOCATED(Work)  ) DEALLOCATE(Work)
        IF (ALLOCATED(AlphaR)) DEALLOCATE(AlphaR)
        IF (ALLOCATED(AlphaI)) DEALLOCATE(AlphaI)
        IF (ALLOCATED(Beta)  ) DEALLOCATE(Beta)
        IF (ALLOCATED(VL)    ) DEALLOCATE(VL)
        IF (ALLOCATED(KEY)   ) DEALLOCATE(KEY)
     END SUBROUTINE CleanupEigen
  