← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/linear/TriDiagonalTransformer.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/linear/TriDiagonalTransformer.java
 
1
/*
 
2
 * Licensed to the Apache Software Foundation (ASF) under one or more
 
3
 * contributor license agreements.  See the NOTICE file distributed with
 
4
 * this work for additional information regarding copyright ownership.
 
5
 * The ASF licenses this file to You under the Apache License, Version 2.0
 
6
 * (the "License"); you may not use this file except in compliance with
 
7
 * the License.  You may obtain a copy 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,
 
13
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 
14
 * See the License for the specific language governing permissions and
 
15
 * limitations under the License.
 
16
 */
 
17
 
 
18
package org.apache.commons.math.linear;
 
19
 
 
20
import java.util.Arrays;
 
21
 
 
22
 
 
23
/**
 
24
 * Class transforming a symmetrical matrix to tridiagonal shape.
 
25
 * <p>A symmetrical m &times; m matrix A can be written as the product of three matrices:
 
26
 * A = Q &times; T &times; Q<sup>T</sup> with Q an orthogonal matrix and T a symmetrical
 
27
 * tridiagonal matrix. Both Q and T are m &times; m matrices.</p>
 
28
 * <p>This implementation only uses the upper part of the matrix, the part below the
 
29
 * diagonal is not accessed at all.</p>
 
30
 * <p>Transformation to tridiagonal shape is often not a goal by itself, but it is
 
31
 * an intermediate step in more general decomposition algorithms like {@link
 
32
 * EigenDecomposition eigen decomposition}. This class is therefore intended for internal
 
33
 * use by the library and is not public. As a consequence of this explicitly limited scope,
 
34
 * many methods directly returns references to internal arrays, not copies.</p>
 
35
 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
 
36
 * @since 2.0
 
37
 */
 
38
class TriDiagonalTransformer {
 
39
 
 
40
    /** Householder vectors. */
 
41
    private final double householderVectors[][];
 
42
 
 
43
    /** Main diagonal. */
 
44
    private final double[] main;
 
45
 
 
46
    /** Secondary diagonal. */
 
47
    private final double[] secondary;
 
48
 
 
49
    /** Cached value of Q. */
 
50
    private RealMatrix cachedQ;
 
51
 
 
52
    /** Cached value of Qt. */
 
53
    private RealMatrix cachedQt;
 
54
 
 
55
    /** Cached value of T. */
 
56
    private RealMatrix cachedT;
 
57
 
 
58
    /**
 
59
     * Build the transformation to tridiagonal shape of a symmetrical matrix.
 
60
     * <p>The specified matrix is assumed to be symmetrical without any check.
 
61
     * Only the upper triangular part of the matrix is used.</p>
 
62
     * @param matrix the symmetrical matrix to transform.
 
63
     * @exception InvalidMatrixException if matrix is not square
 
64
     */
 
65
    public TriDiagonalTransformer(RealMatrix matrix)
 
66
        throws InvalidMatrixException {
 
67
        if (!matrix.isSquare()) {
 
68
            throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
 
69
        }
 
70
 
 
71
        final int m = matrix.getRowDimension();
 
72
        householderVectors = matrix.getData();
 
73
        main      = new double[m];
 
74
        secondary = new double[m - 1];
 
75
        cachedQ   = null;
 
76
        cachedQt  = null;
 
77
        cachedT   = null;
 
78
 
 
79
        // transform matrix
 
80
        transform();
 
81
 
 
82
    }
 
83
 
 
84
    /**
 
85
     * Returns the matrix Q of the transform. 
 
86
     * <p>Q is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
 
87
     * @return the Q matrix
 
88
     */
 
89
    public RealMatrix getQ() {
 
90
        if (cachedQ == null) {
 
91
            cachedQ = getQT().transpose();
 
92
        }
 
93
        return cachedQ;
 
94
    }
 
95
 
 
96
    /**
 
97
     * Returns the transpose of the matrix Q of the transform. 
 
98
     * <p>Q is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
 
99
     * @return the Q matrix
 
100
     */
 
101
    public RealMatrix getQT() {
 
102
 
 
103
        if (cachedQt == null) {
 
104
 
 
105
            final int m = householderVectors.length;
 
106
            cachedQt = MatrixUtils.createRealMatrix(m, m);
 
107
 
 
108
            // build up first part of the matrix by applying Householder transforms
 
109
            for (int k = m - 1; k >= 1; --k) {
 
110
                final double[] hK = householderVectors[k - 1];
 
111
                final double inv = 1.0 / (secondary[k - 1] * hK[k]);
 
112
                cachedQt.setEntry(k, k, 1);
 
113
                if (hK[k] != 0.0) {
 
114
                    double beta = 1.0 / secondary[k - 1];
 
115
                    cachedQt.setEntry(k, k, 1 + beta * hK[k]);
 
116
                    for (int i = k + 1; i < m; ++i) {
 
117
                        cachedQt.setEntry(k, i, beta * hK[i]);
 
118
                    }
 
119
                    for (int j = k + 1; j < m; ++j) {
 
120
                        beta = 0;
 
121
                        for (int i = k + 1; i < m; ++i) {
 
122
                            beta += cachedQt.getEntry(j, i) * hK[i];
 
123
                        }
 
124
                        beta *= inv;
 
125
                        cachedQt.setEntry(j, k, beta * hK[k]);
 
126
                        for (int i = k + 1; i < m; ++i) {
 
127
                            cachedQt.addToEntry(j, i, beta * hK[i]);
 
128
                        }
 
129
                    }
 
130
                }
 
131
            }
 
132
            cachedQt.setEntry(0, 0, 1);
 
133
 
 
134
        }
 
135
 
 
136
        // return the cached matrix
 
137
        return cachedQt;
 
138
 
 
139
    }
 
140
 
 
141
    /**
 
142
     * Returns the tridiagonal matrix T of the transform. 
 
143
     * @return the T matrix
 
144
     */
 
145
    public RealMatrix getT() {
 
146
 
 
147
        if (cachedT == null) {
 
148
 
 
149
            final int m = main.length;
 
150
            cachedT = MatrixUtils.createRealMatrix(m, m);
 
151
            for (int i = 0; i < m; ++i) {
 
152
                cachedT.setEntry(i, i, main[i]);
 
153
                if (i > 0) {
 
154
                    cachedT.setEntry(i, i - 1, secondary[i - 1]);
 
155
                }
 
156
                if (i < main.length - 1) {
 
157
                    cachedT.setEntry(i, i + 1, secondary[i]);
 
158
                }
 
159
            }
 
160
 
 
161
        }
 
162
 
 
163
        // return the cached matrix
 
164
        return cachedT;
 
165
 
 
166
    }
 
167
 
 
168
    /**
 
169
     * Get the Householder vectors of the transform.
 
170
     * <p>Note that since this class is only intended for internal use,
 
171
     * it returns directly a reference to its internal arrays, not a copy.</p>
 
172
     * @return the main diagonal elements of the B matrix
 
173
     */
 
174
    double[][] getHouseholderVectorsRef() {
 
175
        return householderVectors;
 
176
    }
 
177
 
 
178
    /**
 
179
     * Get the main diagonal elements of the matrix T of the transform.
 
180
     * <p>Note that since this class is only intended for internal use,
 
181
     * it returns directly a reference to its internal arrays, not a copy.</p>
 
182
     * @return the main diagonal elements of the T matrix
 
183
     */
 
184
    double[] getMainDiagonalRef() {
 
185
        return main;
 
186
    }
 
187
 
 
188
    /**
 
189
     * Get the secondary diagonal elements of the matrix T of the transform.
 
190
     * <p>Note that since this class is only intended for internal use,
 
191
     * it returns directly a reference to its internal arrays, not a copy.</p>
 
192
     * @return the secondary diagonal elements of the T matrix
 
193
     */
 
194
    double[] getSecondaryDiagonalRef() {
 
195
        return secondary;
 
196
    }
 
197
 
 
198
    /**
 
199
     * Transform original matrix to tridiagonal form.
 
200
     * <p>Transformation is done using Householder transforms.</p>
 
201
     */
 
202
    private void transform() {
 
203
 
 
204
        final int m = householderVectors.length;
 
205
        final double[] z = new double[m];
 
206
        for (int k = 0; k < m - 1; k++) {
 
207
 
 
208
            //zero-out a row and a column simultaneously
 
209
            final double[] hK = householderVectors[k];
 
210
            main[k] = hK[k];
 
211
            double xNormSqr = 0;
 
212
            for (int j = k + 1; j < m; ++j) {
 
213
                final double c = hK[j];
 
214
                xNormSqr += c * c;
 
215
            }
 
216
            final double a = (hK[k + 1] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
 
217
            secondary[k] = a;
 
218
            if (a != 0.0) {
 
219
                // apply Householder transform from left and right simultaneously
 
220
 
 
221
                hK[k + 1] -= a;
 
222
                final double beta = -1 / (a * hK[k + 1]);
 
223
 
 
224
                // compute a = beta A v, where v is the Householder vector
 
225
                // this loop is written in such a way
 
226
                //   1) only the upper triangular part of the matrix is accessed
 
227
                //   2) access is cache-friendly for a matrix stored in rows
 
228
                Arrays.fill(z, k + 1, m, 0);
 
229
                for (int i = k + 1; i < m; ++i) {
 
230
                    final double[] hI = householderVectors[i];
 
231
                    final double hKI = hK[i];
 
232
                    double zI = hI[i] * hKI;
 
233
                    for (int j = i + 1; j < m; ++j) {
 
234
                        final double hIJ = hI[j];
 
235
                        zI   += hIJ * hK[j];
 
236
                        z[j] += hIJ * hKI;
 
237
                    }
 
238
                    z[i] = beta * (z[i] + zI);
 
239
                }
 
240
 
 
241
                // compute gamma = beta vT z / 2
 
242
                double gamma = 0;
 
243
                for (int i = k + 1; i < m; ++i) {
 
244
                    gamma += z[i] * hK[i];
 
245
                }
 
246
                gamma *= beta / 2;
 
247
 
 
248
                // compute z = z - gamma v
 
249
                for (int i = k + 1; i < m; ++i) {
 
250
                    z[i] -= gamma * hK[i];
 
251
                }
 
252
 
 
253
                // update matrix: A = A - v zT - z vT
 
254
                // only the upper triangular part of the matrix is updated
 
255
                for (int i = k + 1; i < m; ++i) {
 
256
                    final double[] hI = householderVectors[i];
 
257
                    for (int j = i; j < m; ++j) {
 
258
                        hI[j] -= hK[i] * z[j] + z[i] * hK[j];
 
259
                    }
 
260
                }
 
261
 
 
262
            }
 
263
 
 
264
        }
 
265
        main[m - 1] = householderVectors[m - 1][m - 1];
 
266
    }
 
267
 
 
268
}