amg.java

另一个功能更强大的矩阵运算软件开源代码

            for (int i = 0; i < A.numRows(); ++i) {
                diagind[i] = no.uib.cipr.matrix.sparse.Arrays.binarySearch(
                        colind, i, rowptr[i], rowptr[i + 1]);
                if (diagind[i] < 0)
                    throw new RuntimeException(
                            "Matrix is missing a diagonal entry on row "
                                    + (i + 1));
            }

            return diagind;
        }

        /**
         * Finds the strongly coupled node neighborhoods
         */
        private List<Set<Integer>> findNodeNeighborhood(CompRowMatrix A,
                int[] diagind, double eps) {

            N = new ArrayList<Set<Integer>>(A.numRows());

            int[] rowptr = A.getRowPointers();
            int[] colind = A.getColumnIndices();
            double[] data = A.getData();

            for (int i = 0; i < A.numRows(); ++i) {
                Set<Integer> Ni = new HashSet<Integer>();

                double aii = data[diagind[i]];
                for (int j = rowptr[i]; j < rowptr[i + 1]; ++j) {
                    double aij = data[j];
                    double ajj = data[diagind[colind[j]]];

                    if (Math.abs(aij) >= eps * Math.sqrt(aii * ajj))
                        Ni.add(colind[j]);
                }

                N.add(Ni);
            }

            return N;
        }

        /**
         * Creates the initial R-set by including only the connected nodes
         */
        private boolean[] createInitialR(CompRowMatrix A) {
            boolean[] R = new boolean[A.numRows()];

            int[] rowptr = A.getRowPointers();
            int[] colind = A.getColumnIndices();
            double[] data = A.getData();

            for (int i = 0; i < A.numRows(); ++i) {
                boolean hasOffDiagonal = false;

                for (int j = rowptr[i]; j < rowptr[i + 1]; ++j)
                    if (colind[j] != i && data[j] != 0) {
                        hasOffDiagonal = true;
                        break;
                    }

                R[i] = hasOffDiagonal;
            }

            return R;
        }

        /**
         * Creates the initial aggregates
         */
        private List<Set<Integer>> createInitialAggregates(
                List<Set<Integer>> N, boolean[] R) {
            C = new ArrayList<Set<Integer>>();

            for (int i = 0; i < R.length; ++i) {

                // Skip non-free nodes
                if (!R[i])
                    continue;

                // See if all nodes in the current N-set are free
                boolean free = true;
                for (int j : N.get(i))
                    free &= R[j];

                // Create an aggregate out of N[i]
                if (free) {
                    C.add(new HashSet<Integer>(N.get(i)));
                    for (int j : N.get(i))
                        R[j] = false;
                }

            }

            return C;
        }

        /**
         * Enlarges the aggregates
         */
        private List<Set<Integer>> enlargeAggregates(List<Set<Integer>> C,
                List<Set<Integer>> N, boolean[] R) {

            // Contains the aggregates each node is coupled to
            List<List<Integer>> belong = new ArrayList<List<Integer>>(R.length);
            for (int i = 0; i < R.length; ++i)
                belong.add(new ArrayList<Integer>());

            // Find which aggregate each node is coupled to. This is used for
            // the intersection between Ni and Ck
            for (int k = 0; k < C.size(); ++k)
                for (int j : C.get(k))
                    belong.get(j).add(k);

            // Number of nodes in the intersection between each C and Ni
            int[] intersect = new int[C.size()];

            for (int i = 0; i < R.length; ++i) {

                // Skip non-free nodes
                if (!R[i])
                    continue;

                // Find the number of nodes intersecting Ni and every C, and
                // keep a track on the largest overlap
                Arrays.fill(intersect, 0);
                int largest = 0, maxValue = 0;
                for (int j : N.get(i))

                    // The k-index is to an aggregate coupled to node j
                    for (int k : belong.get(j)) {
                        intersect[k]++;
                        if (intersect[k] > maxValue) {
                            largest = k;
                            maxValue = intersect[largest];
                        }
                    }

                // Add the node to the proper C-set, and mark it as used
                // Also, check if the node actually does couple to a set
                if (maxValue > 0) {
                    R[i] = false;
                    C.get(largest).add(i);
                }
            }

            return C;
        }

        /**
         * Creates final aggregates from the remaining unallocated nodes
         */
        private List<Set<Integer>> createFinalAggregates(List<Set<Integer>> C,
                List<Set<Integer>> N, boolean[] R) {

            for (int i = 0; i < R.length; ++i) {

                // Skip non-free nodes
                if (!R[i])
                    continue;

                // Create new aggregate from the nodes in N[i] which are free
                Set<Integer> Cn = new HashSet<Integer>();
                for (int j : N.get(i))
                    if (R[j]) {
                        R[j] = false;
                        Cn.add(j);
                    }

                if (!Cn.isEmpty())
                    C.add(Cn);
            }

            return C;
        }
    }

    /**
     * Creates interpolation (prolongation) operators using based on the
     * aggregates. Can optionally smooth the aggregates
     */
    private static class Interpolator {

        /**
         * The Galerkin coarse-space operator
         */
        private CompRowMatrix Ac;

        /**
         * The interpolation (prolongation) matrix
         */
        private CompColMatrix I;

        /**
         * Creates the interpolation (prolongation) and Galerkin operators
         * 
         * @param aggregator
         *            Aggregates
         * @param A
         *            Matrix
         * @param omega
         *            Jacobi damping parameter between zero and one. If zero, no
         *            smoothing is performed, and a faster algorithm for forming
         *            the Galerkin operator will be used.
         */
        public Interpolator(Aggregator aggregator, CompRowMatrix A, double omega) {
            List<Set<Integer>> C = aggregator.getAggregates();
            List<Set<Integer>> N = aggregator.getNodeNeighborhoods();
            int[] diagind = aggregator.getDiagonalIndices();

            // Create the tentative prolongation, in compressed form
            int[] pt = createTentativeProlongation(C, A.numRows());

            /*
             * Apply Jacobi smoothing to the prolongator
             */

            if (omega != 0) {

                // Smooth the operator by a damped Jacobi method
                List<Map<Integer, Double>> P = createSmoothedProlongation(C, N,
                        A, diagind, omega, pt);

                // Form a compressed column storage for the operator
                I = createInterpolationMatrix(P, A.numRows());

                // Create the Galerkin operator using a slow method
                Ac = createGalerkinSlow(I, A);
            }

            /*
             * Use the aggregates as-is
             */

            else {

                // Create the Galerkin operator using a fast method
                Ac = createGalerkinFast(A, pt, C.size());

                // Form an explicit interpolation operator
                I = createInterpolationMatrix(pt, C.size());
            }
        }

        /**
         * Creates the tentative prolongation operator. Since the columns are
         * all disjoint, and its entries are binary, it is possible to store it
         * in a single array. Its length equals the number of fine nodes, and
         * the entries are the indices to the corresponding aggregate (C-set).
         */
        private int[] createTentativeProlongation(List<Set<Integer>> C, int n) {
            int[] pt = new int[n];
            Arrays.fill(pt, -1);

            for (int i = 0; i < C.size(); ++i)
                for (int j : C.get(i))
                    pt[j] = i;

            return pt;
        }

        /**
         * Creates the Galerkin operator using the assumption of disjoint
         * (non-smoothed) aggregates
         */
        private CompRowMatrix createGalerkinFast(CompRowMatrix A, int[] pt,
                int c) {
            int n = pt.length;

            FlexCompRowMatrix Ac = new FlexCompRowMatrix(c, c);

            int[] rowptr = A.getRowPointers();
            int[] colind = A.getColumnIndices();
            double[] data = A.getData();

            for (int i = 0; i < n; ++i)
                if (pt[i] != -1)
                    for (int j = rowptr[i]; j < rowptr[i + 1]; ++j)
                        if (pt[colind[j]] != -1)
                            Ac.add(pt[i], pt[colind[j]], data[j]);

            return new CompRowMatrix(Ac);
        }

        /**
         * Creates the interpolation (prolongation) matrix based on the smoothed
         * aggregates
         */
        private CompColMatrix createInterpolationMatrix(
                List<Map<Integer, Double>> P, int n) {

            // Determine the sparsity pattern of I
            int c = P.size();
            int[][] nz = new int[c][];
            for (int j = 0; j < c; ++j) {

                Map<Integer, Double> Pj = P.get(j);
                nz[j] = new int[Pj.size()];

                int l = 0;
                for (int k : Pj.keySet())
                    nz[j][l++] = k;
            }

            I = new CompColMatrix(n, c, nz);

            // Populate it with numerical entries
            for (int j = 0; j < c; ++j) {

                Map<Integer, Double> Pj = P.get(j);

                for (Map.Entry<Integer, Double> e : Pj.entrySet())
                    I.set(e.getKey(), j, e.getValue());
            }

            return I;
        }

        /**
         * Creates the interpolation (prolongation) matrix based on the
         * non-smoothed aggregates
         */
        private CompColMatrix createInterpolationMatrix(int[] pt, int c) {
            FlexCompColMatrix If = new FlexCompColMatrix(pt.length, c);

            for (int i = 0; i < pt.length; ++i)
                if (pt[i] != -1)
                    If.set(i, pt[i], 1);

            return new CompColMatrix(If);
        }

        /**
         * Gets the interpolation (prolongation) operator
         */
        public CompColMatrix getInterpolationOperator() {
            return I;
        }

        /**
         * Creates the smoothes interpolation (prolongation) operator by a
         * single sweep of the damped Jacobi method
         */
        private List<Map<Integer, Double>> createSmoothedProlongation(
                List<Set<Integer>> C, List<Set<Integer>> N, CompRowMatrix A,
                int[] diagind, double omega, int[] pt) {

            int n = A.numRows(), c = C.size();

            // Allocate the interpolation (prolongation) operator
            // It is stored by columns, so the maps take row-indices as keys
            List<Map<Integer, Double>> P = new ArrayList<Map<Integer, Double>>(
                    c);
            for (int i = 0; i < c; ++i)
                P.add(new HashMap<Integer, Double>());

            int[] rowptr = A.getRowPointers();
            int[] colind = A.getColumnIndices();
            double[] data = A.getData();

            double[] dot = new double[c];

            // Apply the damped Jacobi smoother
            for (int i = 0; i < n; ++i) {

                if (pt[i] == -1)
                    continue;

                Arrays.fill(dot, 0);
                Set<Integer> Ni = N.get(i);

                // Calculate A*Pt, except for the diagonal
                double weakAij = 0;
                for (int j = rowptr[i]; j < rowptr[i + 1]; ++j) {

                    if (pt[colind[j]] == -1)
                        continue;

                    double aij = data[j];

                    // Off-diagonal, include only strong couplings, and add the
                    // weak couplings to the diagonal
                    if (aij != 0 && !Ni.contains(colind[j])) {
                        weakAij += aij;
                        continue;
                    }

                    dot[pt[colind[j]]] += aij;
                }

                // Subtract the weak couplings from the diagonal part of A*Pt
                dot[pt[i]] -= weakAij;

                // Scale by omega and the inverse of the diagonal (damping)
                double scale = -omega / data[diagind[i]];
                for (int j = 0; j < dot.length; ++j)
                    dot[j] *= scale;

                // Set to (I-omega*D^{-1}*A)*Pt
                dot[pt[i]]++;

                // This has formed a whole row of P=(I-omega*D^{-1}*A)*Pt
                // Store the non-zeros into the sparse structure
                for (int j = 0; j < dot.length; ++j)
                    if (dot[j] != 0)
                        P.get(j).put(i, dot[j]);
            }

            return P;
        }

        /**
         * Creates the entries of the Galerkin operator
         * <code>Ac = I<sup>T</sup> A I</code>. This is a very
         * time-consuming operation
         */
        private CompRowMatrix createGalerkinSlow(CompColMatrix I,
                CompRowMatrix A) {
            int n = I.numRows(), c = I.numColumns();
            FlexCompRowMatrix Ac = new FlexCompRowMatrix(c, c);

            double[] aiCol = new double[n];
            double[] iCol = new double[n];
            DenseVector aiV = new DenseVector(aiCol, false);
            DenseVector iV = new DenseVector(iCol, false);
            double[] itaiCol = new double[c];
            DenseVector itaiV = new DenseVector(itaiCol, false);

            int[] colptr = I.getColumnPointers();
            int[] rowind = I.getRowIndices();
            double[] Idata = I.getData();

            for (int k = 0; k < c; ++k) {

                // Expand column 'k' of I to dense storage
                iV.zero();
                for (int i = colptr[k]; i < colptr[k + 1]; ++i)
                    iCol[rowind[i]] = Idata[i];

                // Form column 'k' of A*I
                A.mult(iV, aiV);

                // Form column 'k' of I'*A*I
                I.transMult(aiV, itaiV);

                // Store non-zeros into Ac
                for (int i = 0; i < c; ++i)
                    if (itaiCol[i] != 0)
                        Ac.set(i, k, itaiCol[i]);
            }

            return new CompRowMatrix(Ac);
        }

        /**
         * Gets the Galerkin operator
         */
        public CompRowMatrix getGalerkinOperator() {
            return Ac;
        }

    }
}