    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.linear;
    018
    019    import org.apache.commons.math3.analysis.function.Sqrt;
    020    import org.apache.commons.math3.util.MathArrays;
    021
    022    /**
    023     * This class implements the standard Jacobi (diagonal) preconditioner. For a
    024     * matrix A<sub>ij</sub>, this preconditioner is
    025     * M = diag(1 / A<sub>11</sub>, 1 / A<sub>22</sub>, &hellip;).
    026     *
    027     * @version $Id: JacobiPreconditioner.java 1422195 2012-12-15 06:45:18Z psteitz $
    028     * @since 3.0
    029     */
    030    public class JacobiPreconditioner extends RealLinearOperator {
    031
    032        /** The diagonal coefficients of the preconditioner. */
    033        private final ArrayRealVector diag;
    034
    035        /**
    036         * Creates a new instance of this class.
    037         *
    038         * @param diag the diagonal coefficients of the linear operator to be
    039         * preconditioned
    040         * @param deep {@code true} if a deep copy of the above array should be
    041         * performed
    042         */
    043        public JacobiPreconditioner(final double[] diag, final boolean deep) {
    044            this.diag = new ArrayRealVector(diag, deep);
    045        }
    046
    047        /**
    048         * Creates a new instance of this class. This method extracts the diagonal
    049         * coefficients of the specified linear operator. If {@code a} does not
    050         * extend {@link AbstractRealMatrix}, then the coefficients of the
    051         * underlying matrix are not accessible, coefficient extraction is made by
    052         * matrix-vector products with the basis vectors (and might therefore take
    053         * some time). With matrices, direct entry access is carried out.
    054         *
    055         * @param a the linear operator for which the preconditioner should be built
    056         * @return the diagonal preconditioner made of the inverse of the diagonal
    057         * coefficients of the specified linear operator
    058         * @throws NonSquareOperatorException if {@code a} is not square
    059         */
    060        public static JacobiPreconditioner create(final RealLinearOperator a)
    061            throws NonSquareOperatorException {
    062            final int n = a.getColumnDimension();
    063            if (a.getRowDimension() != n) {
    064                throw new NonSquareOperatorException(a.getRowDimension(), n);
    065            }
    066            final double[] diag = new double[n];
    067            if (a instanceof AbstractRealMatrix) {
    068                final AbstractRealMatrix m = (AbstractRealMatrix) a;
    069                for (int i = 0; i < n; i++) {
    070                    diag[i] = m.getEntry(i, i);
    071                }
    072            } else {
    073                final ArrayRealVector x = new ArrayRealVector(n);
    074                for (int i = 0; i < n; i++) {
    075                    x.set(0.);
    076                    x.setEntry(i, 1.);
    077                    diag[i] = a.operate(x).getEntry(i);
    078                }
    079            }
    080            return new JacobiPreconditioner(diag, false);
    081        }
    082
    083        /** {@inheritDoc} */
    084        @Override
    085        public int getColumnDimension() {
    086            return diag.getDimension();
    087        }
    088
    089        /** {@inheritDoc} */
    090        @Override
    091        public int getRowDimension() {
    092            return diag.getDimension();
    093        }
    094
    095        /** {@inheritDoc} */
    096        @Override
    097        public RealVector operate(final RealVector x) {
    098            // Dimension check is carried out by ebeDivide
    099            return new ArrayRealVector(MathArrays.ebeDivide(x.toArray(),
    100                                                            diag.toArray()),
    101                                       false);
    102        }
    103
    104        /**
    105         * Returns the square root of {@code this} diagonal operator. More
    106         * precisely, this method returns
    107         * P = diag(1 / &radic;A<sub>11</sub>, 1 / &radic;A<sub>22</sub>, &hellip;).
    108         *
    109         * @return the square root of {@code this} preconditioner
    110         * @since 3.1
    111         */
    112        public RealLinearOperator sqrt() {
    113            final RealVector sqrtDiag = diag.map(new Sqrt());
    114            return new RealLinearOperator() {
    115                /** {@inheritDoc} */
    116                @Override
    117                public RealVector operate(final RealVector x) {
    118                    return new ArrayRealVector(MathArrays.ebeDivide(x.toArray(),
    119                                                                    sqrtDiag.toArray()),
    120                                               false);
    121                }
    122
    123                /** {@inheritDoc} */
    124                @Override
    125                public int getRowDimension() {
    126                    return sqrtDiag.getDimension();
    127                }
    128
    129                /** {@inheritDoc} */
    130                @Override
    131                public int getColumnDimension() {
    132                    return sqrtDiag.getDimension();
    133                }
    134            };
    135        }
    136    }