jacobi.cpp

VC++常用数据算法集,里面包含了大量的常用算法程序,很有用的哟!

void jacobi(double a[4][4], int n, double d[], double v[4][4], int& nrot)
{
    double b[100], z[100],g,s,c,h,tau,sss,ddd,t,tresh,sm,theta;
	int ip,iq,i,j;
    for (ip = 1; ip<=n; ip++)
	{
        for (iq = 1; iq<=n; iq++)
		{
            v[ip][iq] = 0.0;
        }
        v[ip][ip] = 1.0;
    }
    for (ip = 1; ip<=n; ip++)
	{
        b[ip] = a[ip][ip];
        d[ip] = b[ip];
        z[ip] = 0.0;
    }
    nrot = 0;
    for (i = 1; i<=50; i++)
	{
        sm = 0.0;
        for (ip = 1; ip<=n - 1; ip++)
		{
            for (iq = ip + 1; iq<=n; iq++)
			{
                sm = sm + fabs(a[ip][iq]);
            }
        }
        if (sm == 0.0)
		{
			return;
		}
        if (i < 4)
		{
            tresh = 0.2 * sm /double(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]);
                sss = fabs(d[ip]) + g;
                ddd = fabs(d[iq]) + g;
                if ((i > 4) && (sss == fabs(d[ip])) && (ddd == fabs(d[iq])))
				{
                    a[ip][iq] = 0.0;
				}
                else 
				{
					if (fabs(a[ip][iq]) > tresh)
					{
						h = d[iq] - d[ip];
						if ((fabs(h) + g) == fabs(h))
						{
							t = a[ip][iq] / h;
						}
                        else
						{
						  theta=0.5 * h / a[ip][iq];
						  t=1.0/(fabs(theta)+sqrt(1.0+pow(theta,2)));
						  if (theta < 0.0 ) 
						  t = -t;
						}
						c = 1.0 / sqrt(1.0 + t *t);
						s = t * c;
						tau = s / (1.0 + c);
						h = t * a[ip][iq];
						z[ip] = z[ip] - h;
						z[iq] = z[iq] + h;
						d[ip] = d[ip] - h;
						d[iq] = d[iq] + h;
						a[ip][iq] = 0.0;
						for (j = 1; j<=ip - 1; j++)
						{
							g = a[j][ip];
							h = a[j][iq];
							a[j][ip] = g - s * (h + g * tau);
							a[j][iq] = h + s * (g - h * tau);
						}
						for (j = ip + 1 ; j<=iq - 1; j++)
						{
							g = a[ip][j];
							h = a[j][iq];
							a[ip][j] = g - s * (h + g * tau);
							a[j][iq] = h + s * (g - h * tau);
						}
						for (j = iq + 1 ;j<=n; j++)
						{
							g = a[ip][j];
							h = a[iq][j];
							a[ip][j] = g - s * (h + g * tau);
							a[iq][j] = h + s * (g - h * tau);
						}
						for (j = 1 ;j<=n; j++)
						{
							g = v[j][ip];
							h = v[j][iq];
							v[j][ip] = g - s * (h + g * tau);
							v[j][iq] = h + s * (g - h * tau);
						}
						nrot = nrot + 1;
					}
				}
			}
        }
        for (ip = 1; ip<=n; ip++)
		{
            b[ip] = b[ip] + z[ip];
            d[ip] = b[ip];
            z[ip] = 0.0;
        }
    }
    cout<< "50 iterations should never happen"<<endl;
}