001// --- BEGIN LICENSE BLOCK ---
002/* 
003 * Copyright (c) 2009, Mikio L. Braun
004 * All rights reserved.
005 * 
006 * Redistribution and use in source and binary forms, with or without
007 * modification, are permitted provided that the following conditions are
008 * met:
009 * 
010 *     * Redistributions of source code must retain the above copyright
011 *       notice, this list of conditions and the following disclaimer.
012 * 
013 *     * Redistributions in binary form must reproduce the above
014 *       copyright notice, this list of conditions and the following
015 *       disclaimer in the documentation and/or other materials provided
016 *       with the distribution.
017 * 
018 *     * Neither the name of the Technische Universit?t Berlin nor the
019 *       names of its contributors may be used to endorse or promote
020 *       products derived from this software without specific prior
021 *       written permission.
022 * 
023 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
024 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
025 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
026 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
027 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
028 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
029 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
030 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
031 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
032 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
033 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
034 */
035// --- END LICENSE BLOCK ---
036
037package org.jblas;
038
039import org.jblas.util.Functions;
040
041/**
042 * Solving linear equations.
043 */
044public class Solve {
045        /** Solves the linear equation A*X = B. */
046        public static DoubleMatrix solve(DoubleMatrix A, DoubleMatrix B) {
047                A.assertSquare();
048                DoubleMatrix X = B.dup();
049                int[] ipiv = new int[B.rows];
050                SimpleBlas.gesv(A.dup(), ipiv, X);
051                return X;
052        }
053
054        /** Solves the linear equation A*X = B for symmetric A. */
055        public static DoubleMatrix solveSymmetric(DoubleMatrix A, DoubleMatrix B) {
056                A.assertSquare();
057                DoubleMatrix X = B.dup();
058                int[] ipiv = new int[B.rows];
059                SimpleBlas.sysv('U', A.dup(), ipiv, X);
060                return X;
061        }
062
063        
064        /** Solves the linear equation A*X = B for symmetric and positive definite A. */
065        public static DoubleMatrix solvePositive(DoubleMatrix A, DoubleMatrix B) {
066                A.assertSquare();
067                DoubleMatrix X = B.dup();
068                SimpleBlas.posv('U', A.dup(), X);
069                return X;
070        }
071
072  /** Computes the Least Squares solution for over or underdetermined
073   * linear equations A*X = B
074   *
075   * In the overdetermined case, when m > n, that is, there are more equations than
076   * variables, it computes the least squares solution of X -> ||A*X - B ||_2.
077   *
078   * In the underdetermined case, when m < n (less equations than variables), there are infinitely
079   * many solutions and it computes the minimum norm solution.
080   *
081   * @param A an (m,n) matrix
082   * @param B a (m,k) matrix
083   * @return either the minimum norm or least squares solution.
084   */
085  public static DoubleMatrix solveLeastSquares(DoubleMatrix A, DoubleMatrix B) {
086    if (B.rows < A.columns) {
087      DoubleMatrix X = DoubleMatrix.concatVertically(B, new DoubleMatrix(A.columns - B.rows, B.columns));
088      SimpleBlas.gelsd(A.dup(), X);
089      return X;
090    } else {
091      DoubleMatrix X = B.dup();
092      SimpleBlas.gelsd(A.dup(), X);
093      return X.getRange(0, A.columns, 0, B.columns);
094    }
095  }
096
097  /**
098   * Computes the pseudo-inverse.
099   *
100   * Note, this function uses the solveLeastSquares and might produce different numerical
101   * solutions for the underdetermined case than matlab.
102   *
103   * @param A rectangular matrix
104   * @return matrix P such that A*P*A = A and P*A*P = P.
105   */
106  public static DoubleMatrix pinv(DoubleMatrix A) {
107    return solveLeastSquares(A, DoubleMatrix.eye(A.rows));
108  }
109
110//BEGIN
111  // The code below has been automatically generated.
112  // DO NOT EDIT!
113        /** Solves the linear equation A*X = B. */
114        public static FloatMatrix solve(FloatMatrix A, FloatMatrix B) {
115                A.assertSquare();
116                FloatMatrix X = B.dup();
117                int[] ipiv = new int[B.rows];
118                SimpleBlas.gesv(A.dup(), ipiv, X);
119                return X;
120        }
121
122        /** Solves the linear equation A*X = B for symmetric A. */
123        public static FloatMatrix solveSymmetric(FloatMatrix A, FloatMatrix B) {
124                A.assertSquare();
125                FloatMatrix X = B.dup();
126                int[] ipiv = new int[B.rows];
127                SimpleBlas.sysv('U', A.dup(), ipiv, X);
128                return X;
129        }
130
131        
132        /** Solves the linear equation A*X = B for symmetric and positive definite A. */
133        public static FloatMatrix solvePositive(FloatMatrix A, FloatMatrix B) {
134                A.assertSquare();
135                FloatMatrix X = B.dup();
136                SimpleBlas.posv('U', A.dup(), X);
137                return X;
138        }
139
140  /** Computes the Least Squares solution for over or underdetermined
141   * linear equations A*X = B
142   *
143   * In the overdetermined case, when m > n, that is, there are more equations than
144   * variables, it computes the least squares solution of X -> ||A*X - B ||_2.
145   *
146   * In the underdetermined case, when m < n (less equations than variables), there are infinitely
147   * many solutions and it computes the minimum norm solution.
148   *
149   * @param A an (m,n) matrix
150   * @param B a (m,k) matrix
151   * @return either the minimum norm or least squares solution.
152   */
153  public static FloatMatrix solveLeastSquares(FloatMatrix A, FloatMatrix B) {
154    if (B.rows < A.columns) {
155      FloatMatrix X = FloatMatrix.concatVertically(B, new FloatMatrix(A.columns - B.rows, B.columns));
156      SimpleBlas.gelsd(A.dup(), X);
157      return X;
158    } else {
159      FloatMatrix X = B.dup();
160      SimpleBlas.gelsd(A.dup(), X);
161      return X.getRange(0, A.columns, 0, B.columns);
162    }
163  }
164
165  /**
166   * Computes the pseudo-inverse.
167   *
168   * Note, this function uses the solveLeastSquares and might produce different numerical
169   * solutions for the underdetermined case than matlab.
170   *
171   * @param A rectangular matrix
172   * @return matrix P such that A*P*A = A and P*A*P = P.
173   */
174  public static FloatMatrix pinv(FloatMatrix A) {
175    return solveLeastSquares(A, FloatMatrix.eye(A.rows));
176  }
177
178//END
179}