要是实在没有java的源码,C的也成。

解决方案 »

  1.   

    public class d1r1F
    {
    void gaussj(double a[], int n, double b[])
    {
            int i, j, k, l, ll, irow, icol;
            icol = 0;
            irow = 0;
            double big, pivinv, dum;
            int ipiv[] = new int[50];
            int indxr[] = new int[50];
            int indxc[] = new int[50];
            for (j = 0; j <= n - 1; j++)
            {            
                ipiv[j] = 0;
            }
            for (i = 0; i <= n - 1; i++)
    {
                big = 0.0;
                for (j = 0; j <= n - 1; j++)
    {
                    if (ipiv[j] != 1)
    {
                        for (k = 0; k <= n - 1; k++)
    {
                            if (ipiv[k] == 0)
    {
                                if (Math.abs(a[j * n + k]) >= big)
    {
                                    big = Math.abs(a[j * n + k]);
                                    irow = j;
                                    icol = k;
    }
                                else if (ipiv[k] > 1)
    {
    System.out.println("singular matrix");
    }
    }
                        }
                    }
                }
                ipiv[icol] = ipiv[icol] + 1;
                if (irow != icol)
    {
                    for (l = 0; l <= n - 1; l++)
    {
                        dum = (a[irow * n + l]);
                        a[irow * n + l] = a[icol * n + l];
                        a[icol * n + l] = dum;
                    }
                    dum = b[irow];
                    b[irow] = b[icol];
                    b[icol] = dum;
    }
                indxr[i] = irow;
                indxc[i] = icol;
                if (a[icol * n + icol] == 0.0)
    {
    System.out.println( "singular matrix.");
                    System.exit(1);
    }
                pivinv = 1.0 / (a[icol * n + icol]);
                a[icol * n + icol] = 1.0;
                for (l = 0; l <= n - 1; l++)
                {
                    a[icol * n + l] = a[icol * n + l] * pivinv;
                }
                b[icol]=b[icol] * pivinv;
                for (ll = 0; ll <= n - 1; ll++)
                {
                    if (ll != icol)
                    {
                        dum = a[ll * n + icol];
                        a[ll * n + icol] = 0.0;
                        for (l = 0; l <= n - 1; l++)
                        {
                            a[ll * n + l] = a[ll * n + l] - a[icol * n + l] * dum;
                        }
                        b[ll] = b[ll] - b[icol] * dum;
                    }
                } 
                for (l = n - 1; l <= 0; l--)
                {
                    if (indxr[l] != indxc[l])
                    {
                        for (k = 0; k <= n - 1; k++)
                        {
                            dum = a[k * n + indxr[l]];
                            a[k * n + indxr[l]] = a[k * n + indxc[l]];
                            a[k * n + indxr[l]] = dum;
                        }
                    }
                }
            }
    }
    }
      

  2.   

    谢谢 xiaopeipei2004(小裴),我要计算的矩阵是对称的,如下,在这种情况下有没有什么比高斯更快的算法。1139.0 4922.0  769.0  2620.0           10032.0 
    4922.0 33050.0 7201.0 15739.0          62027.8 
    769.0  7201.0  2293.0 4628.0           13981.5 
    2620.0 15739.0 4628.0 15062.0          34733.3
      

  3.   

    用LU分解法试试吧;:)
    public class d1r2F
    {
        void ludcmp(double a[], int n, int indx[], int d)
        {
            int nmax = 100;
            double tiny = 1e-20;
            double vv[] = new double[100];
            double sum, dum, aamax;
            int i, j, k, imax = 0;
            d = 1;
            for (i = 1; i <= n; i++)
            {
                aamax = 0.0;
                for (j = 1; j <= n; j++)
                {
                    if (Math.abs(a[(i - 1) * n + j]) > aamax)
                    {
                        aamax = Math.abs(a[(i - 1) * n + j]);
                    }
                }
                if (aamax == 0.0)
                {
                    System.out.println("singular matrix");
                    System.exit(1);
                }    
                vv[i] = 1.0 / aamax;
            }
            for (j = 1; j <= n; j++)
            {
                if (j > 1)
                {
                    for (i = 1; i <= j-1; i++)
                    {
                        sum = a[(i - 1) * n + j];
                        if (i > 1)
                        {
                            for (k = 1; k <= i - 1; k++)
                            {
                                sum = sum - a[(i - 1) * n + k] * a[(k - 1) * n + j];
                            }
                            a[(i - 1) * n + j] = sum;
                        }
                    }
                }
                aamax = 0.0;
                for (i = j; i <= n; i++)
                {
                    sum = a[(i - 1) * n + j];
                    if (j > 1)
                    {
                        for (k = 1; k <= j - 1; k++)
                        {
                             sum = sum - a[(i - 1) * n + k] * a[(k - 1) * n + j];
                        }
                        a[(i - 1) * n + j] = sum;
                    }
                    dum = vv[i] * Math.abs(sum);
                    if (dum >= aamax)
                    {
                        imax = i;
                        aamax = dum;
                    }
                }
                if (j != imax)
                {
                    for (k = 1; k <= n; k++)
                    {
                        dum = a[(imax - 1) * n + k];
                        a[(imax - 1) * n + k] = a[(j - 1) * n + k];
                        a[(j - 1) * n + k] = dum;
                    }
                    d = -d;
                    vv[imax] = vv[j];
                }
                indx[j] = imax;
                if (j != n)
                {
                    if (a[(j - 1) * n + j] == 0.0)
                    {
                        a[(j - 1) * n + j] = tiny;
                    }
                    dum = 1.0 / a[(j - 1) * n + j];
                    for (i = j + 1; i <= n; i++)
                    {
                        a[(i - 1) * n + j] = a[(i - 1) * n + j] * dum;
                    }
                }
            }
            if (a[(n - 1) * n + n] == 0.0)
            {
                a[(n - 1) * n + n] = tiny;
            }
        } 
    }
    /**
    还有一个class
    */
    public class d1r2F
    {
        void lubksb(double a[], int n, int indx[], double b[])
    {
    int i, j, ll;
    double sum;
            int ii = 0;
            for (i = 1; i <= n; i++)
    {
                ll = indx[i];
                sum = b[ll];
                b[ll] = b[i];
                if (ii != 0)
    {
                    for (j = ii; j<=i-1; j++)
    {
                        sum = sum - a[(i - 1) * n + j] * b[j];
                    }
    }
                else
    {
    if (sum != 0.0)
    {
    ii = i;
    }
                }
                b[i] = sum;
            }
            for (i = n; i >= 1; i--)
    {
                sum = b[i];
                if (i < n)
    {
                    for (j = i + 1; j <= n; j++)
    {
                        sum = sum - a[(i - 1) * n + j] * b[j];
                    }
                }
                b[i] = sum / a[(i - 1) * n + i];
            }
    }
    }