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
039/**
040 * Solving linear equations.
041 */
042public class Solve {
043        /** Solves the linear equation A*X = B. */
044        public static DoubleMatrix solve(DoubleMatrix A, DoubleMatrix B) {
045                A.assertSquare();
046                DoubleMatrix X = B.dup();
047                int[] ipiv = new int[B.rows];
048                SimpleBlas.gesv(A.dup(), ipiv, X);
049                return X;
050        }
051
052        /** Solves the linear equation A*X = B for symmetric A. */
053        public static DoubleMatrix solveSymmetric(DoubleMatrix A, DoubleMatrix B) {
054                A.assertSquare();
055                DoubleMatrix X = B.dup();
056                int[] ipiv = new int[B.rows];
057                SimpleBlas.sysv('U', A.dup(), ipiv, X);
058                return X;
059        }
060
061        
062        /** Solves the linear equation A*X = B for symmetric and positive definite A. */
063        public static DoubleMatrix solvePositive(DoubleMatrix A, DoubleMatrix B) {
064                A.assertSquare();
065                DoubleMatrix X = B.dup();
066                SimpleBlas.posv('U', A.dup(), X);
067                return X;
068        }
069
070  /** Computes the Least Squares solution for over or underdetermined
071   * linear equations A*X = B
072   *
073   * In the overdetermined case, when m > n, that is, there are more equations than
074   * variables, it computes the least squares solution of X -> ||A*X - B ||_2.
075   *
076   * In the underdetermined case, when m < n (less equations than variables), there are infinitely
077   * many solutions and it computes the minimum norm solution.
078   *
079   * @param A an (m,n) matrix
080   * @param B a (m,k) matrix
081   * @return either the minimum norm or least squares solution.
082   */
083  public static DoubleMatrix solveLeastSquares(DoubleMatrix A, DoubleMatrix B) {
084    if (B.rows < A.columns) {
085      DoubleMatrix X = DoubleMatrix.concatVertically(B, new DoubleMatrix(A.columns - B.rows, B.columns));
086      SimpleBlas.gelsd(A.dup(), X);
087      return X;
088    } else {
089      DoubleMatrix X = B.dup();
090      SimpleBlas.gelsd(A.dup(), X);
091      return X.getRange(0, A.columns, 0, B.columns);
092    }
093  }
094
095  /**
096   * Computes the pseudo-inverse.
097   *
098   * Note, this function uses the solveLeastSquares and might produce different numerical
099   * solutions for the underdetermined case than matlab.
100   *
101   * @param A rectangular matrix
102   * @return matrix P such that A*P*A = A and P*A*P = P.
103   */
104  public static DoubleMatrix pinv(DoubleMatrix A) {
105    return solveLeastSquares(A, DoubleMatrix.eye(A.rows));
106  }
107
108//BEGIN
109  // The code below has been automatically generated.
110  // DO NOT EDIT!
111        /** Solves the linear equation A*X = B. */
112        public static FloatMatrix solve(FloatMatrix A, FloatMatrix B) {
113                A.assertSquare();
114                FloatMatrix X = B.dup();
115                int[] ipiv = new int[B.rows];
116                SimpleBlas.gesv(A.dup(), ipiv, X);
117                return X;
118        }
119
120        /** Solves the linear equation A*X = B for symmetric A. */
121        public static FloatMatrix solveSymmetric(FloatMatrix A, FloatMatrix B) {
122                A.assertSquare();
123                FloatMatrix X = B.dup();
124                int[] ipiv = new int[B.rows];
125                SimpleBlas.sysv('U', A.dup(), ipiv, X);
126                return X;
127        }
128
129        
130        /** Solves the linear equation A*X = B for symmetric and positive definite A. */
131        public static FloatMatrix solvePositive(FloatMatrix A, FloatMatrix B) {
132                A.assertSquare();
133                FloatMatrix X = B.dup();
134                SimpleBlas.posv('U', A.dup(), X);
135                return X;
136        }
137
138  /** Computes the Least Squares solution for over or underdetermined
139   * linear equations A*X = B
140   *
141   * In the overdetermined case, when m > n, that is, there are more equations than
142   * variables, it computes the least squares solution of X -> ||A*X - B ||_2.
143   *
144   * In the underdetermined case, when m < n (less equations than variables), there are infinitely
145   * many solutions and it computes the minimum norm solution.
146   *
147   * @param A an (m,n) matrix
148   * @param B a (m,k) matrix
149   * @return either the minimum norm or least squares solution.
150   */
151  public static FloatMatrix solveLeastSquares(FloatMatrix A, FloatMatrix B) {
152    if (B.rows < A.columns) {
153      FloatMatrix X = FloatMatrix.concatVertically(B, new FloatMatrix(A.columns - B.rows, B.columns));
154      SimpleBlas.gelsd(A.dup(), X);
155      return X;
156    } else {
157      FloatMatrix X = B.dup();
158      SimpleBlas.gelsd(A.dup(), X);
159      return X.getRange(0, A.columns, 0, B.columns);
160    }
161  }
162
163  /**
164   * Computes the pseudo-inverse.
165   *
166   * Note, this function uses the solveLeastSquares and might produce different numerical
167   * solutions for the underdetermined case than matlab.
168   *
169   * @param A rectangular matrix
170   * @return matrix P such that A*P*A = A and P*A*P = P.
171   */
172  public static FloatMatrix pinv(FloatMatrix A) {
173    return solveLeastSquares(A, FloatMatrix.eye(A.rows));
174  }
175
176//END
177}