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 */
017package org.apache.commons.math4.linear;
018
019import org.apache.commons.math4.analysis.function.Sqrt;
020import org.apache.commons.math4.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 * @since 3.0
028 */
029public class JacobiPreconditioner extends RealLinearOperator {
030
031    /** The diagonal coefficients of the preconditioner. */
032    private final ArrayRealVector diag;
033
034    /**
035     * Creates a new instance of this class.
036     *
037     * @param diag the diagonal coefficients of the linear operator to be
038     * preconditioned
039     * @param deep {@code true} if a deep copy of the above array should be
040     * performed
041     */
042    public JacobiPreconditioner(final double[] diag, final boolean deep) {
043        this.diag = new ArrayRealVector(diag, deep);
044    }
045
046    /**
047     * Creates a new instance of this class. This method extracts the diagonal
048     * coefficients of the specified linear operator. If {@code a} does not
049     * extend {@link AbstractRealMatrix}, then the coefficients of the
050     * underlying matrix are not accessible, coefficient extraction is made by
051     * matrix-vector products with the basis vectors (and might therefore take
052     * some time). With matrices, direct entry access is carried out.
053     *
054     * @param a the linear operator for which the preconditioner should be built
055     * @return the diagonal preconditioner made of the inverse of the diagonal
056     * coefficients of the specified linear operator
057     * @throws NonSquareOperatorException if {@code a} is not square
058     */
059    public static JacobiPreconditioner create(final RealLinearOperator a)
060        throws NonSquareOperatorException {
061        final int n = a.getColumnDimension();
062        if (a.getRowDimension() != n) {
063            throw new NonSquareOperatorException(a.getRowDimension(), n);
064        }
065        final double[] diag = new double[n];
066        if (a instanceof AbstractRealMatrix) {
067            final AbstractRealMatrix m = (AbstractRealMatrix) a;
068            for (int i = 0; i < n; i++) {
069                diag[i] = m.getEntry(i, i);
070            }
071        } else {
072            final ArrayRealVector x = new ArrayRealVector(n);
073            for (int i = 0; i < n; i++) {
074                x.set(0.);
075                x.setEntry(i, 1.);
076                diag[i] = a.operate(x).getEntry(i);
077            }
078        }
079        return new JacobiPreconditioner(diag, false);
080    }
081
082    /** {@inheritDoc} */
083    @Override
084    public int getColumnDimension() {
085        return diag.getDimension();
086    }
087
088    /** {@inheritDoc} */
089    @Override
090    public int getRowDimension() {
091        return diag.getDimension();
092    }
093
094    /** {@inheritDoc} */
095    @Override
096    public RealVector operate(final RealVector x) {
097        // Dimension check is carried out by ebeDivide
098        return new ArrayRealVector(MathArrays.ebeDivide(x.toArray(),
099                                                        diag.toArray()),
100                                   false);
101    }
102
103    /**
104     * Returns the square root of {@code this} diagonal operator. More
105     * precisely, this method returns
106     * P = diag(1 / &radic;A<sub>11</sub>, 1 / &radic;A<sub>22</sub>, &hellip;).
107     *
108     * @return the square root of {@code this} preconditioner
109     * @since 3.1
110     */
111    public RealLinearOperator sqrt() {
112        final RealVector sqrtDiag = diag.map(new Sqrt());
113        return new RealLinearOperator() {
114            /** {@inheritDoc} */
115            @Override
116            public RealVector operate(final RealVector x) {
117                return new ArrayRealVector(MathArrays.ebeDivide(x.toArray(),
118                                                                sqrtDiag.toArray()),
119                                           false);
120            }
121
122            /** {@inheritDoc} */
123            @Override
124            public int getRowDimension() {
125                return sqrtDiag.getDimension();
126            }
127
128            /** {@inheritDoc} */
129            @Override
130            public int getColumnDimension() {
131                return sqrtDiag.getDimension();
132            }
133        };
134    }
135}