SYMMEIGEN - Functions solving the eigenvalue problem for symmetric matrices eig_jacobi - find eigenvalues and corresponding eigenvectors via the jacobi algorithm for symmetric matrices sort_eigenvalues - sort eigenvalues and corresponding eigenvectors in descending order Function Prototypes: void eig_jacobi(float **a, float d[], float **v, int n); void sort_eigenvalues(float d[], float **v, int n); (inspired by Press et. al., 1996) Macro used internally define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\ a[k][l]=h+s*(g-h*tau); void eig_jacobi(float **a, float d[], float **v, int n) eig_jacobi - find eigenvalues and corresponding eigenvectors via the jacobi algorithm for symmetric matrices Function Prototype: void eig_jacobi(float **a, float d[], float **v, int n); (inspired by Press et. al., 1996) { int j,iq,ip,i; float tresh,theta,tau,t,sm,s,h,g,c,*b,*z; /* allocate space temporarily b=alloc1float(n); b-=1; z=alloc1float(n); z-=1; /* initialize v to the identity matrix for (ip=1;ip<=n;ip++) { for (iq=1;iq<=n;iq++) v[ip][iq]=0.0; v[ip][ip]=1.0; } /* initialilize to the diagonal on matrix a for (ip=1;ip<=n;ip++) { b[ip]=d[ip]=a[ip][ip]; z[ip]=0.0; } /* main iteration loop for (i=1;i<=50;i++) { sm=0.0; for (ip=1;ip<=n-1;ip++) { for (iq=ip+1;iq<=n;iq++) sm += fabs(a[ip][iq]); } /* normal return if (sm == 0.0) { z+=1; free1float(z); b+=1; free1float(b); return; } /* tresh values for first 3 sweeps and therafter if (i < 4) tresh=0.2*sm/(n*n); else tresh=0.0; for (ip=1;ip<=n-1;ip++) { for (iq=ip+1;iq<=n;iq++) { g=100.0*fabs(a[ip][iq]); if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) a[ip][iq]=0.0; else if (fabs(a[ip][iq]) > tresh) { h=d[iq]-d[ip]; if ((float)(fabs(h)+g) == (float)fabs(h)) t=(a[ip][iq])/h; else { theta=0.5*h/(a[ip][iq]); t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); if (theta < 0.0) t = -t; } c=1.0/sqrt(1+t*t); s=t*c; tau=s/(1.0+c); h=t*a[ip][iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[ip][iq]=0.0; /* Jacobi rotations for (j=1;j<=ip-1;j++) { ROTATE(a,j,ip,j,iq) } for (j=ip+1;j<=iq-1;j++) { ROTATE(a,ip,j,j,iq) } for (j=iq+1;j<=n;j++) { ROTATE(a,ip,j,iq,j) } for (j=1;j<=n;j++) { ROTATE(v,j,ip,j,iq) } } } } for (ip=1;ip<=n;ip++) { b[ip] += z[ip]; d[ip]=b[ip]; z[ip]=0.0; } } /* this will not happen, hopefully fprintf(stderr,"jacobi iteration does not converge\n"); } void sort_eigenvalues(float d[], float **v, int n) sort_eigenvalues - sort eigenvalues and corresponding eigenvectors in descending order Function Prototypes: void sort_eigenvalues(float d[], float **v, int n); (inspired by Press et. al., 1996) { int k,j,i; float p; for (i=1;i= p) p=d[k=j]; if (k != i) { d[k]=d[i]; d[i]=p; for (j=1;j<=n;j++) { p=v[j][i]; v[j][i]=v[j][k]; v[j][k]=p; } } } }