private static double luDecompose
        (double[][] A,
         int[] p,
         int N)
        {
        // Initialize sign of permutation.
        double signum = 1.0;

        // Initialize permutation.
        for (int i = 0; i < N; ++ i)
            {
            p[i] = i;
            }

        // Do all columns.
        for (int j = 0; j < N-1; ++ j)
            {
            // Find pivot element (maximum element) in the j-th column.
            double max = Math.abs (A[j][j]);
            int i_pivot = j;
            for (int i = j+1; i < N; ++ i)
                {
                double aij = Math.abs (A[i][j]);
                if (aij > max)
                    {
                    max = aij;
                    i_pivot = i;
                    }
                }

            // If the pivot element is not on the diagonal, interchange rows.
            if (i_pivot != j)
                {
                // Swap pivot row with diagonal row.
                double[] swap = A[i_pivot];
                A[i_pivot] = A[j];
                A[j] = swap;

                // Update permutation.
                int swap2 = p[i_pivot];
                p[i_pivot] = p[j];
                p[j] = swap2;
                signum = - signum;
                }

            // Update the decomposition.
            double ajj = A[j][j];
            if (ajj != 0.0)
                {
                for (int i = j+1; i < N; ++ i)
                    {
                    double aij = A[i][j] / ajj;
                    A[i][j] = aij;
                    for (int k = j+1; k < N; ++ k)
                        {
                        A[i][k] -= aij * A[j][k];
                        }
                    }
                }
            else // (ajj == 0.0)
                {
                throw new RuntimeException
                    ("LinearSolve(): Zero pivot encountered");
                }
            }

        // Return sign of permutation.
        return signum;
        }