#include "math.h"
void jcbj(a,n,v,eps)
int n;
double a[],v[],eps;
{ 
	int i,j,p,q,u,w,t,s;
    double ff,fm,cn,sn,omega,x,y,d;
    for (i=0; i<=n-1; i++)
    { 
		v[i*n+i]=1.0;
        for (j=0; j<=n-1; j++)
			if (i!=j) v[i*n+j]=0.0;
      }
    ff=0.0;
    for (i=1; i<=n-1; i++)
    for (j=0; j<=i-1; j++)
    { 
		d=a[i*n+j]; 
		ff=ff+d*d; 
	}
    ff=sqrt(2.0*ff);
    loop0:
    ff=ff/(1.0*n);


    loop1:
		for (i=1; i<=n-1; i++)
        for (j=0; j<=i-1; j++)
        { 
			d=fabs(a[i*n+j]);
            if (d>ff)
			{	p=i; q=j;
                goto loop;
            }
        }
        if (ff<eps) return;
        goto loop0;


	loop: 
		u=p*n+q; w=p*n+p; t=q*n+p; s=q*n+q;
        x=-a[u]; y=(a[s]-a[w])/2.0;
        omega=x/sqrt(x*x+y*y);
        if (y<0.0) omega=-omega;
        sn=1.0+sqrt(1.0-omega*omega);
        sn=omega/sqrt(2.0*sn);
        cn=sqrt(1.0-sn*sn);
        fm=a[w];
        a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega;
        a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega;
        a[u]=0.0; a[t]=0.0;
        for (j=0; j<=n-1; j++)
		{
			if ((j!=p)&&(j!=q))
			{ 
				u=p*n+j; w=q*n+j;
				fm=a[u];
				a[u]=fm*cn+a[w]*sn;
			    a[w]=-fm*sn+a[w]*cn;
			}
		}
        for (i=0; i<=n-1; i++)
		{
			if ((i!=p)&&(i!=q))
            { 
				u=i*n+p; w=i*n+q;
				fm=a[u];
				a[u]=fm*cn+a[w]*sn;
				a[w]=-fm*sn+a[w]*cn;
            }
		}
        for (i=0; i<=n-1; i++)
        { 
			u=i*n+p; w=i*n+q;
            fm=v[u];
            v[u]=fm*cn+v[w]*sn;
            v[w]=-fm*sn+v[w]*cn;
        }
       goto loop1;
}