1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math4.legacy.linear;
18
19 import org.apache.commons.math4.legacy.analysis.function.Sqrt;
20 import org.apache.commons.math4.legacy.core.MathArrays;
21
22 /**
23 * This class implements the standard Jacobi (diagonal) preconditioner. For a
24 * matrix A<sub>ij</sub>, this preconditioner is
25 * M = diag(1 / A<sub>11</sub>, 1 / A<sub>22</sub>, …).
26 *
27 * @since 3.0
28 */
29 public class JacobiPreconditioner extends RealLinearOperator {
30
31 /** The diagonal coefficients of the preconditioner. */
32 private final ArrayRealVector diag;
33
34 /**
35 * Creates a new instance of this class.
36 *
37 * @param diag the diagonal coefficients of the linear operator to be
38 * preconditioned
39 * @param deep {@code true} if a deep copy of the above array should be
40 * performed
41 */
42 public JacobiPreconditioner(final double[] diag, final boolean deep) {
43 this.diag = new ArrayRealVector(diag, deep);
44 }
45
46 /**
47 * Creates a new instance of this class. This method extracts the diagonal
48 * coefficients of the specified linear operator. If {@code a} does not
49 * extend {@link AbstractRealMatrix}, then the coefficients of the
50 * underlying matrix are not accessible, coefficient extraction is made by
51 * matrix-vector products with the basis vectors (and might therefore take
52 * some time). With matrices, direct entry access is carried out.
53 *
54 * @param a the linear operator for which the preconditioner should be built
55 * @return the diagonal preconditioner made of the inverse of the diagonal
56 * coefficients of the specified linear operator
57 * @throws NonSquareOperatorException if {@code a} is not square
58 */
59 public static JacobiPreconditioner create(final RealLinearOperator a)
60 throws NonSquareOperatorException {
61 final int n = a.getColumnDimension();
62 if (a.getRowDimension() != n) {
63 throw new NonSquareOperatorException(a.getRowDimension(), n);
64 }
65 final double[] diag = new double[n];
66 if (a instanceof AbstractRealMatrix) {
67 final AbstractRealMatrix m = (AbstractRealMatrix) a;
68 for (int i = 0; i < n; i++) {
69 diag[i] = m.getEntry(i, i);
70 }
71 } else {
72 final ArrayRealVector x = new ArrayRealVector(n);
73 for (int i = 0; i < n; i++) {
74 x.set(0.);
75 x.setEntry(i, 1.);
76 diag[i] = a.operate(x).getEntry(i);
77 }
78 }
79 return new JacobiPreconditioner(diag, false);
80 }
81
82 /** {@inheritDoc} */
83 @Override
84 public int getColumnDimension() {
85 return diag.getDimension();
86 }
87
88 /** {@inheritDoc} */
89 @Override
90 public int getRowDimension() {
91 return diag.getDimension();
92 }
93
94 /** {@inheritDoc} */
95 @Override
96 public RealVector operate(final RealVector x) {
97 // Dimension check is carried out by ebeDivide
98 return new ArrayRealVector(MathArrays.ebeDivide(x.toArray(),
99 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 / √A<sub>11</sub>, 1 / √A<sub>22</sub>, …).
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 }