2 * $RCSfile: GMatrix.java,v $
4 * Copyright 1997-2008 Sun Microsystems, Inc. All Rights Reserved.
5 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
7 * This code is free software; you can redistribute it and/or modify it
8 * under the terms of the GNU General Public License version 2 only, as
9 * published by the Free Software Foundation. Sun designates this
10 * particular file as subject to the "Classpath" exception as provided
11 * by Sun in the LICENSE file that accompanied this code.
13 * This code is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 * version 2 for more details (a copy is included in the LICENSE file that
17 * accompanied this code).
19 * You should have received a copy of the GNU General Public License version
20 * 2 along with this work; if not, write to the Free Software Foundation,
21 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
23 * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,
24 * CA 95054 USA or visit www.sun.com if you need additional information or
28 * $Date: 2008/02/28 20:18:50 $
32 package javax.vecmath;
34 import java.lang.Math;
37 * A double precision, general, dynamically-resizable,
38 * two-dimensional matrix class. Row and column numbering begins with
39 * zero. The representation is row major.
42 public class GMatrix implements java.io.Serializable, Cloneable {
44 // Compatible with 1.1
45 static final long serialVersionUID = 2777097312029690941L;
46 private static final boolean debug = false;
51 // double dereference is slow
54 private static final double EPS = 1.0E-10;
57 * Constructs an nRow by NCol identity matrix.
58 * Note that because row and column numbering begins with
59 * zero, nRow and nCol will be one larger than the maximum
60 * possible matrix index values.
61 * @param nRow number of rows in this matrix.
62 * @param nCol number of columns in this matrix.
64 public GMatrix(int nRow, int nCol)
66 values = new double[nRow][nCol];
71 for (i = 0; i < nRow; i++) {
72 for (j = 0; j < nCol; j++) {
83 for (i = 0; i < l; i++) {
89 * Constructs an nRow by nCol matrix initialized to the values
90 * in the matrix array. The array values are copied in one row at
91 * a time in row major fashion. The array should be at least
92 * nRow*nCol in length.
93 * Note that because row and column numbering begins with
94 * zero, nRow and nCol will be one larger than the maximum
95 * possible matrix index values.
96 * @param nRow number of rows in this matrix.
97 * @param nCol number of columns in this matrix.
98 * @param matrix a 1D array that specifies a matrix in row major fashion
100 public GMatrix(int nRow, int nCol, double[] matrix)
102 values = new double[nRow][nCol];
107 for (i = 0; i < nRow; i++) {
108 for (j = 0; j < nCol; j++) {
109 values[i][j] = matrix[i*nCol+j];
115 * Constructs a new GMatrix and copies the initial values
116 * from the parameter matrix.
117 * @param matrix the source of the initial values of the new GMatrix
119 public GMatrix(GMatrix matrix)
123 values = new double[nRow][nCol];
126 for (i = 0; i < nRow; i++) {
127 for (j = 0; j < nCol; j++) {
128 values[i][j] = matrix.values[i][j];
134 * Sets the value of this matrix to the result of multiplying itself
135 * with matrix m1 (this = this * m1).
136 * @param m1 the other matrix
138 public final void mul(GMatrix m1)
142 if (nCol != m1.nRow || nCol != m1.nCol)
143 throw new MismatchedSizeException
144 (VecMathI18N.getString("GMatrix0"));
146 double [][] tmp = new double[nRow][nCol];
148 for (i = 0; i < nRow; i++) {
149 for (j = 0; j < nCol; j++) {
151 for (k = 0; k < nCol; k++) {
152 tmp[i][j] += values[i][k]*m1.values[k][j];
161 * Sets the value of this matrix to the result of multiplying
162 * the two argument matrices together (this = m1 * m2).
163 * @param m1 the first matrix
164 * @param m2 the second matrix
166 public final void mul(GMatrix m1, GMatrix m2)
170 if (m1.nCol != m2.nRow || nRow != m1.nRow || nCol != m2.nCol)
171 throw new MismatchedSizeException
172 (VecMathI18N.getString("GMatrix1"));
174 double[][] tmp = new double[nRow][nCol];
176 for (i = 0; i < m1.nRow; i++) {
177 for (j = 0; j < m2.nCol; j++) {
179 for (k = 0; k < m1.nCol; k++) {
180 tmp[i][j] += m1.values[i][k]*m2.values[k][j];
189 * Computes the outer product of the two vectors; multiplies the
190 * the first vector by the transpose of the second vector and places
191 * the matrix result into this matrix. This matrix must be
192 * be as big or bigger than getSize(v1)xgetSize(v2).
193 * @param v1 the first vector, treated as a row vector
194 * @param v2 the second vector, treated as a column vector
196 public final void mul(GVector v1, GVector v2)
200 if (nRow < v1.getSize())
201 throw new MismatchedSizeException
202 (VecMathI18N.getString("GMatrix2"));
204 if (nCol < v2.getSize())
205 throw new MismatchedSizeException
206 (VecMathI18N.getString("GMatrix3"));
208 for (i = 0; i < v1.getSize(); i++ ) {
209 for (j = 0; j < v2.getSize(); j++ ) {
210 values[i][j] = v1.values[i]*v2.values[j];
216 * Sets the value of this matrix to sum of itself and matrix m1.
217 * @param m1 the other matrix
219 public final void add(GMatrix m1)
224 throw new MismatchedSizeException
225 (VecMathI18N.getString("GMatrix4"));
228 throw new MismatchedSizeException
229 (VecMathI18N.getString("GMatrix5"));
231 for (i = 0; i < nRow; i++) {
232 for (j = 0; j < nCol; j++) {
233 values[i][j] = values[i][j] + m1.values[i][j];
239 * Sets the value of this matrix to the matrix sum of matrices m1 and m2.
240 * @param m1 the first matrix
241 * @param m2 the second matrix
243 public final void add(GMatrix m1, GMatrix m2)
247 if (m2.nRow != m1.nRow)
248 throw new MismatchedSizeException
249 (VecMathI18N.getString("GMatrix6"));
251 if (m2.nCol != m1.nCol)
252 throw new MismatchedSizeException
253 (VecMathI18N.getString("GMatrix7"));
255 if (nCol != m1.nCol || nRow != m1.nRow)
256 throw new MismatchedSizeException
257 (VecMathI18N.getString("GMatrix8"));
259 for (i = 0; i < nRow; i++) {
260 for (j = 0; j < nCol; j++) {
261 values[i][j] = m1.values[i][j] + m2.values[i][j];
267 * Sets the value of this matrix to the matrix difference of itself
268 * and matrix m1 (this = this - m1).
269 * @param m1 the other matrix
271 public final void sub(GMatrix m1)
275 throw new MismatchedSizeException
276 (VecMathI18N.getString("GMatrix9"));
279 throw new MismatchedSizeException
280 (VecMathI18N.getString("GMatrix28"));
282 for (i = 0; i < nRow; i++) {
283 for (j = 0; j < nCol; j++) {
284 values[i][j] = values[i][j] - m1.values[i][j];
290 * Sets the value of this matrix to the matrix difference
291 * of matrices m1 and m2 (this = m1 - m2).
292 * @param m1 the first matrix
293 * @param m2 the second matrix
295 public final void sub(GMatrix m1, GMatrix m2)
298 if (m2.nRow != m1.nRow)
299 throw new MismatchedSizeException
300 (VecMathI18N.getString("GMatrix10"));
302 if (m2.nCol != m1.nCol)
303 throw new MismatchedSizeException
304 (VecMathI18N.getString("GMatrix11"));
306 if (nRow != m1.nRow || nCol != m1.nCol)
307 throw new MismatchedSizeException
308 (VecMathI18N.getString("GMatrix12"));
310 for (i = 0; i < nRow; i++) {
311 for (j = 0; j < nCol; j++) {
312 values[i][j] = m1.values[i][j] - m2.values[i][j];
318 * Negates the value of this matrix: this = -this.
320 public final void negate()
323 for (i = 0; i < nRow; i++) {
324 for (j = 0;j < nCol; j++) {
325 values[i][j] = -values[i][j];
331 * Sets the value of this matrix equal to the negation of
332 * of the GMatrix parameter.
333 * @param m1 The source matrix
335 public final void negate(GMatrix m1)
338 if (nRow != m1.nRow || nCol != m1.nCol)
339 throw new MismatchedSizeException
340 (VecMathI18N.getString("GMatrix13"));
342 for (i = 0; i < nRow; i++) {
343 for (j = 0; j < nCol; j++) {
344 values[i][j] = -m1.values[i][j];
350 * Sets this GMatrix to the identity matrix.
352 public final void setIdentity()
355 for (i = 0; i < nRow; i++) {
356 for (j = 0; j < nCol; j++) {
367 for (i = 0; i < l; i++) {
373 * Sets all the values in this matrix to zero.
375 public final void setZero()
378 for (i = 0; i < nRow; i++) {
379 for (j = 0; j < nCol; j++) {
386 * Subtracts this matrix from the identity matrix and puts the values
387 * back into this (this = I - this).
389 public final void identityMinus()
393 for(i = 0; i < nRow; i++) {
394 for(j = 0; j < nCol; j++) {
395 values[i][j] = -values[i][j];
405 for(i = 0; i < l; i++) {
412 * Inverts this matrix in place.
414 public final void invert()
420 * Inverts matrix m1 and places the new values into this matrix. Matrix
421 * m1 is not modified.
422 * @param m1 the matrix to be inverted
424 public final void invert(GMatrix m1)
430 * Copies a sub-matrix derived from this matrix into the target matrix.
431 * The upper left of the sub-matrix is located at (rowSource, colSource);
432 * the lower right of the sub-matrix is located at
433 * (lastRowSource,lastColSource). The sub-matrix is copied into the
434 * the target matrix starting at (rowDest, colDest).
435 * @param rowSource the top-most row of the sub-matrix
436 * @param colSource the left-most column of the sub-matrix
437 * @param numRow the number of rows in the sub-matrix
438 * @param numCol the number of columns in the sub-matrix
439 * @param rowDest the top-most row of the position of the copied
440 * sub-matrix within the target matrix
441 * @param colDest the left-most column of the position of the copied
442 * sub-matrix within the target matrix
443 * @param target the matrix into which the sub-matrix will be copied
445 public final void copySubMatrix(int rowSource, int colSource,
446 int numRow, int numCol, int rowDest,
447 int colDest, GMatrix target)
451 if (this != target) {
452 for (i = 0; i < numRow; i++) {
453 for (j = 0; j < numCol; j++) {
454 target.values[rowDest+i][colDest+j] =
455 values[rowSource+i][colSource+j];
459 double[][] tmp = new double[numRow][numCol];
460 for (i = 0; i < numRow; i++) {
461 for (j = 0; j < numCol; j++) {
462 tmp[i][j] = values[rowSource+i][colSource+j];
465 for (i = 0; i < numRow; i++) {
466 for (j = 0; j < numCol; j++) {
467 target.values[rowDest+i][colDest+j] = tmp[i][j];
474 * Changes the size of this matrix dynamically. If the size is increased
475 * no data values will be lost. If the size is decreased, only those data
476 * values whose matrix positions were eliminated will be lost.
477 * @param nRow number of desired rows in this matrix
478 * @param nCol number of desired columns in this matrix
480 public final void setSize(int nRow, int nCol)
482 double[][] tmp = new double[nRow][nCol];
483 int i, j, maxRow, maxCol;
485 if (this.nRow < nRow)
490 if (this.nCol < nCol)
495 for (i = 0; i < maxRow; i++) {
496 for (j = 0; j < maxCol; j++) {
497 tmp[i][j] = values[i][j];
508 * Sets the value of this matrix to the values found in the array parameter.
509 * The values are copied in one row at a time, in row major
510 * fashion. The array should be at least equal in length to
511 * the number of matrix rows times the number of matrix columns
513 * @param matrix the row major source array
515 public final void set(double[] matrix)
519 for (i = 0; i < nRow; i++) {
520 for (j = 0; j < nCol; j++) {
521 values[i][j] = matrix[nCol*i+j];
527 * Sets the value of this matrix to that of the Matrix3f provided.
528 * @param m1 the matrix
530 public final void set(Matrix3f m1)
534 if (nCol < 3 || nRow < 3) { // expand matrix if too small
537 values = new double[nRow][nCol];
540 values[0][0] = m1.m00;
541 values[0][1] = m1.m01;
542 values[0][2] = m1.m02;
544 values[1][0] = m1.m10;
545 values[1][1] = m1.m11;
546 values[1][2] = m1.m12;
548 values[2][0] = m1.m20;
549 values[2][1] = m1.m21;
550 values[2][2] = m1.m22;
552 for (i = 3; i < nRow; i++) { // pad rest or matrix with zeros
553 for (j = 3; j < nCol; j++) {
560 * Sets the value of this matrix to that of the Matrix3d provided.
561 * @param m1 the matrix
563 public final void set(Matrix3d m1)
565 if (nRow < 3 || nCol < 3) {
566 values = new double[3][3];
571 values[0][0] = m1.m00;
572 values[0][1] = m1.m01;
573 values[0][2] = m1.m02;
575 values[1][0] = m1.m10;
576 values[1][1] = m1.m11;
577 values[1][2] = m1.m12;
579 values[2][0] = m1.m20;
580 values[2][1] = m1.m21;
581 values[2][2] = m1.m22;
583 for (int i = 3; i < nRow; i++) { // pad rest or matrix with zeros
584 for(int j = 3; j < nCol; j++) {
592 * Sets the value of this matrix to that of the Matrix4f provided.
593 * @param m1 the matrix
595 public final void set(Matrix4f m1)
597 if (nRow < 4 || nCol < 4) {
598 values = new double[4][4];
603 values[0][0] = m1.m00;
604 values[0][1] = m1.m01;
605 values[0][2] = m1.m02;
606 values[0][3] = m1.m03;
608 values[1][0] = m1.m10;
609 values[1][1] = m1.m11;
610 values[1][2] = m1.m12;
611 values[1][3] = m1.m13;
613 values[2][0] = m1.m20;
614 values[2][1] = m1.m21;
615 values[2][2] = m1.m22;
616 values[2][3] = m1.m23;
618 values[3][0] = m1.m30;
619 values[3][1] = m1.m31;
620 values[3][2] = m1.m32;
621 values[3][3] = m1.m33;
623 for (int i = 4 ; i < nRow; i++) { // pad rest or matrix with zeros
624 for (int j = 4; j < nCol; j++) {
631 * Sets the value of this matrix to that of the Matrix4d provided.
632 * @param m1 the matrix
634 public final void set(Matrix4d m1)
636 if (nRow < 4 || nCol < 4) {
637 values = new double[4][4];
642 values[0][0] = m1.m00;
643 values[0][1] = m1.m01;
644 values[0][2] = m1.m02;
645 values[0][3] = m1.m03;
647 values[1][0] = m1.m10;
648 values[1][1] = m1.m11;
649 values[1][2] = m1.m12;
650 values[1][3] = m1.m13;
652 values[2][0] = m1.m20;
653 values[2][1] = m1.m21;
654 values[2][2] = m1.m22;
655 values[2][3] = m1.m23;
657 values[3][0] = m1.m30;
658 values[3][1] = m1.m31;
659 values[3][2] = m1.m32;
660 values[3][3] = m1.m33;
662 for (int i = 4; i < nRow; i++) { // pad rest or matrix with zeros
663 for (int j = 4; j < nCol; j++) {
670 * Sets the value of this matrix to the values found in matrix m1.
671 * @param m1 the source matrix
673 public final void set(GMatrix m1)
677 if (nRow < m1.nRow || nCol < m1.nCol) {
680 values = new double[nRow][nCol];
683 for (i = 0; i < Math.min(nRow, m1.nRow); i++) {
684 for (j = 0; j < Math.min(nCol, m1.nCol); j++) {
685 values[i][j] = m1.values[i][j];
689 for (i = m1.nRow; i < nRow; i++) { // pad rest or matrix with zeros
690 for (j = m1.nCol; j < nCol; j++) {
697 * Returns the number of rows in this matrix.
698 * @return number of rows in this matrix
700 public final int getNumRow()
706 * Returns the number of colmuns in this matrix.
707 * @return number of columns in this matrix
709 public final int getNumCol()
715 * Retrieves the value at the specified row and column of this matrix.
716 * @param row the row number to be retrieved (zero indexed)
717 * @param column the column number to be retrieved (zero indexed)
718 * @return the value at the indexed element
720 public final double getElement(int row, int column)
722 return(values[row][column]);
727 * Modifies the value at the specified row and column of this matrix.
728 * @param row the row number to be modified (zero indexed)
729 * @param column the column number to be modified (zero indexed)
730 * @param value the new matrix element value
732 public final void setElement(int row, int column, double value)
734 values[row][column] = value;
738 * Places the values of the specified row into the array parameter.
739 * @param row the target row number
740 * @param array the array into which the row values will be placed
742 public final void getRow(int row, double[] array)
744 for (int i = 0; i < nCol; i++) {
745 array[i] = values[row][i];
750 * Places the values of the specified row into the vector parameter.
751 * @param row the target row number
752 * @param vector the vector into which the row values will be placed
754 public final void getRow(int row, GVector vector)
756 if (vector.getSize() < nCol)
757 vector.setSize(nCol);
759 for (int i = 0; i < nCol; i++) {
760 vector.values[i] = values[row][i];
765 * Places the values of the specified column into the array parameter.
766 * @param col the target column number
767 * @param array the array into which the column values will be placed
769 public final void getColumn(int col, double[] array)
771 for (int i = 0; i < nRow; i++) {
772 array[i] = values[i][col];
778 * Places the values of the specified column into the vector parameter.
779 * @param col the target column number
780 * @param vector the vector into which the column values will be placed
782 public final void getColumn(int col, GVector vector)
784 if (vector.getSize() < nRow)
785 vector.setSize(nRow);
787 for (int i = 0; i < nRow; i++) {
788 vector.values[i] = values[i][col];
793 * Places the values in the upper 3x3 of this GMatrix into
795 * @param m1 The matrix that will hold the new values
797 public final void get(Matrix3d m1)
799 if (nRow < 3 || nCol < 3) {
803 m1.m00 = values[0][0];
805 m1.m10 = values[1][0];
807 m1.m20= values[2][0];
813 m1.m01 = values[0][1];
815 m1.m11 = values[1][1];
817 m1.m21 = values[2][1];
823 m1.m02 = values[0][2];
825 m1.m12 = values[1][2];
827 m1.m22 = values[2][2];
835 m1.m00 = values[0][0];
836 m1.m01 = values[0][1];
837 m1.m02 = values[0][2];
839 m1.m10 = values[1][0];
840 m1.m11 = values[1][1];
841 m1.m12 = values[1][2];
843 m1.m20 = values[2][0];
844 m1.m21 = values[2][1];
845 m1.m22 = values[2][2];
850 * Places the values in the upper 3x3 of this GMatrix into
852 * @param m1 The matrix that will hold the new values
854 public final void get(Matrix3f m1)
857 if (nRow < 3 || nCol < 3) {
861 m1.m00 = (float)values[0][0];
863 m1.m10 = (float)values[1][0];
865 m1.m20 = (float)values[2][0];
871 m1.m01 = (float)values[0][1];
873 m1.m11 = (float)values[1][1];
875 m1.m21 = (float)values[2][1];
881 m1.m02 = (float)values[0][2];
883 m1.m12 = (float)values[1][2];
885 m1.m22 = (float)values[2][2];
893 m1.m00 = (float)values[0][0];
894 m1.m01 = (float)values[0][1];
895 m1.m02 = (float)values[0][2];
897 m1.m10 = (float)values[1][0];
898 m1.m11 = (float)values[1][1];
899 m1.m12 = (float)values[1][2];
901 m1.m20 = (float)values[2][0];
902 m1.m21 = (float)values[2][1];
903 m1.m22 = (float)values[2][2];
908 * Places the values in the upper 4x4 of this GMatrix into
910 * @param m1 The matrix that will hold the new values
912 public final void get(Matrix4d m1)
914 if (nRow < 4 || nCol < 4) {
918 m1.m00 = values[0][0];
920 m1.m10 = values[1][0];
922 m1.m20 = values[2][0];
924 m1.m30 = values[3][0];
931 m1.m01 = values[0][1];
933 m1.m11 = values[1][1];
935 m1.m21 = values[2][1];
937 m1.m31 = values[3][1];
944 m1.m02 = values[0][2];
946 m1.m12 = values[1][2];
948 m1.m22 = values[2][2];
950 m1.m32 = values[3][2];
957 m1.m03 = values[0][3];
959 m1.m13 = values[1][3];
961 m1.m23 = values[2][3];
963 m1.m33 = values[3][3];
973 m1.m00 = values[0][0];
974 m1.m01 = values[0][1];
975 m1.m02 = values[0][2];
976 m1.m03 = values[0][3];
978 m1.m10 = values[1][0];
979 m1.m11 = values[1][1];
980 m1.m12 = values[1][2];
981 m1.m13 = values[1][3];
983 m1.m20 = values[2][0];
984 m1.m21 = values[2][1];
985 m1.m22 = values[2][2];
986 m1.m23 = values[2][3];
988 m1.m30 = values[3][0];
989 m1.m31 = values[3][1];
990 m1.m32 = values[3][2];
991 m1.m33 = values[3][3];
997 * Places the values in the upper 4x4 of this GMatrix into
999 * @param m1 The matrix that will hold the new values
1001 public final void get(Matrix4f m1)
1004 if (nRow < 4 || nCol < 4) {
1008 m1.m00 = (float)values[0][0];
1010 m1.m10 = (float)values[1][0];
1012 m1.m20 = (float)values[2][0];
1014 m1.m30 = (float)values[3][0];
1021 m1.m01 = (float)values[0][1];
1023 m1.m11 = (float)values[1][1];
1025 m1.m21 = (float)values[2][1];
1027 m1.m31 = (float)values[3][1];
1034 m1.m02 = (float)values[0][2];
1036 m1.m12 = (float)values[1][2];
1038 m1.m22 = (float)values[2][2];
1040 m1.m32 = (float)values[3][2];
1047 m1.m03 = (float)values[0][3];
1049 m1.m13 = (float)values[1][3];
1051 m1.m23 = (float)values[2][3];
1053 m1.m33 = (float)values[3][3];
1063 m1.m00 = (float)values[0][0];
1064 m1.m01 = (float)values[0][1];
1065 m1.m02 = (float)values[0][2];
1066 m1.m03 = (float)values[0][3];
1068 m1.m10 = (float)values[1][0];
1069 m1.m11 = (float)values[1][1];
1070 m1.m12 = (float)values[1][2];
1071 m1.m13 = (float)values[1][3];
1073 m1.m20 = (float)values[2][0];
1074 m1.m21 = (float)values[2][1];
1075 m1.m22 = (float)values[2][2];
1076 m1.m23 = (float)values[2][3];
1078 m1.m30 = (float)values[3][0];
1079 m1.m31 = (float)values[3][1];
1080 m1.m32 = (float)values[3][2];
1081 m1.m33 = (float)values[3][3];
1086 * Places the values in the this GMatrix into the matrix m1;
1087 * m1 should be at least as large as this GMatrix.
1088 * @param m1 The matrix that will hold the new values
1090 public final void get(GMatrix m1)
1104 for (i = 0; i < nr; i++) {
1105 for (j = 0; j < nc; j++) {
1106 m1.values[i][j] = values[i][j];
1109 for (i = nr; i < m1.nRow; i++) {
1110 for (j = 0; j < m1.nCol; j++) {
1111 m1.values[i][j] = 0.0;
1114 for (j = nc; j < m1.nCol; j++) {
1115 for (i = 0; i < nr; i++) {
1116 m1.values[i][j] = 0.0;
1122 * Copy the values from the array into the specified row of this
1124 * @param row the row of this matrix into which the array values
1126 * @param array the source array
1128 public final void setRow(int row, double[] array)
1130 for (int i = 0; i < nCol; i++) {
1131 values[row][i] = array[i];
1136 * Copy the values from the vector into the specified row of this
1138 * @param row the row of this matrix into which the array values
1140 * @param vector the source vector
1142 public final void setRow(int row, GVector vector)
1144 for(int i = 0; i < nCol; i++) {
1145 values[row][i] = vector.values[i];
1150 * Copy the values from the array into the specified column of this
1152 * @param col the column of this matrix into which the array values
1154 * @param array the source array
1156 public final void setColumn(int col, double[] array)
1158 for(int i = 0; i < nRow; i++) {
1159 values[i][col] = array[i];
1164 * Copy the values from the vector into the specified column of this
1166 * @param col the column of this matrix into which the array values
1168 * @param vector the source vector
1170 public final void setColumn(int col, GVector vector)
1172 for(int i = 0; i < nRow; i++) {
1173 values[i][col] = vector.values[i];
1179 * Multiplies the transpose of matrix m1 times the transpose of matrix
1180 * m2, and places the result into this.
1181 * @param m1 The matrix on the left hand side of the multiplication
1182 * @param m2 The matrix on the right hand side of the multiplication
1184 public final void mulTransposeBoth(GMatrix m1, GMatrix m2)
1188 if (m1.nRow != m2.nCol || nRow != m1.nCol || nCol != m2.nRow)
1189 throw new MismatchedSizeException
1190 (VecMathI18N.getString("GMatrix14"));
1192 if (m1 == this || m2 == this) {
1193 double[][] tmp = new double[nRow][nCol];
1194 for (i = 0; i < nRow; i++) {
1195 for (j = 0; j < nCol; j++) {
1197 for (k = 0; k < m1.nRow; k++) {
1198 tmp[i][j] += m1.values[k][i]*m2.values[j][k];
1204 for (i = 0; i < nRow; i++) {
1205 for (j = 0; j < nCol; j++) {
1207 for (k = 0; k < m1.nRow; k++) {
1208 values[i][j] += m1.values[k][i]*m2.values[j][k];
1216 * Multiplies matrix m1 times the transpose of matrix m2, and
1217 * places the result into this.
1218 * @param m1 The matrix on the left hand side of the multiplication
1219 * @param m2 The matrix on the right hand side of the multiplication
1221 public final void mulTransposeRight(GMatrix m1, GMatrix m2)
1225 if (m1.nCol != m2.nCol || nCol != m2.nRow || nRow != m1.nRow)
1226 throw new MismatchedSizeException
1227 (VecMathI18N.getString("GMatrix15"));
1229 if (m1 == this || m2 == this) {
1230 double[][] tmp = new double[nRow][nCol];
1231 for (i = 0; i < nRow; i++) {
1232 for (j = 0; j < nCol; j++) {
1234 for (k = 0; k < m1.nCol; k++) {
1235 tmp[i][j] += m1.values[i][k]*m2.values[j][k];
1241 for (i = 0; i < nRow; i++) {
1242 for (j = 0;j < nCol; j++) {
1244 for (k = 0; k < m1.nCol; k++) {
1245 values[i][j] += m1.values[i][k]*m2.values[j][k];
1255 * Multiplies the transpose of matrix m1 times matrix m2, and
1256 * places the result into this.
1257 * @param m1 The matrix on the left hand side of the multiplication
1258 * @param m2 The matrix on the right hand side of the multiplication
1260 public final void mulTransposeLeft(GMatrix m1, GMatrix m2)
1264 if (m1.nRow != m2.nRow || nCol != m2.nCol || nRow != m1.nCol)
1265 throw new MismatchedSizeException
1266 (VecMathI18N.getString("GMatrix16"));
1268 if (m1 == this || m2 == this) {
1269 double[][] tmp = new double[nRow][nCol];
1270 for (i = 0; i < nRow; i++) {
1271 for (j = 0; j < nCol; j++) {
1273 for (k = 0; k < m1.nRow; k++) {
1274 tmp[i][j] += m1.values[k][i]*m2.values[k][j];
1280 for (i = 0; i < nRow; i++) {
1281 for (j = 0; j < nCol; j++) {
1283 for (k = 0; k < m1.nRow; k++) {
1284 values[i][j] += m1.values[k][i]*m2.values[k][j];
1293 * Transposes this matrix in place.
1295 public final void transpose()
1304 tmp = new double[nRow][nCol];
1305 for (i = 0; i < nRow; i++) {
1306 for (j = 0; j < nCol; j++) {
1307 tmp[i][j] = values[j][i];
1313 for (i = 0; i < nRow; i++) {
1314 for (j = 0; j < i; j++) {
1315 swap = values[i][j];
1316 values[i][j] = values[j][i];
1317 values[j][i] = swap;
1324 * Places the matrix values of the transpose of matrix m1 into this matrix.
1325 * @param m1 the matrix to be transposed (but not modified)
1327 public final void transpose(GMatrix m1)
1331 if (nRow != m1.nCol || nCol != m1.nRow)
1332 throw new MismatchedSizeException
1333 (VecMathI18N.getString("GMatrix17"));
1336 for (i = 0; i < nRow; i++) {
1337 for (j = 0;j < nCol; j++) {
1338 values[i][j] = m1.values[j][i];
1347 * Returns a string that contains the values of this GMatrix.
1348 * @return the String representation
1350 public String toString()
1352 StringBuffer buffer = new StringBuffer(nRow*nCol*8);
1356 for (i = 0; i < nRow; i++) {
1357 for (j = 0; j < nCol; j++) {
1358 buffer.append(values[i][j]).append(" ");
1360 buffer.append("\n");
1363 return buffer.toString();
1366 private static void checkMatrix( GMatrix m)
1370 for (i = 0; i < m.nRow; i++) {
1371 for (j = 0; j < m.nCol; j++) {
1372 if (Math.abs(m.values[i][j]) < 0.0000000001) {
1373 System.out.print(" 0.0 ");
1375 System.out.print(" " + m.values[i][j]);
1378 System.out.print("\n");
1384 * Returns a hash code value based on the data values in this
1385 * object. Two different GMatrix objects with identical data
1386 * values (i.e., GMatrix.equals returns true) will return the
1387 * same hash number. Two GMatrix objects with different data
1388 * members may return the same hash value, although this is not
1390 * @return the integer hash code value
1392 public int hashCode() {
1395 bits = 31L * bits + (long)nRow;
1396 bits = 31L * bits + (long)nCol;
1398 for (int i = 0; i < nRow; i++) {
1399 for (int j = 0; j < nCol; j++) {
1400 bits = 31L * bits + VecMathUtil.doubleToLongBits(values[i][j]);
1404 return (int) (bits ^ (bits >> 32));
1409 * Returns true if all of the data members of GMatrix m1 are
1410 * equal to the corresponding data members in this GMatrix.
1411 * @param m1 The matrix with which the comparison is made.
1412 * @return true or false
1414 public boolean equals(GMatrix m1)
1419 if (nRow != m1.nRow || nCol != m1.nCol)
1422 for (i = 0;i < nRow; i++) {
1423 for (j = 0; j < nCol; j++) {
1424 if (values[i][j] != m1.values[i][j])
1430 catch (NullPointerException e2) {
1436 * Returns true if the Object o1 is of type GMatrix and all of the
1437 * data members of o1 are equal to the corresponding data members in
1439 * @param o1 The object with which the comparison is made.
1440 * @return true or false
1442 public boolean equals(Object o1)
1445 GMatrix m2 = (GMatrix) o1;
1447 if (nRow != m2.nRow || nCol != m2.nCol)
1450 for (i = 0; i < nRow; i++) {
1451 for (j = 0; j < nCol; j++) {
1452 if (values[i][j] != m2.values[i][j])
1458 catch (ClassCastException e1) {
1461 catch (NullPointerException e2) {
1467 * @deprecated Use epsilonEquals(GMatrix, double) instead
1469 public boolean epsilonEquals(GMatrix m1, float epsilon) {
1470 return epsilonEquals(m1, (double)epsilon);
1474 * Returns true if the L-infinite distance between this matrix
1475 * and matrix m1 is less than or equal to the epsilon parameter,
1476 * otherwise returns false. The L-infinite
1477 * distance is equal to
1478 * MAX[i=0,1,2, . . .n ; j=0,1,2, . . .n ; abs(this.m(i,j) - m1.m(i,j)]
1479 * @param m1 The matrix to be compared to this matrix
1480 * @param epsilon the threshold value
1482 public boolean epsilonEquals(GMatrix m1, double epsilon)
1486 if (nRow != m1.nRow || nCol != m1.nCol)
1489 for (i = 0; i < nRow; i++) {
1490 for (j = 0; j < nCol; j++) {
1491 diff = values[i][j] - m1.values[i][j];
1492 if ((diff < 0 ? -diff : diff) > epsilon)
1500 * Returns the trace of this matrix.
1501 * @return the trace of this matrix
1503 public final double trace()
1514 for (i = 0; i < l; i++) {
1521 * Finds the singular value decomposition (SVD) of this matrix
1522 * such that this = U*W*transpose(V); and returns the rank of
1523 * this matrix; the values of U,W,V are all overwritten. Note
1524 * that the matrix V is output as V, and
1525 * not transpose(V). If this matrix is mxn, then U is mxm, W
1526 * is a diagonal matrix that is mxn, and V is nxn. Using the
1527 * notation W = diag(w), then the inverse of this matrix is:
1528 * inverse(this) = V*diag(1/w)*tranpose(U), where diag(1/w)
1529 * is the same matrix as W except that the reciprocal of each
1530 * of the diagonal components is used.
1531 * @param U The computed U matrix in the equation this = U*W*transpose(V)
1532 * @param W The computed W matrix in the equation this = U*W*transpose(V)
1533 * @param V The computed V matrix in the equation this = U*W*transpose(V)
1534 * @return The rank of this matrix.
1536 public final int SVD(GMatrix U, GMatrix W, GMatrix V)
1538 // check for consistancy in dimensions
1539 if (nCol != V.nCol || nCol != V.nRow) {
1540 throw new MismatchedSizeException
1541 (VecMathI18N.getString("GMatrix18"));
1544 if (nRow != U.nRow || nRow != U.nCol) {
1545 throw new MismatchedSizeException
1546 (VecMathI18N.getString("GMatrix25"));
1549 if (nRow != W.nRow || nCol != W.nCol) {
1550 throw new MismatchedSizeException
1551 (VecMathI18N.getString("GMatrix26"));
1554 // Fix ArrayIndexOutOfBounds for 2x2 matrices, which partially
1555 // addresses bug 4348562 for J3D 1.2.1.
1557 // Does *not* fix the following problems reported in 4348562,
1558 // which will wait for J3D 1.3:
1560 // 1) no output of W
1561 // 2) wrong transposition of U
1562 // 3) wrong results for 4x4 matrices
1563 // 4) slow performance
1564 if (nRow == 2 && nCol == 2) {
1565 if (values[1][0] == 0.0) {
1569 if (values[0][1] == 0.0) {
1573 double[] sinl = new double[1];
1574 double[] sinr = new double[1];
1575 double[] cosl = new double[1];
1576 double[] cosr = new double[1];
1577 double[] single_values = new double[2];
1579 single_values[0] = values[0][0];
1580 single_values[1] = values[1][1];
1582 compute_2X2(values[0][0], values[0][1], values[1][1],
1583 single_values, sinl, cosl, sinr, cosr, 0);
1585 update_u(0, U, cosl, sinl);
1586 update_v(0, V, cosr, sinr);
1590 // else call computeSVD() and check for 2x2 there
1593 return computeSVD(this, U, W, V);
1597 * LU Decomposition: this matrix must be a square matrix and the
1598 * LU GMatrix parameter must be the same size as this matrix.
1599 * The matrix LU will be overwritten as the combination of a
1600 * lower diagonal and upper diagonal matrix decompostion of this
1601 * matrix; the diagonal
1602 * elements of L (unity) are not stored. The GVector parameter
1603 * records the row permutation effected by the partial pivoting,
1604 * and is used as a parameter to the GVector method LUDBackSolve
1605 * to solve sets of linear equations.
1606 * This method returns +/- 1 depending on whether the number
1607 * of row interchanges was even or odd, respectively.
1608 * @param LU The matrix into which the lower and upper decompositions
1610 * @param permutation The row permutation effected by the partial
1612 * @return +-1 depending on whether the number of row interchanges
1613 * was even or odd respectively
1615 public final int LUD(GMatrix LU, GVector permutation)
1617 int size = LU.nRow*LU.nCol;
1618 double[] temp = new double[size];
1619 int[] even_row_exchange = new int[1];
1620 int[] row_perm = new int[LU.nRow];
1624 throw new MismatchedSizeException
1625 (VecMathI18N.getString("GMatrix19"));
1628 if (nRow != LU.nRow) {
1629 throw new MismatchedSizeException
1630 (VecMathI18N.getString("GMatrix27"));
1633 if (nCol != LU.nCol) {
1634 throw new MismatchedSizeException
1635 (VecMathI18N.getString("GMatrix27"));
1638 if (LU.nRow != permutation.getSize()) {
1639 throw new MismatchedSizeException
1640 (VecMathI18N.getString("GMatrix20"));
1643 for (i = 0; i < nRow; i++) {
1644 for (j = 0; j < nCol; j++) {
1645 temp[i*nCol+j] = values[i][j];
1649 // Calculate LU decomposition: Is the matrix singular?
1650 if (!luDecomposition(LU.nRow, temp, row_perm, even_row_exchange)) {
1651 // Matrix has no inverse
1652 throw new SingularMatrixException
1653 (VecMathI18N.getString("GMatrix21"));
1656 for (i = 0; i < nRow; i++) {
1657 for (j = 0; j < nCol; j++) {
1658 LU.values[i][j] = temp[i*nCol+j];
1662 for (i = 0; i < LU.nRow; i++){
1663 permutation.values[i] = (double)row_perm[i];
1666 return even_row_exchange[0];
1670 * Sets this matrix to a uniform scale matrix; all of the
1672 * @param scale The new scale value
1674 public final void setScale(double scale)
1683 for (i = 0; i < nRow; i++) {
1684 for (j = 0; j < nCol; j++) {
1689 for (i = 0; i < l; i++) {
1690 values[i][i] = scale;
1695 * General invert routine. Inverts m1 and places the result in "this".
1696 * Note that this routine handles both the "this" version and the
1697 * non-"this" version.
1699 * Also note that since this routine is slow anyway, we won't worry
1700 * about allocating a little bit of garbage.
1702 final void invertGeneral(GMatrix m1) {
1703 int size = m1.nRow*m1.nCol;
1704 double temp[] = new double[size];
1705 double result[] = new double[size];
1706 int row_perm[] = new int[m1.nRow];
1707 int[] even_row_exchange = new int[1];
1710 // Use LU decomposition and backsubstitution code specifically
1711 // for floating-point nxn matrices.
1712 if (m1.nRow != m1.nCol) {
1713 // Matrix is either under or over determined
1714 throw new MismatchedSizeException
1715 (VecMathI18N.getString("GMatrix22"));
1718 // Copy source matrix to temp
1719 for (i = 0; i < nRow; i++) {
1720 for (j = 0; j < nCol; j++) {
1721 temp[i*nCol+j] = m1.values[i][j];
1725 // Calculate LU decomposition: Is the matrix singular?
1726 if (!luDecomposition(m1.nRow, temp, row_perm, even_row_exchange)) {
1727 // Matrix has no inverse
1728 throw new SingularMatrixException
1729 (VecMathI18N.getString("GMatrix21"));
1732 // Perform back substitution on the identity matrix
1733 for (i = 0; i < size; i++)
1736 for (i = 0; i < nCol; i++)
1737 result[i+i*nCol] = 1.0;
1739 luBacksubstitution(m1.nRow, temp, row_perm, result);
1741 for (i = 0; i < nRow; i++) {
1742 for (j = 0; j < nCol; j++) {
1743 values[i][j] = result[i*nCol+j];
1749 * Given a nxn array "matrix0", this function replaces it with the
1750 * LU decomposition of a row-wise permutation of itself. The input
1751 * parameters are "matrix0" and "dim". The array "matrix0" is also
1752 * an output parameter. The vector "row_perm[]" is an output
1753 * parameter that contains the row permutations resulting from partial
1754 * pivoting. The output parameter "even_row_xchg" is 1 when the
1755 * number of row exchanges is even, or -1 otherwise. Assumes data
1756 * type is always double.
1758 * @return true if the matrix is nonsingular, or false otherwise.
1761 // Reference: Press, Flannery, Teukolsky, Vetterling,
1762 // _Numerical_Recipes_in_C_, Cambridge University Press,
1765 static boolean luDecomposition(int dim, double[] matrix0,
1766 int[] row_perm, int[] even_row_xchg) {
1768 double row_scale[] = new double[dim];
1770 // Determine implicit scaling information by looping over rows
1777 even_row_xchg[0] = 1;
1784 // For each column, find the largest element in the row
1787 temp = matrix0[ptr++];
1788 temp = Math.abs(temp);
1794 // Is the matrix singular?
1798 row_scale[rs++] = 1.0 / big;
1801 // For all columns, execute Crout's method
1803 for (j = 0; j < dim; j++) {
1808 // Determine elements of upper diagonal matrix U
1809 for (i = 0; i < j; i++) {
1810 target = mtx + (dim*i) + j;
1811 sum = matrix0[target];
1816 sum -= matrix0[p1] * matrix0[p2];
1820 matrix0[target] = sum;
1823 // Search for largest pivot element and calculate
1824 // intermediate elements of lower diagonal matrix L.
1827 for (i = j; i < dim; i++) {
1828 target = mtx + (dim*i) + j;
1829 sum = matrix0[target];
1834 sum -= matrix0[p1] * matrix0[p2];
1838 matrix0[target] = sum;
1840 // Is this the best pivot so far?
1841 if ((temp = row_scale[i] * Math.abs(sum)) >= big) {
1848 throw new RuntimeException(VecMathI18N.getString("GMatrix24"));
1851 // Is a row exchange necessary?
1853 // Yes: exchange rows
1855 p1 = mtx + (dim*imax);
1859 matrix0[p1++] = matrix0[p2];
1860 matrix0[p2++] = temp;
1863 // Record change in scale factor
1864 row_scale[imax] = row_scale[j];
1865 even_row_xchg[0] = -even_row_xchg[0]; // change exchange parity
1868 // Record row permutation
1871 // Is the matrix singular
1872 if (matrix0[(mtx + (dim*j) + j)] == 0.0) {
1876 // Divide elements of lower diagonal matrix L by pivot
1878 temp = 1.0 / (matrix0[(mtx + (dim*j) + j)]);
1879 target = mtx + (dim*(j+1)) + j;
1882 matrix0[target] *= temp;
1893 * Solves a set of linear equations. The input parameters "matrix1",
1894 * and "row_perm" come from luDecompostion and do not change
1895 * here. The parameter "matrix2" is a set of column vectors assembled
1896 * into a nxn matrix of floating-point values. The procedure takes each
1897 * column of "matrix2" in turn and treats it as the right-hand side of the
1898 * matrix equation Ax = LUx = b. The solution vector replaces the
1899 * original column of the matrix.
1901 * If "matrix2" is the identity matrix, the procedure replaces its contents
1902 * with the inverse of the matrix from which "matrix1" was originally
1906 // Reference: Press, Flannery, Teukolsky, Vetterling,
1907 // _Numerical_Recipes_in_C_, Cambridge University Press,
1910 static void luBacksubstitution(int dim, double[] matrix1,
1914 int i, ii, ip, j, k;
1922 // For each column vector of matrix2 ...
1923 for (k = 0; k < dim; k++) {
1924 // cv = &(matrix2[0][k]);
1928 // Forward substitution
1929 for (i = 0; i < dim; i++) {
1932 ip = row_perm[rp+i];
1933 sum = matrix2[cv+dim*ip];
1934 matrix2[cv+dim*ip] = matrix2[cv+dim*i];
1936 // rv = &(matrix1[i][0]);
1938 for (j = ii; j <= i-1; j++) {
1939 sum -= matrix1[rv+j] * matrix2[cv+dim*j];
1942 else if (sum != 0.0) {
1945 matrix2[cv+dim*i] = sum;
1949 for (i = 0; i < dim; i++) {
1954 tt += matrix1[rv+dim-j] * matrix2[cv+dim*(dim-j)];
1956 matrix2[cv+dim*ri]= (matrix2[cv+dim*ri] - tt) / matrix1[rv+ri];
1961 static int computeSVD(GMatrix mat, GMatrix U, GMatrix W, GMatrix V) {
1965 int converged, rank;
1966 double cs, sn, r, mag,scale, t;
1967 int eLength, sLength, vecLength;
1969 GMatrix tmp = new GMatrix(mat.nRow, mat.nCol);
1970 GMatrix u = new GMatrix(mat.nRow, mat.nCol);
1971 GMatrix v = new GMatrix(mat.nRow, mat.nCol);
1972 GMatrix m = new GMatrix(mat);
1974 // compute the number of singular values
1975 if (m.nRow >= m.nCol) {
1983 if (m.nRow > m.nCol)
1988 double[] vec = new double[vecLength];
1989 double[] single_values = new double[sLength];
1990 double[] e = new double[eLength];
1993 System.out.println("input to compute_svd = \n"+m.toString());
2004 // householder reduction
2005 for (si = 0; si < sLength; si++) {
2006 // for each singular value
2012 ("*********************** U ***********************\n");
2014 // compute reflector
2016 for (i = 0; i < nr; i++) {
2017 mag += m.values[i+si][si] * m.values[i+si][si];
2020 ("mag = " + mag + " matrix.dot = " +
2021 m.values[i+si][si] * m.values[i+si][si]);
2024 mag = Math.sqrt(mag);
2025 if (m.values[si][si] == 0.0) {
2028 vec[0] = m.values[si][si] + d_sign(mag, m.values[si][si]);
2031 for (i = 1; i < nr; i++) {
2032 vec[i] = m.values[si+i][si];
2036 for (i = 0; i < nr; i++) {
2038 System.out.println("vec["+i+"]="+vec[i]);
2040 scale += vec[i]*vec[i];
2045 System.out.println("scale = "+scale);
2047 for (j = si; j < m.nRow; j++) {
2048 for (k = si; k < m.nRow; k++) {
2049 u.values[j][k] = -scale * vec[j-si]*vec[k-si];
2053 for (i = si; i < m.nRow; i++){
2054 u.values[i][i] += 1.0;
2059 for (i = si; i < m.nRow; i++){
2060 t += u.values[si][i] * m.values[i][si];
2062 m.values[si][si] = t;
2065 for (j = si; j < m.nRow; j++) {
2066 for (k = si+1; k < m.nCol; k++) {
2067 tmp.values[j][k] = 0.0;
2068 for (i = si; i < m.nCol; i++) {
2069 tmp.values[j][k] += u.values[j][i] * m.values[i][k];
2074 for (j = si; j < m.nRow; j++) {
2075 for (k = si+1; k < m.nCol; k++) {
2076 m.values[j][k] = tmp.values[j][k];
2081 System.out.println("U =\n" + U.toString());
2082 System.out.println("u =\n" + u.toString());
2086 for (j = si; j < m.nRow; j++) {
2087 for (k = 0; k < m.nCol; k++) {
2088 tmp.values[j][k] = 0.0;
2089 for (i = si; i < m.nCol; i++) {
2090 tmp.values[j][k] += u.values[j][i] * U.values[i][k];
2095 for (j = si; j < m.nRow; j++) {
2096 for (k = 0; k < m.nCol; k++) {
2097 U.values[j][k] = tmp.values[j][k];
2102 System.out.println("single_values["+si+"] =\n" +
2104 System.out.println("m =\n" + m.toString());
2105 System.out.println("U =\n" + U.toString());
2115 ("*********************** V ***********************\n");
2118 for (i = 1; i < nc; i++){
2119 mag += m.values[si][si+i] * m.values[si][si+i];
2123 System.out.println("mag = " + mag);
2125 // generate the reflection vector, compute the first entry and
2126 // copy the rest from the row to be zeroed
2127 mag = Math.sqrt(mag);
2128 if (m.values[si][si+1] == 0.0) {
2131 vec[0] = m.values[si][si+1] +
2132 d_sign(mag, m.values[si][si+1]);
2135 for (i = 1; i < nc - 1; i++){
2136 vec[i] = m.values[si][si+i+1];
2139 // use reflection vector to compute v matrix
2141 for (i = 0; i < nc - 1; i++){
2142 if( debug )System.out.println("vec["+i+"]="+vec[i]);
2143 scale += vec[i]*vec[i];
2148 System.out.println("scale = "+scale);
2150 for (j = si + 1; j < nc; j++) {
2151 for (k = si+1; k < m.nCol; k++) {
2152 v.values[j][k] = -scale * vec[j-si-1]*vec[k-si-1];
2156 for (i = si + 1; i < m.nCol; i++){
2157 v.values[i][i] += 1.0;
2161 for (i = si; i < m.nCol; i++){
2162 t += v.values[i][si+1] * m.values[si][i];
2164 m.values[si][si+1]=t;
2167 for (j = si + 1; j < m.nRow; j++) {
2168 for (k = si + 1; k < m.nCol; k++) {
2169 tmp.values[j][k] = 0.0;
2170 for (i = si + 1; i < m.nCol; i++) {
2171 tmp.values[j][k] += v.values[i][k] * m.values[j][i];
2176 for (j = si + 1; j < m.nRow; j++) {
2177 for (k = si + 1; k < m.nCol; k++) {
2178 m.values[j][k] = tmp.values[j][k];
2183 System.out.println("V =\n" + V.toString());
2184 System.out.println("v =\n" + v.toString());
2185 System.out.println("tmp =\n" + tmp.toString());
2189 for (j = 0; j < m.nRow; j++) {
2190 for (k = si + 1; k < m.nCol; k++) {
2191 tmp.values[j][k] = 0.0;
2192 for (i = si + 1; i < m.nCol; i++) {
2193 tmp.values[j][k] += v.values[i][k] * V.values[j][i];
2199 System.out.println("tmp =\n" + tmp.toString());
2201 for (j = 0;j < m.nRow; j++) {
2202 for (k = si + 1; k < m.nCol; k++) {
2203 V.values[j][k] = tmp.values[j][k];
2208 System.out.println("m =\n" + m.toString());
2209 System.out.println("V =\n" + V.toString());
2216 for (i = 0; i < sLength; i++){
2217 single_values[i] = m.values[i][i];
2220 for (i = 0; i < eLength; i++){
2221 e[i] = m.values[i][i+1];
2224 // Fix ArrayIndexOutOfBounds for 2x2 matrices, which partially
2225 // addresses bug 4348562 for J3D 1.2.1.
2227 // Does *not* fix the following problems reported in 4348562,
2228 // which will wait for J3D 1.3:
2230 // 1) no output of W
2231 // 2) wrong transposition of U
2232 // 3) wrong results for 4x4 matrices
2233 // 4) slow performance
2234 if (m.nRow == 2 && m.nCol == 2) {
2235 double[] cosl = new double[1];
2236 double[] cosr = new double[1];
2237 double[] sinl = new double[1];
2238 double[] sinr = new double[1];
2240 compute_2X2(single_values[0], e[0], single_values[1],
2241 single_values, sinl, cosl, sinr, cosr, 0);
2243 update_u(0, U, cosl, sinl);
2244 update_v(0, V, cosr, sinr);
2249 // compute_qr causes ArrayIndexOutOfBounds for 2x2 matrices
2250 compute_qr (0, e.length-1, single_values, e, U, V);
2252 // compute rank = number of non zero singular values
2253 rank = single_values.length;
2255 // sort by order of size of single values
2256 // and check for zero's
2260 static void compute_qr(int start, int end, double[] s, double[] e,
2261 GMatrix u, GMatrix v) {
2265 double shift, r, utemp, vtemp, f, g;
2266 double[] cosl = new double[1];
2267 double[] cosr = new double[1];
2268 double[] sinl = new double[1];
2269 double[] sinr = new double[1];
2270 GMatrix m = new GMatrix(u.nCol, v.nRow);
2272 final int MAX_INTERATIONS = 2;
2273 final double CONVERGE_TOL = 4.89E-15;
2276 System.out.println("start =" + start);
2277 System.out.println("s =\n");
2278 for(i=0;i<s.length;i++) {
2279 System.out.println(s[i]);
2282 System.out.println("\nes =\n");
2283 for (i = 0; i < e.length; i++) {
2284 System.out.println(e[i]);
2287 for (i = 0; i < s.length; i++) {
2288 m.values[i][i] = s[i];
2291 for (i = 0; i < e.length; i++) {
2292 m.values[i][i+1] = e[i];
2294 System.out.println("\nm =\n" + m.toString());
2298 double c_b71 = -1.0;
2302 print_svd(s, e, u, v);
2307 for (k = 0; k < MAX_INTERATIONS && !converged;k++) {
2308 for (i = start; i <= end; i++) {
2310 // if at start of iterfaction compute shift
2312 if (e.length == s.length)
2317 shift = compute_shift(s[sl-1], e[end], s[sl]);
2319 f = (Math.abs(s[i]) - shift) *
2320 (d_sign(c_b48, s[i]) + shift/s[i]);
2324 r = compute_rot(f, g, sinr, cosr);
2328 f = cosr[0] * s[i] + sinr[0] * e[i];
2329 e[i] = cosr[0] * e[i] - sinr[0] * s[i];
2330 g = sinr[0] * s[i+1];
2331 s[i+1] = cosr[0] * s[i+1];
2333 // if (debug) print_se(s,e);
2334 update_v (i, v, cosr, sinr);
2338 r = compute_rot(f, g, sinl, cosl);
2340 f = cosl[0] * e[i] + sinl[0] * s[i+1];
2341 s[i+1] = cosl[0] * s[i+1] - sinl[0] * e[i];
2345 g = sinl[0] * e[i+1];
2346 e[i+1] = cosl[0] * e[i+1];
2348 //if (debug) print_se(s,e);
2350 update_u(i, u, cosl, sinl);
2355 // if extra off diagonal perform one more right side rotation
2356 if (s.length == e.length) {
2357 r = compute_rot(f, g, sinr, cosr);
2358 f = cosr[0] * s[i] + sinr[0] * e[i];
2359 e[i] = cosr[0] * e[i] - sinr[0] * s[i];
2360 s[i+1] = cosr[0] * s[i+1];
2362 update_v(i, v, cosr, sinr);
2369 ("\n*********************** iteration #" + k +
2370 " ***********************\n");
2371 print_svd(s, e, u, v);
2374 // check for convergence on off diagonals and reduce
2375 while ((end-start > 1) && (Math.abs(e[end]) < CONVERGE_TOL)) {
2379 // check if need to split
2380 for (n = end - 2; n > start; n--) {
2381 if (Math.abs(e[n]) < CONVERGE_TOL) { // split
2382 compute_qr(n + 1, end, s, e, u, v); // do lower matrix
2383 end = n - 1; // do upper matrix
2385 // check for convergence on off diagonals and reduce
2386 while ((end - start > 1) &&
2387 (Math.abs(e[end]) < CONVERGE_TOL)) {
2394 System.out.println("start = " + start);
2396 if ((end - start <= 1) && (Math.abs(e[start+1]) < CONVERGE_TOL)) {
2399 // check if zero on the diagonal
2405 System.out.println("\n****call compute_2X2 ********************\n");
2407 if (Math.abs(e[1]) < CONVERGE_TOL) {
2408 compute_2X2(s[start], e[start], s[start+1], s,
2409 sinl, cosl, sinr, cosr, 0);
2416 update_u(i, u, cosl, sinl);
2417 update_v(i, v, cosr, sinr);
2421 ("\n*******after call compute_2X2 **********************\n");
2422 print_svd(s, e, u, v);
2428 private static void print_se(double[] s, double[] e) {
2429 System.out.println("\ns =" + s[0] + " " + s[1] + " " + s[2]);
2430 System.out.println("e =" + e[0] + " " + e[1]);
2433 private static void update_v(int index, GMatrix v,
2434 double[] cosr, double[] sinr) {
2438 for (j = 0; j < v.nRow; j++) {
2439 vtemp = v.values[j][index];
2440 v.values[j][index] =
2441 cosr[0]*vtemp + sinr[0]*v.values[j][index+1];
2442 v.values[j][index+1] =
2443 -sinr[0]*vtemp + cosr[0]*v.values[j][index+1];
2447 private static void chase_up(double[] s, double[] e, int k, GMatrix v) {
2449 double[] cosr = new double[1];
2450 double[] sinr = new double[1];
2452 GMatrix t = new GMatrix(v.nRow, v.nCol);
2453 GMatrix m = new GMatrix(v.nRow, v.nCol);
2457 for (i = 0; i < s.length; i++) {
2458 m.values[i][i] = s[i];
2460 for (i = 0; i < e.length; i++) {
2461 m.values[i][i+1] = e[i];
2468 for (i = k; i > 0; i--) {
2469 r = compute_rot(f, g, sinr, cosr);
2470 f = -e[i-1] * sinr[0];
2473 e[i-1] = e[i-1] * cosr[0];
2474 update_v_split(i, k+1, v, cosr, sinr, t, m);
2477 s[i+1] = compute_rot(f, g, sinr, cosr);
2478 update_v_split(i, k+1, v, cosr, sinr, t, m);
2481 private static void chase_across(double[] s, double[] e, int k, GMatrix u) {
2483 double[] cosl = new double[1];
2484 double[] sinl = new double[1];
2486 GMatrix t = new GMatrix(u.nRow, u.nCol);
2487 GMatrix m = new GMatrix(u.nRow, u.nCol);
2491 for (i = 0; i < s.length; i++) {
2492 m.values[i][i] = s[i];
2494 for (i = 0; i < e.length; i++) {
2495 m.values[i][i+1] = e[i];
2502 for (i = k; i < u.nCol-2; i++){
2503 r = compute_rot(f, g, sinl, cosl);
2504 g = -e[i+1] * sinl[0];
2507 e[i+1] = e[i+1] * cosl[0];
2508 update_u_split(k, i + 1, u, cosl, sinl, t, m);
2511 s[i+1] = compute_rot(f, g, sinl, cosl);
2512 update_u_split(k, i + 1, u, cosl, sinl, t, m);
2515 private static void update_v_split(int topr, int bottomr, GMatrix v,
2516 double[] cosr, double[] sinr,
2517 GMatrix t, GMatrix m) {
2521 for (j = 0; j < v.nRow; j++) {
2522 vtemp = v.values[j][topr];
2523 v.values[j][topr] = cosr[0]*vtemp - sinr[0]*v.values[j][bottomr];
2524 v.values[j][bottomr] = sinr[0]*vtemp + cosr[0]*v.values[j][bottomr];
2529 for (j = 0; j < v.nRow; j++) {
2530 vtemp = t.values[j][topr];
2532 cosr[0]*vtemp - sinr[0]*t.values[j][bottomr];
2533 t.values[j][bottomr] =
2534 sinr[0]*vtemp + cosr[0]*t.values[j][bottomr];
2538 System.out.println("topr =" + topr);
2539 System.out.println("bottomr =" + bottomr);
2540 System.out.println("cosr =" + cosr[0]);
2541 System.out.println("sinr =" + sinr[0]);
2542 System.out.println("\nm =");
2544 System.out.println("\nv =");
2547 System.out.println("\nt*m =");
2551 private static void update_u_split(int topr, int bottomr, GMatrix u,
2552 double[] cosl, double[] sinl,
2553 GMatrix t, GMatrix m) {
2557 for (j = 0; j < u.nCol; j++) {
2558 utemp = u.values[topr][j];
2559 u.values[topr][j] = cosl[0]*utemp - sinl[0]*u.values[bottomr][j];
2560 u.values[bottomr][j] = sinl[0]*utemp + cosl[0]*u.values[bottomr][j];
2565 for (j = 0;j < u.nCol; j++) {
2566 utemp = t.values[topr][j];
2568 cosl[0]*utemp - sinl[0]*t.values[bottomr][j];
2569 t.values[bottomr][j] =
2570 sinl[0]*utemp + cosl[0]*t.values[bottomr][j];
2573 System.out.println("\nm=");
2575 System.out.println("\nu=");
2578 System.out.println("\nt*m=");
2582 private static void update_u(int index, GMatrix u,
2583 double[] cosl, double[] sinl) {
2587 for (j = 0; j < u.nCol; j++) {
2588 utemp = u.values[index][j];
2589 u.values[index][j] =
2590 cosl[0]*utemp + sinl[0]*u.values[index+1][j];
2591 u.values[index+1][j] =
2592 -sinl[0]*utemp + cosl[0]*u.values[index+1][j];
2596 private static void print_m(GMatrix m, GMatrix u, GMatrix v) {
2597 GMatrix mtmp = new GMatrix(m.nCol, m.nRow);
2601 System.out.println("\n m = \n" + mtmp.toString(mtmp));
2605 private static String toString(GMatrix m)
2607 StringBuffer buffer = new StringBuffer(m.nRow * m.nCol * 8);
2610 for (i = 0; i < m.nRow; i++) {
2611 for(j = 0; j < m.nCol; j++) {
2612 if (Math.abs(m.values[i][j]) < .000000001) {
2613 buffer.append("0.0000 ");
2615 buffer.append(m.values[i][j]).append(" ");
2618 buffer.append("\n");
2620 return buffer.toString();
2623 private static void print_svd(double[] s, double[] e,
2624 GMatrix u, GMatrix v) {
2626 GMatrix mtmp = new GMatrix(u.nCol, v.nRow);
2628 System.out.println(" \ns = ");
2629 for (i = 0; i < s.length; i++) {
2630 System.out.println(" " + s[i]);
2633 System.out.println(" \ne = ");
2634 for (i = 0; i < e.length; i++) {
2635 System.out.println(" " + e[i]);
2638 System.out.println(" \nu = \n" + u.toString());
2639 System.out.println(" \nv = \n" + v.toString());
2642 for (i = 0; i < s.length; i++) {
2643 mtmp.values[i][i] = s[i];
2645 for (i = 0; i < e.length; i++) {
2646 mtmp.values[i][i+1] = e[i];
2648 System.out.println(" \nm = \n"+mtmp.toString());
2650 mtmp.mulTransposeLeft(u, mtmp);
2651 mtmp.mulTransposeRight(mtmp, v);
2653 System.out.println(" \n u.transpose*m*v.transpose = \n" +
2657 static double max(double a, double b) {
2664 static double min(double a, double b) {
2671 static double compute_shift(double f, double g, double h) {
2673 double fhmn, fhmx, c, fa, ga, ha, as, at, au;
2686 d__1 = min(fhmx,ga) / max(fhmx,ga);
2690 as = fhmn / fhmx + 1.0;
2691 at = (fhmx - fhmn) / fhmx;
2694 c = 2.0 / (Math.sqrt(as * as + au) + Math.sqrt(at * at + au));
2699 ssmin = fhmn * fhmx / ga;
2701 as = fhmn / fhmx + 1.0;
2702 at = (fhmx - fhmn) / fhmx;
2705 c = 1.0 / (Math.sqrt(d__1 * d__1 + 1.0) +
2706 Math.sqrt(d__2 * d__2 + 1.0));
2707 ssmin = fhmn * c * au;
2716 static int compute_2X2(double f, double g, double h,
2717 double[] single_values, double[] snl, double[] csl,
2718 double[] snr, double[] csr, int index) {
2727 double a, d, l, m, r, s, t, tsign, fa, ga, ha;
2728 double ft, gt, ht, mm;
2730 double tt, clt, crt, slt, srt;
2733 ssmax = single_values[0];
2734 ssmin = single_values[1];
2766 single_values[1] = ha;
2767 single_values[0] = fa;
2776 if (fa / ga < EPS) {
2781 ssmin = fa / (ga / ha);
2783 ssmin = fa / ga * ha;
2804 s = Math.sqrt(tt + mm);
2809 r = Math.sqrt(l * l + mm);
2815 if (fa / ga < EPS) {
2819 ssmin = fa / (ga / ha);
2821 ssmin = fa / ga * ha;
2842 s = Math.sqrt(tt + mm);
2847 r = Math.sqrt(l * l + mm);
2856 t = d_sign(c_b3, ft) * d_sign(c_b4, gt);
2858 t = gt / d_sign(d, ft) + m / t;
2861 t = (m / (s + t) + m / (r + l)) * (a + 1.0);
2864 l = Math.sqrt(t * t + 4.0);
2867 clt = (crt + srt * m) / a;
2868 slt = ht / ft * srt / a;
2884 tsign = d_sign(c_b4, csr[0]) *
2885 d_sign(c_b4, csl[0]) * d_sign(c_b4, f);
2888 tsign = d_sign(c_b4, snr[0]) *
2889 d_sign(c_b4, csl[0]) * d_sign(c_b4, g);
2892 tsign = d_sign(c_b4, snr[0]) *
2893 d_sign(c_b4, snl[0]) * d_sign(c_b4, h);
2896 single_values[index] = d_sign(ssmax, tsign);
2897 d__1 = tsign * d_sign(c_b4, f) * d_sign(c_b4, h);
2898 single_values[index+1] = d_sign(ssmin, d__1);
2904 static double compute_rot(double f, double g, double[] sin, double[] cos) {
2913 final double safmn2 = 2.002083095183101E-146;
2914 final double safmx2 = 4.994797680505588E+145;
2920 } else if (f == 0.0) {
2927 scale = max(Math.abs(f1),Math.abs(g1));
2928 if (scale >= safmx2) {
2930 while(scale >= safmx2) {
2934 scale = max(Math.abs(f1), Math.abs(g1));
2936 r = Math.sqrt(f1*f1 + g1*g1);
2940 for (i = 1; i <= count; ++i) {
2943 } else if (scale <= safmn2) {
2945 while(scale <= safmn2) {
2949 scale = max(Math.abs(f1), Math.abs(g1));
2951 r = Math.sqrt(f1*f1 + g1*g1);
2955 for (i = 1; i <= count; ++i) {
2959 r = Math.sqrt(f1*f1 + g1*g1);
2963 if (Math.abs(f) > Math.abs(g) && cs < 0.0) {
2974 static double d_sign(double a, double b) {
2976 x = (a >= 0 ? a : - a);
2977 return (b >= 0 ? x : -x);
2981 * Creates a new object of the same class as this object.
2983 * @return a clone of this instance.
2984 * @exception OutOfMemoryError if there is not enough memory.
2985 * @see java.lang.Cloneable
2986 * @since vecmath 1.3
2988 public Object clone() {
2991 m1 = (GMatrix)super.clone();
2992 } catch (CloneNotSupportedException e) {
2993 // this shouldn't happen, since we are Cloneable
2994 throw new InternalError();
2997 // Also need to clone array of values
2998 m1.values = new double[nRow][nCol];
2999 for (int i = 0; i < nRow; i++) {
3000 for(int j = 0; j < nCol; j++) {
3001 m1.values[i][j] = values[i][j];