d9r4f.java

解具有导数的多元函数的无约束极小值点或具有导数的非线性方程组的近似解

                    v[(j - 1) * n + i] = 0.0;
                }
            }
            v[(i - 1) * n + i] = 1.0;
            g = rv1[i];
            l = i;
        }
        for (i = n; i >= 1; i--)
        {
            l = i + 1;
            g = w[i];
            if (i < n) 
            {
                for (j = l; j <= n; j++)
                {
                    a[(i - 1) * n + j] = 0.0;
                }
            }
            if (g != 0.0) 
            {
                g = 1.0 / g;
                if (i != n)
                {
                    for (j = l; j <= n; j++)
                    {
                        s = 0.0;
                        for (k = l; k <= m; k++)
                        {
                            s = s + a[(k - 1) * n + i] * a[(k - 1) * n + j];
                        }
                        f = (s / a[(i - 1) * n + i]) * g;
                        for (k = i; k <= m; k++)
                        {    
                            a[(k - 1) * n + j] = a[(k - 1) * n + j] + f * a[(k - 1) * n + i];
                        }
                    }
                }
                for (j = i; j <= m; j++)
                {
                    a[(j - 1) * n + i] = a[(j - 1) * n + i] * g;
                }
            }
            else
            {
                for(j = i; j <= m; j++)
                {
                    a[(j - 1) * n + i] = 0.0;
                }
            }
            a[(i - 1) * n + i] = a[(i - 1) * n + i] + 1.0;
        }
        for(k = n; k >= 1; k--)
        {
            for (int its = 1; its <= 30; its++)
            {
                for (l = k; l >= 1; l--)
                {
                    nm = l - 1;
                    if ((Math.abs(rv1[l]) + anorm) == anorm)
                    {
                        break;
                    }
                    if ((Math.abs(w[nm]) + anorm) == anorm)
                    {
                        break;
                    }
                }
                if (Math.abs(rv1[l]) + anorm != anorm)
                {
                    c = 0.0;
                    s = 1.0;
                    for(i = l; i <= k; i++)
                    {
                        f = s * rv1[i];
                        if ((Math.abs(f) + anorm) != anorm)
                        {
                            g = w[i];
                            h = Math.sqrt(f * f + g * g);
                            w[i] = h;
                            h = 1.0 / h;
                            c = (g * h);
                            s = -(f * h);
                            for (j = 1;j <= m; j++)
                            {
                                y = a[(j - 1) * n + nm];
                                z = a[(j - 1) * n + i];
                                a[(j - 1) * n + nm] = (y * c) + (z * s);
                                a[(j - 1) * n + i] = -(y * s) + (z * c);
                            }
                        }
                    }
                }
                z = w[k];
                if (l == k)
                {
                    if (z < 0.0)
                    {
                        w[k] = -z;
                        for (j = 1; j <= n; j++)
                        {
                            v[(j - 1) * n + k] = -v[(j - 1) * n + k];
                        }
                    }
                    break; 
                }
                if (its == 30)
                {
                    System.out.println("no convergence in 30 iterations");
                    System.exit(1);                
                }
                x = w[l];
                nm = k - 1;
                y = w[nm];
                g = rv1[nm];
                h = rv1[k];
                f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
                g = Math.sqrt(f * f + 1.0);
                f = ((x - z) * (x + z) + h * ((y / (f + Math.abs(g) * sgn(f))) - h)) / x;
                c = 1.0;
                s = 1.0;
                for (j = l; j <= nm; j++)
                {
                    i = j + 1;
                    g = rv1[i];
                    y = w[i];
                    h = s * g;
                    g = g * c;
                    z = Math.sqrt(f * f + h * h);
                    rv1[j] = z;
                    c = f / z;
                    s = h / z;
                    f = (x * c) + (g * s);
                    g = -(x * s) + (g * c);
                    h = y * s;
                    y = y * c;
                    for (int pp = 1; pp <= n; pp++)
                    {
                        x = v[(pp - 1) * n + j];
                        z = v[(pp - 1) * n + i];
                        v[(pp - 1) * n + j] = (x * c) + (z * s);
                        v[(pp - 1) * n + i] = -(x * s) + (z * c);
                    }
                    z = Math.sqrt(f * f + h * h);
                    w[j] = z;
                    if (z != 0.0)
                    {
                        z = 1.0 / z;
                        c = f * z;
                        s = h * z;
                    }
                    f = (c * g) + (s * y);
                    x = -(s * g) + (c * y);
                    for (int qq = 1; qq <= m; qq++)
                    {
                        y = a[(qq - 1) * n + j];
                        z = a[(qq - 1) * n + i];
                        a[(qq - 1) * n + j] = (y * c) + (z * s);
                        a[(qq - 1) * n + i] = -(y * s) + (z * c);
                    }
                }
                rv1[l] = 0.0;
                rv1[k] = f;
                w[k] = x;
            } 
        }
    }

    static int iset;
    static double gset;
    int gasdev_idum;
    double gasdev(int idum)
    {     
        double t, v1, v2, fac, r;
        if (iset == 0.0)
        {
            do
            {
                v1 = 2.0 * ran1(idum) - 1.0;
                idum=ran1_idum;
                v2 = 2.0 * ran1(idum) - 1.0;
                idum=ran1_idum;
                r = v1 * v1 + v2 * v2;
            }while ((r >= 1.0) || (r == 0.0));
            fac = Math.sqrt(-2.0 * Math.log(r) / r);
            gset = v1 * fac;
            t = v2 * fac;
            iset = 1;
            gasdev_idum = idum;
            return t;
        }
        else
        {
            t = gset;
            gasdev_idum = idum;
            iset = 0;
            return t;
        }
    }

	void fleg(double x, double pl[], int nl)
	{
		pl[1] = 1.0;
		pl[2] = x;
		if (nl > 2)
		{
			double twox = 2.0 * x;
			double f2 = x;
			double d = 1.0;
			for (int j = 3; j <= nl; j++)
			{
				double f1 = d;
				f2 = f2 + twox;
				d = d + 1.0;
				pl[j] = (f2 * pl[j - 1] - f1 * pl[j - 2]) / d;
			}
		}
	}

	int ran1_idum;
	static double r[] = new double[98];
    static int ix1, ix2, ix3;
	static int iff;
	double ran1(int idum)
	{
        double rm1, rm2 ,t;
        int m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, j;         
        m1 = 259200; ia1 = 7141; ic1 = 54773; rm1 = 0.0000038580247;
        m2 = 134456; ia2 = 8121; ic2 = 28411; rm2 = 0.0000074373773;
        m3 = 243000; ia3 = 4561; ic3 = 51349;
        if ((idum < 0) || (iff == 0))
		{
            iff = 1;
            ix1 =  (ic1 - idum) % m1;
            ix1 = (ia1 * ix1 + ic1) % m1;
            ix2 = ix1 % m2;
            ix1 = (ia1 * ix1 + ic1) % m1;
            ix3 = ix1 % m3;
            for (j = 1; j <= 97; j++)
			{
                ix1 = (ia1 * ix1 + ic1) % m1;
                ix2 = (ia2 * ix2 + ic2) % m2;
                r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
            }
            idum = 1;
        }
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = (ia2 * ix2 + ic2) % m2;
        ix3 = (ia3 * ix3 + ic3) % m3;
        j = 1 + (int)((97 * ix3) / m3);
        if ((j > 97) || (j < 1)) 
        {
            System.out.println( "abnormal exit");
            System.exit(1);
        }
        t = r[j];
        r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
        ran1_idum = idum;
        return t;
    }

	int sgn(double pa)
	{
		if (pa > 0.0)
		{
			return 1;
		}
		else
		{
			if (pa < 0.0)
			{
				return -1;
			}
		}
		return 0;
	}
}