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    }