diff --git a/code/math.f90 b/code/math.f90
index a3ef37ce5..43ea66e6f 100644
--- a/code/math.f90
+++ b/code/math.f90
@@ -150,8 +150,8 @@ module math
    math_sampleFiberOri, &
    math_sampleGaussVar, &
    math_symmetricEulers, &
-   math_spectralDecompositionSym33, &
-   math_spectralDecompositionSym, &
+   math_eigenvectorBasisSym33, &
+   math_eigenvectorBasisSym, &
    math_eigenValuesVectorsSym33, &
    math_eigenValuesVectorsSym, &
    math_rotationalPart33, &
@@ -1955,84 +1955,112 @@ end subroutine math_eigenValuesVectorsSym
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m
+!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m using an analytical expression
+!> and the general LAPACK powered version for arbritrary sized matrices as fallback
+!> @author Joachim Kopp, Max–Planck–Institut für Kernphysik, Heidelberg (Copyright (C) 2006)
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
 !--------------------------------------------------------------------------------------------------
 subroutine math_eigenValuesVectorsSym33(m,values,vectors)
-
+ 
  implicit none
- real(pReal), dimension(3,3), intent(in)  :: m
- real(pReal), dimension(3),   intent(out) :: values
- real(pReal), dimension(3,3), intent(out) :: vectors
+ real(pReal), dimension(3,3),intent(in)  :: m
+ real(pReal), dimension(3),  intent(out) :: values
+ real(pReal), dimension(3,3),intent(out) :: vectors
+ real(pReal) :: T, U, norm, threshold
  logical :: error
- integer(pInt) :: info
- real(pReal), dimension((64+2)*3) :: work                                                            ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
 
- vectors = m                                                                                         ! copy matrix to input (doubles as output) array
-#if(FLOAT==8)
- call dsyev('V','U',3,vectors,3,values,work,(64+2)*3,info)
-#elif(FLOAT==4)
- call ssyev('V','U',3,vectors,3,values,work,(64+2)*3,info)
-#endif
- error = (info == 0_pInt)
+ values = math_eigenvaluesSym33(m)
 
+ vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), &
+                    m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), &
+                    m(1, 2)**2_pInt]
+
+ T = maxval(abs(values))
+ U = max(T, T**2_pInt)
+ threshold = sqrt(5.68e-14_pReal * U**2_pInt)
+
+! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2
+ vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), &
+                    vectors(2,2) + m(2, 3) * values(1), &
+                    (m(1,1) - values(1)) * (m(2,2) - values(1)) - vectors(3,2)]
+ norm = norm2(vectors(1:3, 1))
+
+ fallback1: if(norm < threshold) then
+   call math_eigenValuesVectorsSym(m,values,vectors,error)
+   return
+ endif fallback1
+
+ vectors(1:3,1) = vectors(1:3, 1) / norm 
+ 
+! Calculate second eigenvector by the formula v[1] = (m - lambda[1]).e1 x (m - lambda[1]).e2
+ vectors(1:3,2) = [ vectors(1,2) + m(1, 3) * values(2), &
+                    vectors(2,2) + m(2, 3) * values(2), &
+                    (m(1,1) - values(2)) * (m(2,2) - values(2)) - vectors(3,2)]
+ norm = norm2(vectors(1:3, 2))
+
+ fallback2: if(norm < threshold) then
+   call math_eigenValuesVectorsSym(m,values,vectors,error)
+   return
+ endif fallback2
+ vectors(1:3,2) = vectors(1:3, 2) / norm
+
+! Calculate third eigenvector according to  v[2] = v[0] x v[1]
+ vectors(1:3,3) = math_crossproduct(vectors(1:3,1),vectors(1:3,2))
 end subroutine math_eigenValuesVectorsSym33
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief eigenvalues and eigenvectors of symmetric matrix m
+!> @brief eigenvector basis of symmetric matrix m
 !--------------------------------------------------------------------------------------------------
-function math_spectralDecompositionSym(m)
+function math_eigenvectorBasisSym(m)
 
  implicit none
  real(pReal), dimension(:,:),      intent(in)  :: m
  real(pReal), dimension(size(m,1))             :: values
  real(pReal), dimension(size(m,1),size(m,1))   :: vectors
- real(pReal), dimension(size(m,1),size(m,1))   :: math_spectralDecompositionSym
+ real(pReal), dimension(size(m,1),size(m,1))   :: math_eigenvectorBasisSym
  logical :: error
  integer(pInt) :: i
 
- math_spectralDecompositionSym = 0.0_pReal
+ math_eigenvectorBasisSym = 0.0_pReal
  call math_eigenValuesVectorsSym(m,values,vectors,error)
  if(error) return
  
  do i=1_pInt, size(m,1)
-   math_spectralDecompositionSym = math_spectralDecompositionSym &
-                                 +  sqrt(values(i)) * math_tensorproduct(vectors(:,i),vectors(:,i))
+   math_eigenvectorBasisSym = math_eigenvectorBasisSym &
+                            +  sqrt(values(i)) * math_tensorproduct(vectors(:,i),vectors(:,i))
  enddo
 
-end function math_spectralDecompositionSym
+end function math_eigenvectorBasisSym
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m
+!> @brief eigenvector basis of symmetric 33 matrix m
 !--------------------------------------------------------------------------------------------------
-function math_spectralDecompositionSym33(m)
+function math_eigenvectorBasisSym33(m)
 
  implicit none
- real(pReal), dimension(3,3)              :: math_spectralDecompositionSym33
- real(pReal), dimension(3)                :: invariants, values,C
+ real(pReal), dimension(3,3)              :: math_eigenvectorBasisSym33
+ real(pReal), dimension(3)                :: invariants, values
  real(pReal), dimension(3,3), intent(in) :: m
  real(pReal) :: P, Q, rho, phi
  real(pReal), parameter :: TOL=1.e-14_pReal
- real(pReal), dimension(3,3) :: M1, M2, M3,EB1, EB2, EB3
- real(pReal) :: D1, D2, D3
+ real(pReal), dimension(3,3,3) :: N, EB
 
  invariants = math_invariantsSym33(m)
+ EB = 0.0_pReal
 
- P=invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
- Q=-2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+invariants(1)*invariants(2)/3.0_pReal-invariants(3)
+ P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
+ Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
 
- EB1=0.0_pReal
- EB2=0.0_pReal
- EB3=0.0_pReal
- if((ABS(P) < TOL).AND.(ABS(Q) < TOL)) then ! EV_2 = EV_1 = EV_3
+ threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then
    values = invariants(1)/3.0_pReal
-!   this is not really correct, but this way U is calculated
-!   correctly in PDECOMPOSITION (correct is EB?=I)
-   EB1(1,1)=1.0_pReal
-   EB2(2,2)=1.0_pReal
-   EB3(3,3)=1.0_pReal
- else
+!   this is not really correct, but at least the basis is correct
+   EB(1,1,1)=1.0_pReal
+   EB(2,2,2)=1.0_pReal
+   EB(3,3,3)=1.0_pReal
+ else threeSimilarEigenvalues
    rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
    phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
    values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
@@ -2040,35 +2068,36 @@ function math_spectralDecompositionSym33(m)
                               cos((phi+2.0_pReal*PI)/3.0_pReal), &
                               cos((phi+4.0_pReal*PI)/3.0_pReal) &
                              ] + invariants(1)/3.0_pReal
-   C = abs([values(1)-values(2),values(2)-values(3),values(3)-values(1)])
-   M1=M-values(1)*math_I3
-   M2=M-values(2)*math_I3
-   M3=M-values(3)*math_I3
-   if (C(1) < TOL) then ! EV_2 = EV_1, no contribution from EV_2
-     D3=1.0_pReal/(values(3)-values(1))/(values(3)-values(2))
-     EB3=math_mul33x33(M1,M2)*D3
-     EB1=math_I3-EB3
-   elseif (C(2) < TOL) then ! EV_2 = EV_3, no contribution from EV_3
-     D1=1.0_pReal/(values(1)-values(2))/(values(1)-values(3))
-     EB1=math_mul33x33(M2,M3)*D1
-     EB2=math_I3-EB1
-   elseif(C(3) < TOL) then ! EV_1 = EV_3, no contribution from EV_3
-     D2=1.0_pReal/(values(2)-values(1))/(values(2)-values(3))
-     EB2=math_mul33x33(M1,M3)*D2
-     EB1=math_I3-EB2
-   else !  all three eigenvectors are different
-     D1=1.0_pReal/(values(1)-values(2))/(values(1)-values(3))
-     D2=1.0_pReal/(values(2)-values(1))/(values(2)-values(3))
-     D3=1.0_pReal/(values(3)-values(1))/(values(3)-values(2))
-     EB1=math_mul33x33(M2,M3)*D1
-     EB2=math_mul33x33(M1,M3)*D2
-     EB3=math_mul33x33(M1,M2)*D3
-   endif
- endif
+   N(1:3,1:3,1) = m-values(1)*math_I3
+   N(1:3,1:3,2) = m-values(2)*math_I3
+   N(1:3,1:3,3) = m-values(3)*math_I3
+   twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then 
+     EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
+                                               ((values(3)-values(1))*(values(3)-values(2)))
+     EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
+   elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
+     EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
+                                               ((values(1)-values(2))*(values(1)-values(3)))
+     EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
+   elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues 
+     EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
+                                               ((values(2)-values(1))*(values(2)-values(3)))
+     EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
+   else twoSimilarEigenvalues
+     EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
+                                               ((values(1)-values(2))*(values(1)-values(3)))
+     EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
+                                               ((values(2)-values(1))*(values(2)-values(3)))
+     EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
+                                               ((values(3)-values(1))*(values(3)-values(2)))
+   endif twoSimilarEigenvalues
+ endif threeSimilarEigenvalues
 
- math_spectralDecompositionSym33 = sqrt(values(1)) * EB1 + sqrt(values(2)) * EB2 + sqrt(values(3)) * EB3
+ math_eigenvectorBasisSym33 = sqrt(values(1)) * EB(1:3,1:3,1) &
+                            + sqrt(values(2)) * EB(1:3,1:3,2) &
+                            + sqrt(values(3)) * EB(1:3,1:3,3)
 
-end function math_spectralDecompositionSym33
+end function math_eigenvectorBasisSym33
 
 
 !--------------------------------------------------------------------------------------------------
@@ -2083,7 +2112,7 @@ function math_rotationalPart33(m)
  real(pReal), dimension(3,3) :: math_rotationalPart33
  real(pReal), dimension(3,3) :: U , Uinv
 
- U = math_spectralDecompositionSym33(math_mul33x33(transpose(m),m))
+ U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m))
  Uinv = math_inv33(U)
 
  if (all(abs(Uinv) <= tiny(Uinv))) then                                                             ! math_inv33 returns zero when failed, avoid floating point equality comparison
@@ -2143,7 +2172,7 @@ function math_eigenvaluesSym33(m)
  P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
  Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
 
- if(any(abs([p,q]) < TOL)) then
+ if(all(abs([P,Q]) < TOL)) then
    math_eigenvaluesSym33 = math_eigenvaluesSym(m)
  else
    rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
diff --git a/code/plastic_disloUCLA.f90 b/code/plastic_disloUCLA.f90
index d95a5e6a4..c9e3bc846 100644
--- a/code/plastic_disloUCLA.f90
+++ b/code/plastic_disloUCLA.f90
@@ -1207,13 +1207,11 @@ subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
    math_Plain3333to99, &
    math_Mandel6to33, &
    math_Mandel33to6, &
-   math_spectralDecompositionSym33, &
    math_symmetric33, &
    math_mul33x3
  use material, only: &
    material_phase, &
    phase_plasticityInstance, &
-   !plasticState, &
    phaseAt, phasememberAt
  use lattice, only: &
    lattice_Sslip, &