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>, …).
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 / √A<sub>11</sub>, 1 / √A<sub>22</sub>, …).
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 }