eigenvalue for 3x3 matrix

#include <stdio.h> 

#include <math.h> 

void getEigenValue(double a[][3], double eigen[3]) { 

     double b,c,d,f,g,h; 

     double i,j,k,l,m,n, sign, eq1, eq2, eq3; 

     double p,q,r,s,t,u,x1,x2,x3; 

     b=a[0][0]+a[1][1]+a[2][2]; 

     b=-b; 

     c=a[0][0]*(a[1][1]+a[2][2])+a[1][1]*a[2][2]; 

     c=c-a[1][2]*a[2][1]-a[0][1]*a[1][0]-a[0][2]*a[2][0]; 

     d=a[0][0]*(a[1][1]*a[2][2]-a[1][2]*a[2][1]) 

      +a[2][0]*(a[0][1]*a[1][2]-a[1][1]*a[0][2]) 

      +a[1][0]*(a[0][2]*a[2][1]-a[2][2]*a[0][1]); 

     d=-d; 

     f=(3.0*c-b*b)/3.0; 

     g=(2.0*b*b*b-9.0*b*c+27.0*d)/27.0; 

     h=g*g/4.0+f*f*f/27.0; 

     printf("b=%g c=%g d=%g\n", b, c, d); 

     //printf("f=%g g=%g h=%g\n", f, g, h); 

     if(h>0) { 

            r=-g/2.0+sqrt(h); 

            if(r<0) { r=-r; sign=-1.0; } else sign=1.0; 

            s=sign*pow(r, 1.0/3.0); 

            t=-g/2.0-sqrt(h); 

            if(t<0) { t=-t; sign=-1.0; } else sign=1.0; 

            u=sign*pow(t, 1.0/3.0); 

            x1=(s+u)-b/3.0; 

            eigen[0]=x1; eigen[1]=0.0; eigen[2]=0.0; 

     } 

     else if(h==0) { 

            if(d<0) { d=-d; sign=-1.0; } sign=1.0; 

            x1=-1.0*sign*pow(d, 1.0/3.0); 

            eigen[0]=x1;  eigen[1]=0.0; eigen[2]=0.0; 

     } 

     else { 

            i=g*g/4.0-h; 

            i=sqrt(i); 

            j=pow(i, 1.0/3.0); 

            k=-g/(2.0*i); 

            k=acos(k); 

            l=-j; 

            m=cos(k/3.0); 

            n=sqrt(3.0)*sin(k/3.0); 

            p=-b/3.0; 

            x1=2*j*m+p; 

            x2=l*(m+n)+p; 

            x3=l*(m-n)+p; 

            eq1=x1*x1*x1+b*x1*x1+c*x1+d; 

            eq2=x2*x2*x2+b*x2*x2+c*x2+d; 

            eq3=x3*x3*x3+b*x3*x3+c*x3+d; 

            if(fabs(eq1)<0.00000000001) eq1=0.0; 

            if(fabs(eq2)<0.00000000001) eq2=0.0; 

            if(fabs(eq3)<0.00000000001) eq3=0.0; 

            //printf("i=%g j=%g k=%g l=%g m=%g\n", i, j, k, l, m); 

            printf("x1=%.14f \teq1=%g\n", x1, eq1); 

            printf("x2=%.14f \teq2=%g\n", x2, eq2); 

            printf("x3=%.14f \teq3=%g\n", x3, eq3); 

            eigen[0]=x1; eigen[1]=x2; eigen[2]=x3; 

     } 

} 

int main() { 

     double a[3][3]={{-2,2,3},{-1,1,3},{2,0,-1}}; 

     double eigen[3]; 

     getEigenValue(a, eigen); 

     printf("eigen values are %.14f %.14f %.14f", 

         eigen[0], eigen[1], eigen[2]); 

     return 0; 

}