/* Cholesky.java copyright claim: Copyright (C) 1996 by Steve Verrill. This class is free software; you can redistribute it and/or modify it under the terms of version 2 of the GNU General Public License as published by the Free Software Foundation. This class is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with these programs in the file COPYING; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. The author's mail address is: Steve Verrill USDA Forest Products Laboratory 1 Gifford Pinchot Drive Madison, Wisconsin 53705 The author's e-mail address is: steve@swst.org *********************************************************************** DISCLAIMER OF WARRANTIES: THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND. THE AUTHOR DOES NOT WARRANT, GUARANTEE OR MAKE ANY REPRESENTATIONS REGARDING THE SOFTWARE OR DOCUMENTATION IN TERMS OF THEIR CORRECTNESS, RELIABILITY, CURRENTNESS, OR OTHERWISE. THE ENTIRE RISK AS TO THE RESULTS AND PERFORMANCE OF THE SOFTWARE IS ASSUMED BY YOU. IN NO CASE WILL ANY PARTY INVOLVED WITH THE CREATION OR DISTRIBUTION OF THE SOFTWARE BE LIABLE FOR ANY DAMAGE THAT MAY RESULT FROM THE USE OF THIS SOFTWARE. Sorry about that. *********************************************************************** History: Date Author Changes 11/30/96 Steve Verrill Created */ package linear_algebra; import linear_algebra.*; /** * **This class contains: *
*
* *- a method that obtains the Cholesky *factorization RR´, where R is a lower triangular matrix, *of a symmetric positive definite matrix A. *
- a method to invert a symmetric positive *definite matrix. *
- a method to solve Ax = b where A is a symmetric *positive definite matrix. *
*This class was written by a statistician rather than *a numerical analyst. I have tried to check the code carefully, *but it may still contain bugs. Further, its stability and *efficiency do not meet the standards of high quality *numerical analysis software. When public domain Java numerical *analysis routines become available from numerical analysts *(e.g., the people who produce LAPACK), then THE CODE PRODUCED *BY THE NUMERICAL ANALYSTS SHOULD BE USED. * *
*Meanwhile, if you have suggestions for improving this *code, please contact Steve Verrill at *. * *@author Steve Verrill *@version .5 --- November 30, 1996 * */ public class Cholesky extends Object { /** * *
*The method proceeds by first factoring A as RR´ where R is a lower *triangular matrix. Thus Ax = b is equivalent to R(R´x) = b. Then the *method performs two additional operations. First it solves Ry = b for y. *Then it solves R´x = y for x. It stores x in b. * *@param a[ ][ ] On entrance, if the boolean variable factored is false, *the lower triangle of a[ ][ ] should contain the lower triangle of A. *If factored is true, the lower triangle of a[ ][ ] should *contain a lower triangular matrix R such that RR´ = A. *@param b On entrance b must contain the known b of Ax = b. On exit *it contains the solution x to Ax = b. *@param y A work vector of order at least n. *@param n The order of A and b. *@param factored On entrance, factored should be set to true if A *already has been factored, false *if A has not yet been factored. *@throws NotPosDefException if A cannot be factored as RR´ for a *full rank lower triangular matrix R. * */ public void solvePosDef (double a[][], double b[], double y[], int n, boolean factored) throws NotPosDefException { Cholesky cholesky = new Cholesky(); Triangular triangular = new Triangular(); int i,j; if (factored == false) { /* A has not yet been factored. Factor it. The lower triangle of a[][] will contain R where RR' = A. */ try { cholesky.factorPosDef(a,n); } catch (NotPosDefException npd) { System.out.print("\nCholesky.solvePosDef error:" + " The matrix was not positive definite\n" + "so it was not possible" + "to use the Cholesky decomposition\n" + " to solve Ax = b.\n\n"); throw npd; } } /* Solve Ry = b. */ try { triangular.solveLower(a, y, b, n); } catch (NotFullRankException nfr) { System.out.print("\nTriangular.solveLower error:" + " The lower triangular factor of A\n" + "was not of full rank.\n\n"); throw new NotPosDefException(); } /* Fill the upper triangle of a[][] with R' where RR' = A. */ for (i = 0; i < n; i++) { for (j = i; j < n; j++) { a[i][j] = a[j][i]; } } /* Solve R'x = y. Put x in the b vector. */ try { triangular.solveUpper(a, b, y, n); } catch (NotFullRankException nfr) { System.out.print("\nTriangular.solveUpper error:" + " The upper triangular factor of A\n" + "was not of full rank.\n\n"); throw new NotPosDefException(); } return; } /** * *
*This method obtains the inverse of an n by n
*symmetric positive definite matrix A.
*