public void invert
        (double[][] Ainv)
        {
        // Verify preconditions.
        if (Ainv.length != N)
            {
            throw new IllegalArgumentException
                ("LinearSolve.invert(): Ainv is not a " + N + "-by-" + N +
                 " matrix");
            }
        for (int i = 0; i < N; ++ i)
            {
            if (Ainv[i].length != N)
                {
                throw new IllegalArgumentException
                    ("LinearSolve.invert(): Ainv is not a " + N + "-by-" + N +
                     " matrix");
                }
            }

        // Allocate temporary storage.
        double[] x = new double [N];

        // Compute and store the columns of the inverse.
        for (int i = 0; i < N; ++ i)
            {
            // Apply the permutation to column i of an identity matrix.
            for (int j = 0; j < N; ++ j)
                {
                x[j] = i == p[j] ? 1.0 : 0.0;
                }
            luSolve (LU, p, x, N);
            for (int j = 0; j < N; ++ j)
                {
                Ainv[j][i] = x[j];
                }
            }
        }