KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > org > apache > commons > math > linear > CholeskySolver


1 /*
2  *
3  * Copyright (c) 2003-2004 The Apache Software Foundation. All rights reserved.
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License"); you may not
6  * use this file except in compliance with the License. You may obtain a copy
7  * 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, WITHOUT
13  * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
14  * License for the specific language governing permissions and limitations
15  * under the License.
16  *
17  */
18
19 package org.apache.commons.math.linear;
20
21
22 /**
23  * Solves a linear equitation with symmetrical, positiv definit
24  * coefficient matrix by Cholesky decomposition.
25  * <p>
26  * For every symmetric, positiv definit matrix <code>M</code> there is a
27  * lower triangular matrix <code>L</code> so that <code>L*L^T=M</code>.
28  * <code>L</code> is called the <i>Cholesky decomposition</i> of <code>M</code>.
29  * For any constant vector <code>c</code> it can be used to solve
30  * the linear equitation <code>M*x=L*(L^T*x)=c</code>.<br>
31  * Compared to the LU-decompoistion the Cholesky methods requires only half
32  * the number of operations.
33  * <p>
34  * @author Stefan Koeberle, 11/2003
35  */
36 public class CholeskySolver {
37     
38     private double numericalZero = 10E-12;
39     
40     /** The lower triangular matrix */
41     private RealMatrixImpl decompMatrix;
42     
43    
44     /**
45      * Creates a new instance of CholeskySolver
46      */
47     public CholeskySolver() {
48     }//constructor CholeskySolver
49
50     
51     /**
52      * Every double <code>d</code> satisfying
53      * <code>java.lang.Math.abs(d) <= numericalZero</code>
54      * is considered equal to <code>0.0d.</code>
55      */
56     public void setNumericalZero(double numericalZero) {
57         this.numericalZero = numericalZero;
58     }//setNumericalZero
59
60     /**
61      * See <code>setNumericalZero</code>
62      */
63     public double getNumericalZero() {
64         return numericalZero;
65     }//getNumericalZero
66
67     
68     /**
69      * Calculates the Cholesky-decomposition of the symmetrical, positiv definit
70      * matrix <code>M</code>.
71      * <p>
72      * The decomposition matrix is internally stored.
73      * <p>
74      * @throws IllegalArgumentException if <code>M</code> ist not square or
75      * not positiv definit
76      */
77     public void decompose(RealMatrix m)
78     throws IllegalArgumentException JavaDoc {
79        
80        decompMatrix = null;
81        double[][] mval = m.getData();
82        int numRows = m.getRowDimension();
83        int numCols = m.getColumnDimension();
84        if (numRows != numCols)
85            throw new IllegalArgumentException JavaDoc("matrix is not square");
86        double[][] decomp = new double[numRows][numCols];
87        double sum;
88        
89        //for all columns
90 for (int col=0; col<numCols; col++) {
91           
92            //diagonal element
93 sum = mval[col][col];
94            for (int k=0; k<col; k++)
95                sum = sum - decomp[col][k]*decomp[col][k];
96            if (sum <= numericalZero) {
97                throw new IllegalArgumentException JavaDoc(
98                              "Matrix is not positiv definit");
99            }
100            decomp[col][col] += Math.sqrt(sum);
101            
102            //column below diagonal
103 for (int row=col+1; row<numRows; row++) {
104                sum = mval[row][col];
105                for (int k=0; k<col; k++)
106                    sum = sum - decomp[col][k]*decomp[row][k];
107                decomp[row][col] = sum/decomp[col][col];
108            }//for
109
110        }//for all columns
111
112        decompMatrix = new RealMatrixImpl(decomp);
113        
114     }//decompose
115
116     
117     /**
118      * Returns the last calculated decomposition matrix.
119      * <p>
120      * Caution: Every call of this Method will return the same object.
121      * Decomposing another matrix will generate a new one.
122      */
123     public RealMatrixImpl getDecomposition() {
124         return decompMatrix;
125     }//getDecomposition
126
127     
128     /**
129      * Returns the solution for a linear system with constant vector <code>c</code>.
130      * <p>
131      * This method solves a linear system <code>M*x=c</code> for a symmetrical,
132      * positiv definit coefficient matrix <code>M</code>. Before using this
133      * method the matrix <code>M</code> must have been decomposed.
134      * <p>
135      * @throws IllegalStateException if this methode is called before
136      * a matrix was decomposed
137      * @throws IllegalArgumentException if the dimension of <code>c</code> doesn't
138      * match the row dimension of <code>M</code>
139      */
140     public double[] solve(double[] c)
141     throws IllegalStateException JavaDoc, IllegalArgumentException JavaDoc {
142       
143         if (decompMatrix == null) {
144             throw new IllegalStateException JavaDoc("no decomposed matrix available");
145         }//if
146 if (decompMatrix.getColumnDimension() != c.length)
147            throw new IllegalArgumentException JavaDoc("matrix dimension mismatch");
148        
149         double[][] decomp = decompMatrix.getData();
150         double[] x = new double[decomp.length];
151         double sum;
152         
153         //forward elimination
154 for (int i=0; i<x.length; i++) {
155             sum = c[i];
156             for (int k=0; k<i; k++)
157                 sum = sum - decomp[i][k]*x[k];
158             x[i] = sum / decomp[i][i];
159         }//forward elimination
160
161         //backward elimination
162 for (int i=x.length-1; i>=0; i--) {
163             sum = x[i];
164             for (int k=i+1; k<x.length; k++)
165                 sum = sum - decomp[k][i]*x[k];
166             x[i] = sum / decomp[i][i];
167         }//backward elimination
168
169         return x;
170     }//solve
171
172     
173     /**
174      * Returns the solution for a linear system with a symmetrical,
175      * positiv definit coefficient matrix <code>M</code> and
176      * constant vector <code>c</code>.
177      * <p>
178      * As a side effect, the Cholesky-decomposition <code>L*L^T=M</code> is
179      * calculated and internally stored.
180      * <p>
181      * This is a convenience method for <code><pre>
182      * solver.decompose(m);
183      * solver.solve(c);
184      * </pre></code>
185      * @throws IllegalArgumentException if M ist not square, not positive definit
186      * or the dimensions of <code>M</code> and
187      * <code>c</code> don't match.
188      */
189     public double[] solve(RealMatrix m, double[] c)
190     throws IllegalArgumentException JavaDoc {
191         decompose(m);
192         return solve(c);
193     }//solve
194
195     
196     /**
197      * Returns the determinant of the a matrix <code>M</code>.
198      * <p>
199      * Before using this method the matrix <code>M</code> must
200      * have been decomposed.
201      * <p>
202      * @throws IllegalStateException if this method is called before
203      * a matrix was decomposed
204      */
205     public double getDeterminant() {
206         
207         if (decompMatrix == null) {
208             throw new IllegalStateException JavaDoc("no decomposed matrix available");
209         }//if
210
211         double[][] data = decompMatrix.getData();
212         double res = 1.0d;
213         for (int i=0; i<data.length; i++) {
214             res *= data[i][i];
215         }//for
216 res = res*res;
217         
218         return res;
219     }//getDeterminant
220
221 }//class CholeskySolver
222
Popular Tags