diff --git a/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java b/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java index ab0e1ac38..b9a2df97f 100644 --- a/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java +++ b/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java @@ -187,6 +187,8 @@ public enum LocalizedFormats implements Localizable { NOT_POSITIVE_COLUMNDIMENSION("invalid column dimension: {0} (must be positive)"), NOT_POSITIVE_DEFINITE_MATRIX("not positive definite matrix"), NON_POSITIVE_DEFINITE_MATRIX("not positive definite matrix: diagonal element at ({0},{0}) is larger than {2}"), /* keep */ + NON_POSITIVE_DEFINITE_LINEAR_OPERATOR("non positive definite linear operator: x' A x <= 0 when x is {0}"), + NON_SELF_ADJOINT_LINEAR_OPERATOR("non self-adjoint linear operator: |x' A y - y' A x| > {0} when x is {1} and y is {2}"), DEGREES_OF_FREEDOM("degrees of freedom ({0})"), /* keep */ NOT_POSITIVE_DEGREES_OF_FREEDOM("degrees of freedom must be positive ({0})"), NOT_POSITIVE_ELEMENT_AT_INDEX("element {0} is not positive: {1}"), diff --git a/src/main/java/org/apache/commons/math/linear/NonPositiveDefiniteLinearOperatorException.java b/src/main/java/org/apache/commons/math/linear/NonPositiveDefiniteLinearOperatorException.java new file mode 100644 index 000000000..43d1b5c0c --- /dev/null +++ b/src/main/java/org/apache/commons/math/linear/NonPositiveDefiniteLinearOperatorException.java @@ -0,0 +1,86 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.linear; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; + +/** + * Exception to be thrown when a symmetric, definite positive + * {@link RealLinearOperator} is expected. + * Since the coefficients of the matrix are not accessible, the most + * general definition is used to check that {@code A} is not positive + * definite, i.e. there exists {@code x} such that {@code x' A x <= 0}. + * In the terminology of this exception, {@code A} is the "offending" + * linear operator and {@code x} the "offending" vector. + * + * @version $Id$ + */ +public class NonPositiveDefiniteLinearOperatorException + extends MathIllegalArgumentException { + /** The offending linear operator.*/ + private final RealLinearOperator a; + /** A reference to the offending vector. */ + private final RealVector x; + + /** + * Creates a new instance of this class. + * + * @param a Offending linear operator. + * @param x Offending vector. + */ + public NonPositiveDefiniteLinearOperatorException(final RealLinearOperator a, + final double[] x) { + this(a, new ArrayRealVector(x, false)); + } + + /** + * Creates a new instance of this class. + * + * @param a Offending linear operator. + * @param x Offending vector. + */ + public NonPositiveDefiniteLinearOperatorException(final RealLinearOperator a, + final RealVector x) { + super(LocalizedFormats.NON_POSITIVE_DEFINITE_LINEAR_OPERATOR, x); + this.a = a; + this.x = x; + } + + /** + * Returns a reference to the offending vector. + * If the exception was raised by a call to + * {@link #NonPositiveDefiniteLinearOperatorException(RealLinearOperator, + * double[])}, then a new {@link ArrayRealVector} holding a reference to + * the actual {@code double[]} is returned. + * + * @return the offending vector. + */ + public RealVector copyOffendingVector() { + return x; + } + + /** + * Returns a reference to the offending linear operator. + * + * @return the offending linear operator. + */ + public RealLinearOperator getOffendingLinearOperator() { + return a; + } +} diff --git a/src/main/java/org/apache/commons/math/linear/NonSelfAdjointLinearOperatorException.java b/src/main/java/org/apache/commons/math/linear/NonSelfAdjointLinearOperatorException.java new file mode 100644 index 000000000..578bfcd6f --- /dev/null +++ b/src/main/java/org/apache/commons/math/linear/NonSelfAdjointLinearOperatorException.java @@ -0,0 +1,128 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.linear; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; + +/** + * Exception to be thrown when a self-adjoint {@link RealLinearOperator} + * is expected. + * Since the coefficients of the matrix are not accessible, the most + * general definition is used to check that A is not self-adjoint, i.e. + * there exist x and y such as {@code | x' A y - y' A x | >= eps}, + * where {@code eps} is a user-specified tolerance, and {@code x'} + * denotes the transpose of {@code x}. + * In the terminology of this exception, {@code A} is the "offending" + * linear operator, {@code x} and {@code y} are the first and second + * "offending" vectors, respectively. + * + * @version $Id$ + */ +public class NonSelfAdjointLinearOperatorException + extends MathIllegalArgumentException { + /** The offending linear operator, A. */ + private final RealLinearOperator a; + /** The threshold. */ + private final double threshold; + /** A reference to the first offending vector*/ + private final RealVector x; + /** A reference to the second offending vector*/ + private final RealVector y; + + /** + * Creates a new instance of this class. + * + * @param a Offending linear operator. + * @param x First offending vector. + * @param y Second offending vector. + * @param threshold Threshold. + */ + public NonSelfAdjointLinearOperatorException(final RealLinearOperator a, + final double[] x, + final double[] y, + final double threshold) { + this(a, + new ArrayRealVector(x, false), + new ArrayRealVector(y, false), + threshold); + } + + /** + * Creates a new instance of this class. + * + * @param a Offending linear operator. + * @param x First offending vector. + * @param y Second offending vector. + * @param threshold Threshold. + */ + public NonSelfAdjointLinearOperatorException(final RealLinearOperator a, + final RealVector x, + final RealVector y, + final double threshold) { + super(LocalizedFormats.NON_SELF_ADJOINT_LINEAR_OPERATOR, threshold, x, y); + this.a = a; + this.x = x; + this.y = y; + this.threshold = threshold; + } + + /** + * Returns a reference to the first offending vector. + * If the exception was raised by a call to + * {@link #NonSelfAdjointLinearOperatorException(RealLinearOperator, + * double[], double[], double)}, then a new {@link ArrayRealVector} + * holding a reference to the actual {@code double[]} is returned. + * + * @return the first offending vector. + */ + public RealVector getFirstOffendingVector() { + return x; + } + + /** + * Returns a reference to the offending linear operator. + * + * @return the offending linear operator. + */ + public RealLinearOperator getOffendingLinearOperator() { + return a; + } + + /** + * Returns a copy of the second offending vector. + * If the exception was raised by a call to + * {@link #NonSelfAdjointLinearOperatorException(RealLinearOperator, + * double[], double[], double)}, then a new {@link ArrayRealVector} + * holding a reference to the actual {@code double[]} is returned. + * + * @return the second offending vector. + */ + public RealVector getSecondOffendingVector() { + return y; + } + + /** + * Returns the threshold. + * + * @return the threshold. + */ + public double getThreshold() { + return threshold; + } +} diff --git a/src/main/java/org/apache/commons/math/linear/RealLinearOperator.java b/src/main/java/org/apache/commons/math/linear/RealLinearOperator.java new file mode 100644 index 000000000..9ed3ed2d5 --- /dev/null +++ b/src/main/java/org/apache/commons/math/linear/RealLinearOperator.java @@ -0,0 +1,113 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.linear; + +import org.apache.commons.math.exception.DimensionMismatchException; + +/** + * This class defines a linear operator operating on real ({@code double}) + * vector spaces. + * No direct access to the coefficients of the underlying matrix is provided. + * + * The motivation for such an interface is well stated by + * Barrett et al. (1994): + *
+ * We restrict ourselves to iterative methods, which work by repeatedly + * improving an approximate solution until it is accurate enough. These + * methods access the coefficient matrix {@code A} of the linear system + * only via the matrix-vector product {@code y = A x} (and perhaps + * {@code z} = {@code A}T {@code x}). Thus the user need only + * supply a subroutine for computing {@code y} (and perhaps {@code z}) + * given {@code x}, which permits full exploitation of the sparsity or + * other special structure of A. + *
+ *
+ * + *
+ *
Barret et al. (1994)
+ *
+ * R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra, V. + * Eijkhout, R. Pozo, C. Romine and H. Van der Vorst, + * Templates for the Solution of Linear Systems: Building Blocks for + * Iterative Methods, SIAM
+ *
+ * + * @version $Id$ + */ +public abstract class RealLinearOperator { + /** The dimension of the codomain. */ + private final int rowDimension; + /** The dimension of the domain. */ + private final int columnDimension; + + /** + * Creates a new instance of this class, with specified dimensions + * of the domain and codomain. + * + * @param rowDimension Dimension of the codomain (number of rows). + * @param columnDimension Dimension of the domain (number of columns). + */ + public RealLinearOperator(final int rowDimension, final int columnDimension) { + this.columnDimension = columnDimension; + this.rowDimension = rowDimension; + } + + /** + * Returns the dimension of the codomain of this operator. + * + * @return the number of rows of the underlying matrix. + */ + public final int getRowDimension() { + return rowDimension; + } + + /** + * Returns the dimension of the domain of this operator. + * + * @return the number of columns of the underlying matrix. + */ + public final int getColumnDimension() { + return columnDimension; + } + + /** + * Returns the result of multiplying {@code this} by the vector {@code x}. + * + * @param x Vector to operate on. + * @return the product of {@code this} instance with {@code x}. + */ + public double[] operate(final double[] x) { + if (x.length != getColumnDimension()) { + throw new DimensionMismatchException(x.length, getColumnDimension()); + } + final RealVector y = operate(new ArrayRealVector(x, false)); + if (y instanceof ArrayRealVector) { + return ((ArrayRealVector) y).getDataRef(); + } else { + return y.getData(); + } + } + + /** + * Returns the result of multiplying {@code this} by the vector {@code x}. + * + * @param x Vector to operate on. + * @return the product of {@code this} instance with {@code x}. + */ + public abstract RealVector operate(final RealVector x); +} diff --git a/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties b/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties index e59503b93..1b34bd61b 100644 --- a/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties +++ b/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties @@ -157,6 +157,8 @@ BETA = beta NOT_POSITIVE_COLUMNDIMENSION = nombre de colonnes invalide : {0} (doit \u00eatre positif) NOT_POSITIVE_DEFINITE_MATRIX = matrice non d\u00e9finie positive NON_POSITIVE_DEFINITE_MATRIX = matrice non d\u00e9finie positive: l''\u00e9l\u00e9ment diagonal ({0},{0}) est sup\u00e9rieur \u0061 {1} +NON_POSITIVE_DEFINITE_LINEAR_OPERATOR = op\u00e9 lin\u00e9aire non d\u00e9fini positif: x'' A x <= 0 quand x est {0} +NON_SELF_ADJOINT_LINEAR_OPERATOR = op\u00e9 lin\u00e9aire non auto-adjoint: |x'' A y - y'' A x| > {0} quand x est {1} et y est {2} NOT_POSITIVE_DEGREES_OF_FREEDOM = les degr\u00e9s de libert\u00e9 doivent \u00eatre positifs ({0}) DEGREES_OF_FREEDOM = degr\u00e9s de libert\u00e9 ({0}) NOT_POSITIVE_ELEMENT_AT_INDEX = l''\u00e9l\u00e9ment {0} n''est pas positif : {1} diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 28a9091ed..34c63692e 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,9 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Framework for iterative linear solvers. + Improved efficiency in RandomDataImpl, LaguerreSolver, FastMath and OutlineExtractor by moving conditional code into blocks where it is needed.