javax.vecmath/src/javax/vecmath/GMatrix.java

   1 /*
   2  * $RCSfile: GMatrix.java,v $
   3  *
   4  * Copyright 1997-2008 Sun Microsystems, Inc.  All Rights Reserved.
   5  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   6  *
   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.
  12  *
  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).
  18  *
  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.
  22  *
  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
  25  * have any questions.
  26  *
  27  * $Revision: 1.7 $
  28  * $Date: 2008/02/28 20:18:50 $
  29  * $State: Exp $
  30  */
  31
  32 package javax.vecmath;
  33
  34 import java.lang.Math;
  35
  36 /**
  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.
  40  */
  41
  42 public class GMatrix implements java.io.Serializable, Cloneable {
  43
  44     // Compatible with 1.1
  45     static final long serialVersionUID = 2777097312029690941L;
  46     private static final boolean debug = false;
  47
  48     int nRow;
  49     int nCol;
  50
  51     // double dereference is slow
  52     double[][] values;
  53
  54     private static final double EPS = 1.0E-10;
  55
  56     /**
  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.
  63      */
  64     public GMatrix(int nRow, int nCol)
  65     {
  66         values = new double[nRow][nCol];
  67         this.nRow = nRow;
  68         this.nCol = nCol;
  69
  70         int i, j;
  71         for (i = 0; i < nRow; i++) {
  72             for (j = 0; j < nCol; j++) {
  73                 values[i][j] = 0.0;
  74             }
  75         }
  76
  77         int l;
  78         if (nRow < nCol)
  79             l = nRow;
  80         else
  81             l = nCol;
  82
  83         for (i = 0; i < l; i++) {
  84             values[i][i] = 1.0;
  85         }
  86     }
  87
  88     /**
  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
  99      */
 100     public GMatrix(int nRow, int nCol, double[] matrix)
 101     {
 102         values = new double[nRow][nCol];
 103         this.nRow = nRow;
 104         this.nCol = nCol;
 105
 106         int i, j;
 107         for (i = 0; i < nRow; i++) {
 108             for (j = 0; j < nCol; j++) {
 109                 values[i][j] = matrix[i*nCol+j];
 110             }
 111         }
 112     }
 113
 114     /**
 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
 118      */
 119     public GMatrix(GMatrix matrix)
 120     {
 121         nRow = matrix.nRow;
 122         nCol = matrix.nCol;
 123         values = new double[nRow][nCol];
 124
 125         int i, j;
 126         for (i = 0; i < nRow; i++) {
 127             for (j = 0; j < nCol; j++) {
 128                 values[i][j] = matrix.values[i][j];
 129             }
 130         }
 131     }
 132
 133     /**
 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
 137      */
 138     public final void mul(GMatrix m1)
 139     {
 140         int i, j, k;
 141
 142         if (nCol != m1.nRow ||  nCol != m1.nCol)
 143             throw new MismatchedSizeException
 144                 (VecMathI18N.getString("GMatrix0"));
 145
 146         double [][] tmp = new double[nRow][nCol];
 147
 148         for (i = 0; i < nRow; i++) {
 149             for (j = 0; j < nCol; j++) {
 150                 tmp[i][j] = 0.0;
 151                 for (k = 0; k < nCol; k++) {
 152                     tmp[i][j] += values[i][k]*m1.values[k][j];
 153                 }
 154             }
 155         }
 156
 157         values = tmp;
 158     }
 159
 160     /**
 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
 165      */
 166     public final void mul(GMatrix m1, GMatrix m2)
 167     {
 168         int i, j, k;
 169
 170         if (m1.nCol != m2.nRow || nRow != m1.nRow || nCol != m2.nCol)
 171             throw new MismatchedSizeException
 172                 (VecMathI18N.getString("GMatrix1"));
 173
 174         double[][] tmp = new double[nRow][nCol];
 175
 176         for (i = 0; i < m1.nRow; i++) {
 177             for (j = 0; j < m2.nCol; j++) {
 178                 tmp[i][j] = 0.0;
 179                 for (k = 0; k < m1.nCol; k++) {
 180                     tmp[i][j] += m1.values[i][k]*m2.values[k][j];
 181                 }
 182             }
 183         }
 184
 185         values = tmp;
 186     }
 187
 188     /**
 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
 195      */
 196     public final void mul(GVector v1, GVector v2)
 197     {
 198         int i, j;
 199
 200         if (nRow < v1.getSize())
 201             throw new MismatchedSizeException
 202                 (VecMathI18N.getString("GMatrix2"));
 203
 204         if (nCol < v2.getSize())
 205             throw new MismatchedSizeException
 206                 (VecMathI18N.getString("GMatrix3"));
 207
 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];
 211             }
 212         }
 213     }
 214
 215     /**
 216      * Sets the value of this matrix to sum of itself and matrix m1.
 217      * @param m1 the other matrix
 218      */
 219     public final void add(GMatrix m1)
 220     {
 221         int i, j;
 222
 223         if (nRow != m1.nRow)
 224             throw new MismatchedSizeException
 225                 (VecMathI18N.getString("GMatrix4"));
 226
 227         if (nCol != m1.nCol)
 228             throw new MismatchedSizeException
 229                 (VecMathI18N.getString("GMatrix5"));
 230
 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];
 234             }
 235         }
 236     }
 237
 238     /**
 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
 242      */
 243     public final void add(GMatrix m1, GMatrix m2)
 244     {
 245         int i, j;
 246
 247         if (m2.nRow != m1.nRow)
 248             throw new MismatchedSizeException
 249                 (VecMathI18N.getString("GMatrix6"));
 250
 251         if (m2.nCol != m1.nCol)
 252             throw new MismatchedSizeException
 253                 (VecMathI18N.getString("GMatrix7"));
 254
 255         if (nCol != m1.nCol  || nRow != m1.nRow)
 256             throw new MismatchedSizeException
 257                 (VecMathI18N.getString("GMatrix8"));
 258
 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];
 262             }
 263         }
 264     }
 265
 266     /**
 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
 270      */
 271     public final void sub(GMatrix m1)
 272     {
 273         int i, j;
 274         if (nRow != m1.nRow)
 275             throw new MismatchedSizeException
 276                 (VecMathI18N.getString("GMatrix9"));
 277
 278         if (nCol != m1.nCol)
 279             throw new MismatchedSizeException
 280                 (VecMathI18N.getString("GMatrix28"));
 281
 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];
 285             }
 286         }
 287     }
 288
 289     /**
 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
 294      */
 295     public final void sub(GMatrix m1, GMatrix m2)
 296     {
 297         int i, j;
 298         if (m2.nRow != m1.nRow)
 299             throw new MismatchedSizeException
 300                 (VecMathI18N.getString("GMatrix10"));
 301
 302         if (m2.nCol != m1.nCol)
 303             throw new MismatchedSizeException
 304                 (VecMathI18N.getString("GMatrix11"));
 305
 306         if (nRow !=  m1.nRow || nCol != m1.nCol)
 307             throw new MismatchedSizeException
 308                 (VecMathI18N.getString("GMatrix12"));
 309
 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];
 313             }
 314         }
 315     }
 316
 317     /**
 318      * Negates the value of this matrix: this = -this.
 319      */
 320     public final void negate()
 321     {
 322         int i, j;
 323         for (i = 0; i < nRow; i++) {
 324             for (j = 0;j < nCol; j++) {
 325                 values[i][j] = -values[i][j];
 326             }
 327         }
 328     }
 329
 330     /**
 331      *  Sets the value of this matrix equal to the negation of
 332      *  of the GMatrix parameter.
 333      *  @param m1  The source matrix
 334      */
 335     public final void negate(GMatrix m1)
 336     {
 337         int i, j;
 338         if (nRow != m1.nRow  || nCol != m1.nCol)
 339             throw new MismatchedSizeException
 340                 (VecMathI18N.getString("GMatrix13"));
 341
 342         for (i = 0; i < nRow; i++) {
 343             for (j = 0; j < nCol; j++) {
 344                 values[i][j] =  -m1.values[i][j];
 345             }
 346         }
 347     }
 348
 349     /**
 350      * Sets this GMatrix to the identity matrix.
 351      */
 352     public final void setIdentity()
 353     {
 354         int i, j;
 355         for (i = 0; i < nRow; i++) {
 356             for (j = 0; j < nCol; j++) {
 357                 values[i][j] = 0.0;
 358             }
 359         }
 360
 361         int l;
 362         if (nRow < nCol)
 363             l = nRow;
 364         else
 365             l = nCol;
 366
 367         for (i = 0; i < l; i++) {
 368             values[i][i] = 1.0;
 369         }
 370     }
 371
 372     /**
 373      * Sets all the values in this matrix to zero.
 374      */
 375     public final void setZero()
 376     {
 377         int i, j;
 378         for (i = 0; i < nRow; i++) {
 379             for (j = 0; j < nCol; j++) {
 380                 values[i][j] = 0.0;
 381             }
 382         }
 383     }
 384
 385     /**
 386      * Subtracts this matrix from the identity matrix and puts the values
 387      * back into this (this = I - this).
 388      */
 389     public final void identityMinus()
 390     {
 391         int i, j;
 392
 393         for(i = 0; i < nRow; i++) {
 394             for(j = 0; j < nCol; j++) {
 395                 values[i][j] = -values[i][j];
 396             }
 397         }
 398
 399         int l;
 400         if( nRow < nCol)
 401             l = nRow;
 402         else
 403             l = nCol;
 404
 405         for(i = 0; i < l; i++) {
 406             values[i][i] += 1.0;
 407         }
 408     }
 409
 410
 411     /**
 412      * Inverts this matrix in place.
 413      */
 414     public final void invert()
 415     {
 416         invertGeneral(this);
 417     }
 418
 419     /**
 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
 423      */
 424     public final void invert(GMatrix m1)
 425     {
 426         invertGeneral(m1);
 427     }
 428
 429     /**
 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
 444      */
 445     public final void copySubMatrix(int rowSource, int colSource,
 446                                     int numRow, int numCol, int rowDest,
 447                                     int colDest, GMatrix target)
 448     {
 449         int i, j;
 450
 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];
 456                 }
 457             }
 458         } else {
 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];
 463                 }
 464             }
 465             for (i = 0; i < numRow; i++) {
 466                 for (j = 0; j < numCol; j++) {
 467                     target.values[rowDest+i][colDest+j] = tmp[i][j];
 468                 }
 469             }
 470         }
 471     }
 472
 473     /**
 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
 479      */
 480     public final void setSize(int nRow, int nCol)
 481     {
 482         double[][] tmp = new double[nRow][nCol];
 483         int i, j, maxRow, maxCol;
 484
 485         if (this.nRow < nRow)
 486             maxRow = this.nRow;
 487         else
 488             maxRow = nRow;
 489
 490         if (this.nCol < nCol)
 491             maxCol = this.nCol;
 492         else
 493             maxCol = nCol;
 494
 495         for (i = 0; i < maxRow; i++) {
 496             for (j = 0; j < maxCol; j++) {
 497                 tmp[i][j] = values[i][j];
 498             }
 499         }
 500
 501         this.nRow = nRow;
 502         this.nCol = nCol;
 503
 504         values = tmp;
 505     }
 506
 507     /**
 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
 512      * in this matrix.
 513      * @param matrix  the row major source array
 514      */
 515     public final void set(double[] matrix)
 516     {
 517         int i, j;
 518
 519         for (i = 0; i < nRow; i++) {
 520             for (j = 0; j < nCol; j++) {
 521                 values[i][j] = matrix[nCol*i+j];
 522             }
 523         }
 524     }
 525
 526     /**
 527      * Sets the value of this matrix to that of the Matrix3f provided.
 528      * @param m1 the matrix
 529      */
 530     public final void set(Matrix3f m1)
 531     {
 532         int i, j;
 533
 534         if (nCol < 3 || nRow < 3) { // expand matrix if too small
 535             nCol = 3;
 536             nRow = 3;
 537             values = new double[nRow][nCol];
 538         }
 539
 540         values[0][0] = m1.m00;
 541         values[0][1] = m1.m01;
 542         values[0][2] = m1.m02;
 543
 544         values[1][0] = m1.m10;
 545         values[1][1] = m1.m11;
 546         values[1][2] = m1.m12;
 547
 548         values[2][0] = m1.m20;
 549         values[2][1] = m1.m21;
 550         values[2][2] = m1.m22;
 551
 552         for (i = 3; i < nRow; i++) {   // pad rest or matrix with zeros
 553             for (j = 3; j < nCol; j++) {
 554                 values[i][j] = 0.0;
 555             }
 556         }
 557     }
 558
 559     /**
 560      * Sets the value of this matrix to that of the Matrix3d provided.
 561      * @param m1 the matrix
 562      */
 563     public final void set(Matrix3d m1)
 564     {
 565         if (nRow < 3 || nCol < 3) {
 566             values = new double[3][3];
 567             nRow = 3;
 568             nCol = 3;
 569         }
 570
 571         values[0][0] = m1.m00;
 572         values[0][1] = m1.m01;
 573         values[0][2] = m1.m02;
 574
 575         values[1][0] = m1.m10;
 576         values[1][1] = m1.m11;
 577         values[1][2] = m1.m12;
 578
 579         values[2][0] = m1.m20;
 580         values[2][1] = m1.m21;
 581         values[2][2] = m1.m22;
 582
 583         for (int i = 3; i < nRow; i++) {   // pad rest or matrix with zeros
 584             for(int j = 3; j < nCol; j++) {
 585                 values[i][j] = 0.0;
 586             }
 587         }
 588
 589     }
 590
 591     /**
 592      * Sets the value of this matrix to that of the Matrix4f provided.
 593      * @param m1 the matrix
 594      */
 595     public final void set(Matrix4f m1)
 596     {
 597         if (nRow < 4 || nCol < 4) {
 598             values = new double[4][4];
 599             nRow = 4;
 600             nCol = 4;
 601         }
 602
 603         values[0][0] = m1.m00;
 604         values[0][1] = m1.m01;
 605         values[0][2] = m1.m02;
 606         values[0][3] = m1.m03;
 607
 608         values[1][0] = m1.m10;
 609         values[1][1] = m1.m11;
 610         values[1][2] = m1.m12;
 611         values[1][3] = m1.m13;
 612
 613         values[2][0] = m1.m20;
 614         values[2][1] = m1.m21;
 615         values[2][2] = m1.m22;
 616         values[2][3] = m1.m23;
 617
 618         values[3][0] = m1.m30;
 619         values[3][1] = m1.m31;
 620         values[3][2] = m1.m32;
 621         values[3][3] = m1.m33;
 622
 623         for (int i = 4 ; i < nRow; i++) {   // pad rest or matrix with zeros
 624             for (int j = 4; j < nCol; j++) {
 625                 values[i][j] = 0.0;
 626             }
 627         }
 628     }
 629
 630     /**
 631      * Sets the value of this matrix to that of the Matrix4d provided.
 632      * @param m1 the matrix
 633      */
 634     public final void set(Matrix4d m1)
 635     {
 636         if (nRow < 4 || nCol < 4) {
 637             values = new double[4][4];
 638             nRow = 4;
 639             nCol = 4;
 640         }
 641
 642         values[0][0] = m1.m00;
 643         values[0][1] = m1.m01;
 644         values[0][2] = m1.m02;
 645         values[0][3] = m1.m03;
 646
 647         values[1][0] = m1.m10;
 648         values[1][1] = m1.m11;
 649         values[1][2] = m1.m12;
 650         values[1][3] = m1.m13;
 651
 652         values[2][0] = m1.m20;
 653         values[2][1] = m1.m21;
 654         values[2][2] = m1.m22;
 655         values[2][3] = m1.m23;
 656
 657         values[3][0] = m1.m30;
 658         values[3][1] = m1.m31;
 659         values[3][2] = m1.m32;
 660         values[3][3] = m1.m33;
 661
 662         for (int i = 4; i < nRow; i++) {   // pad rest or matrix with zeros
 663             for (int j = 4; j < nCol; j++) {
 664                 values[i][j] = 0.0;
 665             }
 666         }
 667     }
 668
 669     /**
 670      * Sets the value of this matrix to the values found in matrix m1.
 671      * @param m1  the source matrix
 672      */
 673     public final void set(GMatrix m1)
 674     {
 675         int i, j;
 676
 677         if (nRow < m1.nRow || nCol < m1.nCol) {
 678             nRow = m1.nRow;
 679             nCol = m1.nCol;
 680             values = new double[nRow][nCol];
 681         }
 682
 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];
 686             }
 687         }
 688
 689         for (i = m1.nRow; i < nRow; i++) {   // pad rest or matrix with zeros
 690             for (j = m1.nCol; j < nCol; j++) {
 691                 values[i][j] = 0.0;
 692             }
 693         }
 694     }
 695
 696     /**
 697      * Returns the number of rows in this matrix.
 698      * @return  number of rows in this matrix
 699      */
 700     public final int getNumRow()
 701     {
 702         return(nRow);
 703     }
 704
 705     /**
 706      * Returns the number of colmuns in this matrix.
 707      * @return  number of columns in this matrix
 708      */
 709     public final int getNumCol()
 710     {
 711         return(nCol);
 712     }
 713
 714     /**
 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
 719      */
 720     public final double getElement(int row, int column)
 721     {
 722         return(values[row][column]);
 723     }
 724
 725
 726     /**
 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
 731      */
 732     public final void setElement(int row, int column, double value)
 733     {
 734         values[row][column] = value;
 735     }
 736
 737     /**
 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
 741      */
 742     public final void getRow(int row, double[] array)
 743     {
 744         for (int i = 0; i < nCol; i++) {
 745             array[i] = values[row][i];
 746         }
 747     }
 748
 749     /**
 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
 753      */
 754     public final void getRow(int row, GVector vector)
 755     {
 756         if (vector.getSize() < nCol)
 757             vector.setSize(nCol);
 758
 759         for (int i = 0; i < nCol; i++) {
 760             vector.values[i] = values[row][i];
 761         }
 762     }
 763
 764     /**
 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
 768      */
 769     public final void getColumn(int col, double[] array)
 770     {
 771         for (int i = 0; i < nRow; i++) {
 772             array[i] = values[i][col];
 773         }
 774
 775     }
 776
 777     /**
 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
 781      */
 782     public final void getColumn(int col, GVector vector)
 783     {
 784         if (vector.getSize() < nRow)
 785             vector.setSize(nRow);
 786
 787         for (int i = 0; i < nRow; i++) {
 788             vector.values[i] = values[i][col];
 789         }
 790     }
 791
 792     /**
 793      * Places the values in the upper 3x3 of this GMatrix into
 794      * the matrix m1.
 795      * @param m1  The matrix that will hold the new values
 796      */
 797     public final void get(Matrix3d m1)
 798     {
 799         if (nRow < 3 || nCol < 3) {
 800             m1.setZero();
 801             if (nCol > 0) {
 802                 if (nRow > 0){
 803                     m1.m00 = values[0][0];
 804                     if (nRow > 1){
 805                         m1.m10 = values[1][0];
 806                         if( nRow > 2 ){
 807                             m1.m20= values[2][0];
 808                         }
 809                     }
 810                 }
 811                 if (nCol > 1) {
 812                     if (nRow > 0) {
 813                         m1.m01 = values[0][1];
 814                         if (nRow > 1){
 815                             m1.m11 = values[1][1];
 816                             if (nRow >  2){
 817                                 m1.m21 = values[2][1];
 818                             }
 819                         }
 820                     }
 821                     if (nCol > 2) {
 822                         if (nRow > 0) {
 823                             m1.m02 = values[0][2];
 824                             if (nRow > 1) {
 825                                 m1.m12 = values[1][2];
 826                                 if (nRow > 2) {
 827                                     m1.m22 = values[2][2];
 828                                 }
 829                             }
 830                         }
 831                     }
 832                 }
 833             }
 834         } else {
 835             m1.m00 = values[0][0];
 836             m1.m01 = values[0][1];
 837             m1.m02 = values[0][2];
 838
 839             m1.m10 = values[1][0];
 840             m1.m11 = values[1][1];
 841             m1.m12 = values[1][2];
 842
 843             m1.m20 = values[2][0];
 844             m1.m21 = values[2][1];
 845             m1.m22 = values[2][2];
 846         }
 847     }
 848
 849     /**
 850      * Places the values in the upper 3x3 of this GMatrix into
 851      * the matrix m1.
 852      * @param m1  The matrix that will hold the new values
 853      */
 854     public final void get(Matrix3f m1)
 855     {
 856
 857         if (nRow < 3 || nCol < 3) {
 858             m1.setZero();
 859             if (nCol > 0) {
 860                 if (nRow > 0) {
 861                     m1.m00 = (float)values[0][0];
 862                     if (nRow > 1) {
 863                         m1.m10 = (float)values[1][0];
 864                         if (nRow > 2) {
 865                             m1.m20 = (float)values[2][0];
 866                         }
 867                     }
 868                 }
 869                 if (nCol > 1) {
 870                     if (nRow > 0) {
 871                         m1.m01 = (float)values[0][1];
 872                         if (nRow >  1){
 873                             m1.m11 = (float)values[1][1];
 874                             if (nRow >  2){
 875                                 m1.m21 = (float)values[2][1];
 876                             }
 877                         }
 878                     }
 879                     if (nCol > 2) {
 880                         if (nRow > 0) {
 881                             m1.m02 = (float)values[0][2];
 882                             if (nRow > 1) {
 883                                 m1.m12 = (float)values[1][2];
 884                                 if (nRow > 2) {
 885                                     m1.m22 = (float)values[2][2];
 886                                 }
 887                             }
 888                         }
 889                     }
 890                 }
 891             }
 892         } else {
 893             m1.m00 = (float)values[0][0];
 894             m1.m01 = (float)values[0][1];
 895             m1.m02 = (float)values[0][2];
 896
 897             m1.m10 = (float)values[1][0];
 898             m1.m11 = (float)values[1][1];
 899             m1.m12 = (float)values[1][2];
 900
 901             m1.m20 = (float)values[2][0];
 902             m1.m21 = (float)values[2][1];
 903             m1.m22 = (float)values[2][2];
 904         }
 905     }
 906
 907     /**
 908      * Places the values in the upper 4x4 of this GMatrix into
 909      * the matrix m1.
 910      * @param m1  The matrix that will hold the new values
 911      */
 912     public final void get(Matrix4d m1)
 913     {
 914         if (nRow < 4 || nCol < 4) {
 915             m1.setZero();
 916             if (nCol > 0) {
 917                 if (nRow > 0) {
 918                     m1.m00 = values[0][0];
 919                     if (nRow > 1) {
 920                         m1.m10 = values[1][0];
 921                         if (nRow > 2) {
 922                             m1.m20 = values[2][0];
 923                             if (nRow > 3) {
 924                                 m1.m30 = values[3][0];
 925                             }
 926                         }
 927                     }
 928                 }
 929                 if (nCol > 1) {
 930                     if (nRow > 0) {
 931                         m1.m01 = values[0][1];
 932                         if (nRow > 1) {
 933                             m1.m11 = values[1][1];
 934                             if (nRow > 2) {
 935                                 m1.m21 = values[2][1];
 936                                 if (nRow > 3) {
 937                                     m1.m31 = values[3][1];
 938                                 }
 939                             }
 940                         }
 941                     }
 942                     if (nCol > 2) {
 943                         if (nRow > 0) {
 944                             m1.m02 = values[0][2];
 945                             if (nRow > 1) {
 946                                 m1.m12 = values[1][2];
 947                                 if (nRow > 2) {
 948                                     m1.m22 = values[2][2];
 949                                     if (nRow > 3) {
 950                                         m1.m32 = values[3][2];
 951                                     }
 952                                 }
 953                             }
 954                         }
 955                         if (nCol > 3) {
 956                             if (nRow > 0) {
 957                                 m1.m03 = values[0][3];
 958                                 if (nRow > 1) {
 959                                     m1.m13 = values[1][3];
 960                                     if (nRow > 2) {
 961                                         m1.m23 = values[2][3];
 962                                         if (nRow > 3) {
 963                                             m1.m33 = values[3][3];
 964                                         }
 965                                     }
 966                                 }
 967                             }
 968                         }
 969                     }
 970                 }
 971             }
 972         } else {
 973             m1.m00 = values[0][0];
 974             m1.m01 = values[0][1];
 975             m1.m02 = values[0][2];
 976             m1.m03 = values[0][3];
 977
 978             m1.m10 = values[1][0];
 979             m1.m11 = values[1][1];
 980             m1.m12 = values[1][2];
 981             m1.m13 = values[1][3];
 982
 983             m1.m20 = values[2][0];
 984             m1.m21 = values[2][1];
 985             m1.m22 = values[2][2];
 986             m1.m23 = values[2][3];
 987
 988             m1.m30 = values[3][0];
 989             m1.m31 = values[3][1];
 990             m1.m32 = values[3][2];
 991             m1.m33 = values[3][3];
 992         }
 993
 994     }
 995
 996     /**
 997      * Places the values in the upper 4x4 of this GMatrix into
 998      * the matrix m1.
 999      * @param m1  The matrix that will hold the new values
1000      */
1001     public final void get(Matrix4f m1)
1002     {
1003
1004         if (nRow < 4 || nCol < 4) {
1005             m1.setZero();
1006             if (nCol > 0) {
1007                 if (nRow > 0) {
1008                     m1.m00 = (float)values[0][0];
1009                     if (nRow > 1) {
1010                         m1.m10 = (float)values[1][0];
1011                         if (nRow > 2) {
1012                             m1.m20 = (float)values[2][0];
1013                             if (nRow > 3) {
1014                                 m1.m30 = (float)values[3][0];
1015                             }
1016                         }
1017                     }
1018                 }
1019                 if (nCol > 1) {
1020                     if (nRow > 0) {
1021                         m1.m01 = (float)values[0][1];
1022                         if (nRow > 1) {
1023                             m1.m11 = (float)values[1][1];
1024                             if (nRow > 2) {
1025                                 m1.m21 = (float)values[2][1];
1026                                 if (nRow > 3) {
1027                                     m1.m31 = (float)values[3][1];
1028                                 }
1029                             }
1030                         }
1031                     }
1032                     if (nCol > 2) {
1033                         if (nRow > 0) {
1034                             m1.m02 = (float)values[0][2];
1035                             if (nRow > 1) {
1036                                 m1.m12 = (float)values[1][2];
1037                                 if (nRow > 2) {
1038                                     m1.m22 = (float)values[2][2];
1039                                     if (nRow > 3) {
1040                                         m1.m32 = (float)values[3][2];
1041                                     }
1042                                 }
1043                             }
1044                         }
1045                         if (nCol > 3) {
1046                             if (nRow > 0) {
1047                                 m1.m03 = (float)values[0][3];
1048                                 if (nRow > 1) {
1049                                     m1.m13 = (float)values[1][3];
1050                                     if (nRow > 2) {
1051                                         m1.m23 = (float)values[2][3];
1052                                         if (nRow > 3) {
1053                                             m1.m33 = (float)values[3][3];
1054                                         }
1055                                     }
1056                                 }
1057                             }
1058                         }
1059                     }
1060                 }
1061             }
1062         } else {
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];
1067
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];
1072
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];
1077
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];
1082         }
1083     }
1084
1085     /**
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
1089      */
1090     public final void get(GMatrix m1)
1091     {
1092         int i, j, nc, nr;
1093
1094         if (nCol < m1.nCol)
1095             nc = nCol;
1096         else
1097             nc = m1.nCol;
1098
1099         if (nRow < m1.nRow)
1100             nr = nRow;
1101         else
1102             nr = m1.nRow;
1103
1104         for (i = 0; i < nr; i++) {
1105             for (j = 0; j < nc; j++) {
1106                 m1.values[i][j] = values[i][j];
1107             }
1108         }
1109         for (i = nr; i < m1.nRow; i++) {
1110             for (j = 0; j < m1.nCol; j++) {
1111                 m1.values[i][j] = 0.0;
1112             }
1113         }
1114         for (j = nc; j < m1.nCol; j++) {
1115             for (i = 0; i < nr; i++) {
1116                 m1.values[i][j] = 0.0;
1117             }
1118         }
1119     }
1120
1121     /**
1122      * Copy the values from the array into the specified row of this
1123      * matrix.
1124      * @param row  the row of this matrix into which the array values
1125      *             will be copied.
1126      * @param array  the source array
1127      */
1128     public final void setRow(int row, double[] array)
1129     {
1130         for (int i = 0; i < nCol; i++) {
1131             values[row][i] = array[i];
1132         }
1133     }
1134
1135     /**
1136      * Copy the values from the vector into the specified row of this
1137      * matrix.
1138      * @param row  the row of this matrix into which the array values
1139      *             will be copied
1140      * @param vector  the source vector
1141      */
1142     public final void setRow(int row, GVector vector)
1143     {
1144         for(int i = 0; i < nCol; i++) {
1145             values[row][i] = vector.values[i];
1146         }
1147     }
1148
1149     /**
1150      * Copy the values from the array into the specified column of this
1151      * matrix.
1152      * @param col  the column of this matrix into which the array values
1153      *             will be copied
1154      * @param array  the source array
1155      */
1156     public final void setColumn(int col, double[] array)
1157     {
1158         for(int i = 0; i < nRow; i++) {
1159             values[i][col] = array[i];
1160         }
1161     }
1162
1163     /**
1164      * Copy the values from the vector into the specified column of this
1165      * matrix.
1166      * @param col  the column of this matrix into which the array values
1167      *             will be copied
1168      * @param vector  the source vector
1169      */
1170     public final void setColumn(int col, GVector vector)
1171     {
1172         for(int i = 0; i < nRow; i++) {
1173             values[i][col] = vector.values[i];
1174         }
1175
1176     }
1177
1178     /**
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
1183      */
1184     public final void mulTransposeBoth(GMatrix m1, GMatrix m2)
1185     {
1186         int i, j, k;
1187
1188         if (m1.nRow != m2.nCol || nRow != m1.nCol || nCol != m2.nRow)
1189             throw new MismatchedSizeException
1190                 (VecMathI18N.getString("GMatrix14"));
1191
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++) {
1196                     tmp[i][j] = 0.0;
1197                     for (k = 0; k < m1.nRow; k++) {
1198                         tmp[i][j] += m1.values[k][i]*m2.values[j][k];
1199                     }
1200                 }
1201             }
1202             values = tmp;
1203         } else {
1204             for (i = 0; i < nRow; i++) {
1205                 for (j = 0; j < nCol; j++) {
1206                     values[i][j] = 0.0;
1207                     for (k = 0; k < m1.nRow; k++) {
1208                         values[i][j] += m1.values[k][i]*m2.values[j][k];
1209                     }
1210                 }
1211             }
1212         }
1213     }
1214
1215     /**
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
1220      */
1221     public final void mulTransposeRight(GMatrix m1, GMatrix m2)
1222     {
1223         int i, j, k;
1224
1225         if (m1.nCol != m2.nCol || nCol != m2.nRow || nRow != m1.nRow)
1226             throw new MismatchedSizeException
1227                 (VecMathI18N.getString("GMatrix15"));
1228
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++) {
1233                     tmp[i][j] = 0.0;
1234                     for (k = 0; k < m1.nCol; k++) {
1235                         tmp[i][j] += m1.values[i][k]*m2.values[j][k];
1236                     }
1237                 }
1238             }
1239             values = tmp;
1240         } else {
1241             for (i = 0; i < nRow; i++) {
1242                 for (j = 0;j < nCol; j++) {
1243                     values[i][j] = 0.0;
1244                     for (k = 0; k < m1.nCol; k++) {
1245                         values[i][j] += m1.values[i][k]*m2.values[j][k];
1246                     }
1247                 }
1248             }
1249         }
1250
1251     }
1252
1253
1254     /**
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
1259      */
1260     public final void mulTransposeLeft(GMatrix m1, GMatrix m2)
1261     {
1262         int i, j, k;
1263
1264         if (m1.nRow != m2.nRow || nCol != m2.nCol || nRow != m1.nCol)
1265             throw new MismatchedSizeException
1266                 (VecMathI18N.getString("GMatrix16"));
1267
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++) {
1272                     tmp[i][j] = 0.0;
1273                     for (k = 0; k < m1.nRow; k++) {
1274                         tmp[i][j] += m1.values[k][i]*m2.values[k][j];
1275                     }
1276                 }
1277             }
1278             values = tmp;
1279         } else {
1280             for (i = 0; i < nRow; i++) {
1281                 for (j = 0; j < nCol; j++) {
1282                     values[i][j] = 0.0;
1283                     for (k = 0; k < m1.nRow; k++) {
1284                         values[i][j] += m1.values[k][i]*m2.values[k][j];
1285                     }
1286                 }
1287             }
1288         }
1289     }
1290
1291
1292     /**
1293      * Transposes this matrix in place.
1294      */
1295     public final void transpose()
1296     {
1297         int i, j;
1298
1299         if (nRow != nCol) {
1300             double[][] tmp;
1301             i=nRow;
1302             nRow = nCol;
1303             nCol = i;
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];
1308                 }
1309             }
1310             values = tmp;
1311         } else {
1312             double swap;
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;
1318                 }
1319             }
1320         }
1321     }
1322
1323     /**
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)
1326      */
1327     public final void transpose(GMatrix m1)
1328     {
1329         int i, j;
1330
1331         if (nRow != m1.nCol || nCol != m1.nRow)
1332             throw new MismatchedSizeException
1333                 (VecMathI18N.getString("GMatrix17"));
1334
1335         if (m1 != this) {
1336             for (i = 0; i < nRow; i++) {
1337                 for (j = 0;j < nCol; j++) {
1338                     values[i][j] = m1.values[j][i];
1339                 }
1340             }
1341         } else {
1342             transpose();
1343         }
1344     }
1345
1346     /**
1347      * Returns a string that contains the values of this GMatrix.
1348      * @return the String representation
1349      */
1350     public String toString()
1351     {
1352         StringBuffer buffer = new StringBuffer(nRow*nCol*8);
1353
1354         int i, j;
1355
1356         for (i = 0; i < nRow; i++) {
1357             for (j = 0; j < nCol; j++) {
1358                 buffer.append(values[i][j]).append(" ");
1359             }
1360             buffer.append("\n");
1361         }
1362
1363         return buffer.toString();
1364     }
1365
1366     private static void checkMatrix( GMatrix m)
1367     {
1368         int i, j;
1369
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     ");
1374                 } else {
1375                     System.out.print(" " + m.values[i][j]);
1376                 }
1377             }
1378             System.out.print("\n");
1379         }
1380     }
1381
1382
1383     /**
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
1389      * likely.
1390      * @return the integer hash code value
1391      */
1392     public int hashCode() {
1393         long bits = 1L;
1394
1395         bits = 31L * bits + (long)nRow;
1396         bits = 31L * bits + (long)nCol;
1397
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]);
1401             }
1402         }
1403
1404         return (int) (bits ^ (bits >> 32));
1405     }
1406
1407
1408     /**
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
1413      */
1414     public boolean equals(GMatrix m1)
1415     {
1416         try {
1417             int i, j;
1418
1419             if (nRow != m1.nRow || nCol != m1.nCol)
1420                 return false;
1421
1422             for (i = 0;i < nRow; i++) {
1423                 for (j = 0; j < nCol; j++) {
1424                     if (values[i][j] != m1.values[i][j])
1425                         return false;
1426                 }
1427             }
1428             return true;
1429         }
1430         catch (NullPointerException e2) {
1431             return false;
1432         }
1433     }
1434
1435     /**
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
1438      * this GMatrix.
1439      * @param o1  The object with which the comparison is made.
1440      * @return  true or false
1441      */
1442     public boolean equals(Object o1)
1443     {
1444         try {
1445             GMatrix m2 = (GMatrix) o1;
1446             int i, j;
1447             if (nRow != m2.nRow || nCol != m2.nCol)
1448                 return false;
1449
1450             for (i = 0; i < nRow; i++) {
1451                 for (j = 0; j < nCol; j++) {
1452                     if (values[i][j] != m2.values[i][j])
1453                         return false;
1454                 }
1455             }
1456             return true;
1457         }
1458         catch (ClassCastException e1) {
1459             return false;
1460         }
1461         catch (NullPointerException e2) {
1462             return false;
1463         }
1464     }
1465
1466     /**
1467      * @deprecated Use epsilonEquals(GMatrix, double) instead
1468      */
1469     public boolean epsilonEquals(GMatrix m1, float epsilon) {
1470         return epsilonEquals(m1, (double)epsilon);
1471     }
1472
1473     /**
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
1481      */
1482     public boolean epsilonEquals(GMatrix m1, double epsilon)
1483     {
1484         int i, j;
1485         double diff;
1486         if (nRow != m1.nRow || nCol != m1.nCol)
1487             return false;
1488
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)
1493                     return false;
1494             }
1495         }
1496         return true;
1497     }
1498
1499     /**
1500      * Returns the trace of this matrix.
1501      * @return  the trace of this matrix
1502      */
1503     public final double trace()
1504     {
1505         int i, l;
1506         double t;
1507
1508         if (nRow < nCol)
1509             l = nRow;
1510         else
1511             l = nCol;
1512
1513         t = 0.0;
1514         for (i = 0; i < l; i++) {
1515             t += values[i][i];
1516         }
1517         return t;
1518     }
1519
1520     /**
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.
1535      */
1536     public final int SVD(GMatrix U, GMatrix W, GMatrix V)
1537     {
1538         // check for consistancy in dimensions
1539         if (nCol != V.nCol || nCol != V.nRow) {
1540             throw new MismatchedSizeException
1541                 (VecMathI18N.getString("GMatrix18"));
1542         }
1543
1544         if (nRow != U.nRow || nRow != U.nCol) {
1545             throw new MismatchedSizeException
1546                 (VecMathI18N.getString("GMatrix25"));
1547         }
1548
1549         if (nRow != W.nRow || nCol != W.nCol) {
1550             throw new MismatchedSizeException
1551                 (VecMathI18N.getString("GMatrix26"));
1552         }
1553
1554         // Fix ArrayIndexOutOfBounds for 2x2 matrices, which partially
1555         // addresses bug 4348562 for J3D 1.2.1.
1556         //
1557         // Does *not* fix the following problems reported in 4348562,
1558         // which will wait for J3D 1.3:
1559         //
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) {
1566                 U.setIdentity();
1567                 V.setIdentity();
1568
1569                 if (values[0][1] == 0.0) {
1570                     return 2;
1571                 }
1572
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];
1578
1579                 single_values[0] = values[0][0];
1580                 single_values[1] = values[1][1];
1581
1582                 compute_2X2(values[0][0], values[0][1], values[1][1],
1583                             single_values, sinl, cosl, sinr, cosr, 0);
1584
1585                 update_u(0, U, cosl, sinl);
1586                 update_v(0, V, cosr, sinr);
1587
1588                 return 2;
1589             }
1590             // else call computeSVD() and check for 2x2 there
1591         }
1592
1593         return computeSVD(this, U, W, V);
1594     }
1595
1596     /**
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
1609      * will be placed.
1610      * @param permutation  The row permutation effected by the partial
1611      * pivoting
1612      * @return  +-1 depending on whether the number of row interchanges
1613      * was even or odd respectively
1614      */
1615     public final int LUD(GMatrix LU, GVector permutation)
1616     {
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];
1621         int i, j;
1622
1623         if (nRow != nCol) {
1624             throw new MismatchedSizeException
1625                 (VecMathI18N.getString("GMatrix19"));
1626         }
1627
1628         if (nRow != LU.nRow) {
1629             throw new MismatchedSizeException
1630                 (VecMathI18N.getString("GMatrix27"));
1631         }
1632
1633         if (nCol != LU.nCol) {
1634             throw new MismatchedSizeException
1635                 (VecMathI18N.getString("GMatrix27"));
1636         }
1637
1638         if (LU.nRow != permutation.getSize()) {
1639             throw new MismatchedSizeException
1640                 (VecMathI18N.getString("GMatrix20"));
1641         }
1642
1643         for (i = 0; i < nRow; i++) {
1644             for (j = 0; j < nCol; j++) {
1645                 temp[i*nCol+j] = values[i][j];
1646             }
1647         }
1648
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"));
1654         }
1655
1656         for (i = 0; i < nRow; i++) {
1657             for (j = 0; j < nCol; j++) {
1658                 LU.values[i][j] = temp[i*nCol+j];
1659             }
1660         }
1661
1662         for (i = 0; i < LU.nRow; i++){
1663             permutation.values[i] = (double)row_perm[i];
1664         }
1665
1666         return even_row_exchange[0];
1667     }
1668
1669     /**
1670      *  Sets this matrix to a uniform scale matrix; all of the
1671      *  values are reset.
1672      *  @param scale  The new scale value
1673      */
1674     public final void setScale(double scale)
1675     {
1676         int i, j, l;
1677
1678         if (nRow < nCol)
1679             l = nRow;
1680         else
1681             l = nCol;
1682
1683         for (i = 0; i < nRow; i++) {
1684             for (j = 0; j < nCol; j++) {
1685                 values[i][j] = 0.0;
1686             }
1687         }
1688
1689         for (i = 0; i < l; i++) {
1690             values[i][i] = scale;
1691         }
1692     }
1693
1694     /**
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.
1698      *
1699      * Also note that since this routine is slow anyway, we won't worry
1700      * about allocating a little bit of garbage.
1701      */
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];
1708         int i, j;
1709
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"));
1716         }
1717
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];
1722             }
1723         }
1724
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"));
1730         }
1731
1732         // Perform back substitution on the identity matrix
1733         for (i = 0; i < size; i++)
1734             result[i] = 0.0;
1735
1736         for (i = 0; i < nCol; i++)
1737             result[i+i*nCol] = 1.0;
1738
1739         luBacksubstitution(m1.nRow, temp, row_perm, result);
1740
1741         for (i = 0; i < nRow; i++) {
1742             for (j = 0; j < nCol; j++) {
1743                 values[i][j] =  result[i*nCol+j];
1744             }
1745         }
1746     }
1747
1748     /**
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.
1757      *
1758      * @return true if the matrix is nonsingular, or false otherwise.
1759      */
1760     //
1761     // Reference: Press, Flannery, Teukolsky, Vetterling,
1762     //        _Numerical_Recipes_in_C_, Cambridge University Press,
1763     //        1988, pp 40-45.
1764     //
1765     static boolean luDecomposition(int dim, double[] matrix0,
1766                                    int[] row_perm, int[] even_row_xchg) {
1767
1768         double row_scale[] = new double[dim];
1769
1770         // Determine implicit scaling information by looping over rows
1771         int i, j;
1772         int ptr, rs, mtx;
1773         double big, temp;
1774
1775         ptr = 0;
1776         rs = 0;
1777         even_row_xchg[0] = 1;
1778
1779         // For each row ...
1780         i = dim;
1781         while (i-- != 0) {
1782             big = 0.0;
1783
1784             // For each column, find the largest element in the row
1785             j = dim;
1786             while (j-- != 0) {
1787                 temp = matrix0[ptr++];
1788                 temp = Math.abs(temp);
1789                 if (temp > big) {
1790                     big = temp;
1791                 }
1792             }
1793
1794             // Is the matrix singular?
1795             if (big == 0.0) {
1796                 return false;
1797             }
1798             row_scale[rs++] = 1.0 / big;
1799         }
1800
1801         // For all columns, execute Crout's method
1802         mtx = 0;
1803         for (j = 0; j < dim; j++) {
1804             int imax, k;
1805             int target, p1, p2;
1806             double sum;
1807
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];
1812                 k = i;
1813                 p1 = mtx + (dim*i);
1814                 p2 = mtx + j;
1815                 while (k-- != 0) {
1816                     sum -= matrix0[p1] * matrix0[p2];
1817                     p1++;
1818                     p2 += dim;
1819                 }
1820                 matrix0[target] = sum;
1821             }
1822
1823             // Search for largest pivot element and calculate
1824             // intermediate elements of lower diagonal matrix L.
1825             big = 0.0;
1826             imax = -1;
1827             for (i = j; i < dim; i++) {
1828                 target = mtx + (dim*i) + j;
1829                 sum = matrix0[target];
1830                 k = j;
1831                 p1 = mtx + (dim*i);
1832                 p2 = mtx + j;
1833                 while (k-- != 0) {
1834                     sum -= matrix0[p1] * matrix0[p2];
1835                     p1++;
1836                     p2 += dim;
1837                 }
1838                 matrix0[target] = sum;
1839
1840                 // Is this the best pivot so far?
1841                 if ((temp = row_scale[i] * Math.abs(sum)) >= big) {
1842                     big = temp;
1843                     imax = i;
1844                 }
1845             }
1846
1847             if (imax < 0) {
1848                 throw new RuntimeException(VecMathI18N.getString("GMatrix24"));
1849             }
1850
1851             // Is a row exchange necessary?
1852             if (j != imax) {
1853                 // Yes: exchange rows
1854                 k = dim;
1855                 p1 = mtx + (dim*imax);
1856                 p2 = mtx + (dim*j);
1857                 while (k-- != 0) {
1858                     temp = matrix0[p1];
1859                     matrix0[p1++] = matrix0[p2];
1860                     matrix0[p2++] = temp;
1861                 }
1862
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
1866             }
1867
1868             // Record row permutation
1869             row_perm[j] = imax;
1870
1871             // Is the matrix singular
1872             if (matrix0[(mtx + (dim*j) + j)] == 0.0) {
1873                 return false;
1874             }
1875
1876             // Divide elements of lower diagonal matrix L by pivot
1877             if (j != (dim-1)) {
1878                 temp = 1.0 / (matrix0[(mtx + (dim*j) + j)]);
1879                 target = mtx + (dim*(j+1)) + j;
1880                 i = (dim-1) - j;
1881                 while (i-- != 0) {
1882                     matrix0[target] *= temp;
1883                     target += dim;
1884                 }
1885             }
1886
1887         }
1888
1889         return true;
1890     }
1891
1892     /**
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.
1900      *
1901      * If "matrix2" is the identity matrix, the procedure replaces its contents
1902      * with the inverse of the matrix from which "matrix1" was originally
1903      * derived.
1904      */
1905     //
1906     // Reference: Press, Flannery, Teukolsky, Vetterling,
1907     //        _Numerical_Recipes_in_C_, Cambridge University Press,
1908     //        1988, pp 44-45.
1909     //
1910     static void luBacksubstitution(int dim, double[] matrix1,
1911                                    int[] row_perm,
1912                                    double[] matrix2) {
1913
1914         int i, ii, ip, j, k;
1915         int rp;
1916         int cv, rv, ri;
1917         double tt;
1918
1919         // rp = row_perm;
1920         rp = 0;
1921
1922         // For each column vector of matrix2 ...
1923         for (k = 0; k < dim; k++) {
1924             // cv = &(matrix2[0][k]);
1925             cv = k;
1926             ii = -1;
1927
1928             // Forward substitution
1929             for (i = 0; i < dim; i++) {
1930                 double sum;
1931
1932                 ip = row_perm[rp+i];
1933                 sum = matrix2[cv+dim*ip];
1934                 matrix2[cv+dim*ip] = matrix2[cv+dim*i];
1935                 if (ii >= 0) {
1936                     // rv = &(matrix1[i][0]);
1937                     rv = i*dim;
1938                     for (j = ii; j <= i-1; j++) {
1939                         sum -= matrix1[rv+j] * matrix2[cv+dim*j];
1940                     }
1941                 }
1942                 else if (sum != 0.0) {
1943                     ii = i;
1944                 }
1945                 matrix2[cv+dim*i] = sum;
1946             }
1947
1948             // Backsubstitution
1949             for (i = 0; i < dim; i++) {
1950                 ri = (dim-1-i);
1951                 rv = dim*(ri);
1952                 tt = 0.0;
1953                 for(j=1;j<=i;j++) {
1954                     tt += matrix1[rv+dim-j] * matrix2[cv+dim*(dim-j)];
1955                 }
1956                 matrix2[cv+dim*ri]= (matrix2[cv+dim*ri] - tt) / matrix1[rv+ri];
1957             }
1958         }
1959     }
1960
1961     static int computeSVD(GMatrix mat, GMatrix U, GMatrix W, GMatrix V) {
1962         int i, j, k;
1963         int nr, nc, si;
1964
1965         int converged, rank;
1966         double cs, sn, r, mag,scale, t;
1967         int eLength, sLength, vecLength;
1968
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);
1973
1974         // compute the number of singular values
1975         if (m.nRow >= m.nCol) {
1976             sLength = m.nCol;
1977             eLength = m.nCol-1;
1978         }else {
1979             sLength = m.nRow;
1980             eLength = m.nRow;
1981         }
1982
1983         if (m.nRow > m.nCol)
1984             vecLength = m.nRow;
1985         else
1986             vecLength = m.nCol;
1987
1988         double[] vec = new double[vecLength];
1989         double[] single_values = new double[sLength];
1990         double[] e = new double[eLength];
1991
1992         if(debug) {
1993             System.out.println("input to compute_svd = \n"+m.toString());
1994         }
1995
1996         rank = 0;
1997
1998         U.setIdentity();
1999         V.setIdentity();
2000
2001         nr = m.nRow;
2002         nc = m.nCol;
2003
2004         // householder reduction
2005         for (si = 0; si < sLength; si++) {
2006             // for each singular value
2007
2008             if (nr > 1) {
2009                 // zero out column
2010                 if (debug)
2011                     System.out.println
2012                         ("*********************** U ***********************\n");
2013
2014                 // compute reflector
2015                 mag = 0.0;
2016                 for (i = 0; i < nr; i++) {
2017                     mag += m.values[i+si][si] * m.values[i+si][si];
2018                     if (debug)
2019                         System.out.println
2020                             ("mag = " + mag + " matrix.dot = " +
2021                              m.values[i+si][si] * m.values[i+si][si]);
2022                 }
2023
2024                 mag = Math.sqrt(mag);
2025                 if (m.values[si][si] == 0.0) {
2026                     vec[0] = mag;
2027                 } else {
2028                     vec[0] = m.values[si][si] + d_sign(mag, m.values[si][si]);
2029                 }
2030
2031                 for (i = 1; i < nr; i++) {
2032                     vec[i] =  m.values[si+i][si];
2033                 }
2034
2035                 scale = 0.0;
2036                 for (i = 0; i < nr; i++) {
2037                     if (debug)
2038                         System.out.println("vec["+i+"]="+vec[i]);
2039
2040                     scale += vec[i]*vec[i];
2041                 }
2042
2043                 scale = 2.0/scale;
2044                 if (debug)
2045                     System.out.println("scale = "+scale);
2046
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];
2050                     }
2051                 }
2052
2053                 for (i = si; i < m.nRow; i++){
2054                     u.values[i][i] +=  1.0;
2055                 }
2056
2057                 // compute s
2058                 t = 0.0;
2059                 for (i = si; i < m.nRow; i++){
2060                     t += u.values[si][i] * m.values[i][si];
2061                 }
2062                 m.values[si][si] = t;
2063
2064                 // apply reflector
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];
2070                         }
2071                     }
2072                 }
2073
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];
2077                     }
2078                 }
2079
2080                 if (debug) {
2081                     System.out.println("U =\n" + U.toString());
2082                     System.out.println("u =\n" + u.toString());
2083                 }
2084
2085                 // update U matrix
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];
2091                         }
2092                     }
2093                 }
2094
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];
2098                     }
2099                 }
2100
2101                 if (debug) {
2102                     System.out.println("single_values["+si+"] =\n" +
2103                                        single_values[si]);
2104                     System.out.println("m =\n" + m.toString());
2105                     System.out.println("U =\n" + U.toString());
2106                 }
2107
2108                 nr--;
2109             }
2110
2111             if( nc > 2 ) {
2112                 // zero out row
2113                 if (debug)
2114                     System.out.println
2115                         ("*********************** V ***********************\n");
2116
2117                 mag = 0.0;
2118                 for (i = 1; i < nc; i++){
2119                     mag += m.values[si][si+i] * m.values[si][si+i];
2120                 }
2121
2122                 if (debug)
2123                     System.out.println("mag = " + mag);
2124
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) {
2129                     vec[0] = mag;
2130                 } else {
2131                     vec[0] = m.values[si][si+1] +
2132                         d_sign(mag, m.values[si][si+1]);
2133                 }
2134
2135                 for (i = 1; i < nc - 1; i++){
2136                     vec[i] =  m.values[si][si+i+1];
2137                 }
2138
2139                 // use reflection vector to compute v matrix
2140                 scale = 0.0;
2141                 for (i = 0; i < nc - 1; i++){
2142                     if( debug )System.out.println("vec["+i+"]="+vec[i]);
2143                     scale += vec[i]*vec[i];
2144                 }
2145
2146                 scale = 2.0/scale;
2147                 if (debug)
2148                     System.out.println("scale = "+scale);
2149
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];
2153                     }
2154                 }
2155
2156                 for (i = si + 1; i < m.nCol; i++){
2157                     v.values[i][i] +=  1.0;
2158                 }
2159
2160                 t=0.0;
2161                 for (i = si; i < m.nCol; i++){
2162                     t += v.values[i][si+1] * m.values[si][i];
2163                 }
2164                 m.values[si][si+1]=t;
2165
2166                 // apply reflector
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];
2172                         }
2173                     }
2174                 }
2175
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];
2179                     }
2180                 }
2181
2182                 if (debug) {
2183                     System.out.println("V =\n" + V.toString());
2184                     System.out.println("v =\n" + v.toString());
2185                     System.out.println("tmp =\n" + tmp.toString());
2186                 }
2187
2188                 // update V matrix
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];
2194                         }
2195                     }
2196                 }
2197
2198                 if (debug)
2199                     System.out.println("tmp =\n" + tmp.toString());
2200
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];
2204                     }
2205                 }
2206
2207                 if (debug) {
2208                     System.out.println("m =\n" + m.toString());
2209                     System.out.println("V =\n" + V.toString());
2210                 }
2211
2212                 nc--;
2213             }
2214         }
2215
2216         for (i = 0; i < sLength; i++){
2217             single_values[i] = m.values[i][i];
2218         }
2219
2220         for (i = 0; i < eLength; i++){
2221             e[i] = m.values[i][i+1];
2222         }
2223
2224         // Fix ArrayIndexOutOfBounds for 2x2 matrices, which partially
2225         // addresses bug 4348562 for J3D 1.2.1.
2226         //
2227         // Does *not* fix the following problems reported in 4348562,
2228         // which will wait for J3D 1.3:
2229         //
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];
2239
2240             compute_2X2(single_values[0], e[0], single_values[1],
2241                         single_values, sinl, cosl, sinr, cosr, 0);
2242
2243             update_u(0, U, cosl, sinl);
2244             update_v(0, V, cosr, sinr);
2245
2246             return 2;
2247         }
2248
2249         // compute_qr causes ArrayIndexOutOfBounds for 2x2 matrices
2250         compute_qr (0, e.length-1, single_values, e, U, V);
2251
2252         // compute rank = number of non zero singular values
2253         rank = single_values.length;
2254
2255         // sort by order of size of single values
2256         // and check for zero's
2257         return rank;
2258     }
2259
2260     static void compute_qr(int start, int end, double[] s, double[] e,
2261                            GMatrix u, GMatrix v) {
2262
2263         int i, j, k, n, sl;
2264         boolean converged;
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);
2271
2272         final int MAX_INTERATIONS = 2;
2273         final double CONVERGE_TOL = 4.89E-15;
2274
2275         if (debug) {
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]);
2280             }
2281
2282             System.out.println("\nes =\n");
2283             for (i = 0; i < e.length; i++) {
2284                 System.out.println(e[i]);
2285             }
2286
2287             for (i = 0; i < s.length; i++) {
2288                 m.values[i][i] = s[i];
2289             }
2290
2291             for (i = 0; i < e.length; i++) {
2292                 m.values[i][i+1] = e[i];
2293             }
2294             System.out.println("\nm =\n" + m.toString());
2295         }
2296
2297         double c_b48 =  1.0;
2298         double c_b71 = -1.0;
2299         converged = false;
2300
2301         if (debug)
2302             print_svd(s, e, u, v);
2303
2304         f = 0.0;
2305         g = 0.0;
2306
2307         for (k = 0; k < MAX_INTERATIONS && !converged;k++) {
2308             for (i = start; i <= end; i++) {
2309
2310                 // if at start of iterfaction compute shift
2311                 if (i == start) {
2312                     if (e.length == s.length)
2313                         sl = end;
2314                     else
2315                         sl = end + 1;
2316
2317                     shift = compute_shift(s[sl-1], e[end], s[sl]);
2318
2319                     f = (Math.abs(s[i]) - shift) *
2320                         (d_sign(c_b48, s[i]) + shift/s[i]);
2321                     g = e[i];
2322                 }
2323
2324                 r = compute_rot(f, g, sinr, cosr);
2325                 if (i != start)
2326                     e[i-1] = r;
2327
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];
2332
2333                 // if (debug) print_se(s,e);
2334                 update_v (i, v, cosr, sinr);
2335                 if (debug)
2336                     print_m(m,u,v);
2337
2338                 r = compute_rot(f, g, sinl, cosl);
2339                 s[i] = r;
2340                 f = cosl[0] * e[i] + sinl[0] * s[i+1];
2341                 s[i+1] = cosl[0] * s[i+1] - sinl[0] * e[i];
2342
2343                 if( i < end) {
2344                     // if not last
2345                     g = sinl[0] * e[i+1];
2346                     e[i+1] =  cosl[0] * e[i+1];
2347                 }
2348                 //if (debug) print_se(s,e);
2349
2350                 update_u(i, u, cosl, sinl);
2351                 if (debug)
2352                     print_m(m,u,v);
2353             }
2354
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];
2361
2362                 update_v(i, v, cosr, sinr);
2363                 if (debug)
2364                     print_m(m,u,v);
2365             }
2366
2367             if (debug) {
2368                 System.out.println
2369                     ("\n*********************** iteration #" + k +
2370                      " ***********************\n");
2371                 print_svd(s, e, u, v);
2372             }
2373
2374             // check for convergence on off diagonals and reduce
2375             while ((end-start > 1) && (Math.abs(e[end]) < CONVERGE_TOL)) {
2376                 end--;
2377             }
2378
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
2384
2385                     // check for convergence on off diagonals and reduce
2386                     while ((end - start > 1) &&
2387                            (Math.abs(e[end]) < CONVERGE_TOL)) {
2388                         end--;
2389                     }
2390                 }
2391             }
2392
2393             if (debug)
2394                 System.out.println("start = " + start);
2395
2396             if ((end - start <= 1) && (Math.abs(e[start+1]) < CONVERGE_TOL)) {
2397                 converged = true;
2398             } else {
2399                 // check if zero on the diagonal
2400             }
2401
2402         }
2403
2404         if (debug)
2405             System.out.println("\n****call compute_2X2 ********************\n");
2406
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);
2410             e[start] = 0.0;
2411             e[start+1] = 0.0;
2412         } else {
2413         }
2414
2415         i = start;
2416         update_u(i, u, cosl, sinl);
2417         update_v(i, v, cosr, sinr);
2418
2419         if(debug) {
2420             System.out.println
2421                 ("\n*******after call compute_2X2 **********************\n");
2422             print_svd(s, e, u, v);
2423         }
2424
2425         return;
2426     }
2427
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]);
2431     }
2432
2433     private static void update_v(int index, GMatrix v,
2434                                  double[] cosr, double[] sinr) {
2435         int j;
2436         double vtemp;
2437
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];
2444         }
2445     }
2446
2447     private static void chase_up(double[] s, double[] e, int k, GMatrix v) {
2448         double f, g, r;
2449         double[] cosr = new double[1];
2450         double[] sinr = new double[1];
2451         int i;
2452         GMatrix t = new GMatrix(v.nRow, v.nCol);
2453         GMatrix m = new GMatrix(v.nRow, v.nCol);
2454
2455         if (debug) {
2456             m.setIdentity();
2457             for (i = 0; i < s.length; i++) {
2458                 m.values[i][i] = s[i];
2459             }
2460             for (i = 0; i < e.length; i++) {
2461                 m.values[i][i+1] = e[i];
2462             }
2463         }
2464
2465         f = e[k];
2466         g = s[k];
2467
2468         for (i = k; i > 0; i--) {
2469             r = compute_rot(f, g, sinr, cosr);
2470             f = -e[i-1] * sinr[0];
2471             g = s[i-1];
2472             s[i] = r;
2473             e[i-1] = e[i-1] * cosr[0];
2474             update_v_split(i, k+1, v, cosr, sinr, t, m);
2475         }
2476
2477         s[i+1] = compute_rot(f, g, sinr, cosr);
2478         update_v_split(i, k+1, v, cosr, sinr, t, m);
2479     }
2480
2481     private static void chase_across(double[] s, double[] e, int k, GMatrix u) {
2482         double f, g, r;
2483         double[] cosl = new double[1];
2484         double[] sinl = new double[1];
2485         int i;
2486         GMatrix t = new GMatrix(u.nRow, u.nCol);
2487         GMatrix m = new GMatrix(u.nRow, u.nCol);
2488
2489         if (debug) {
2490             m.setIdentity();
2491             for (i = 0; i < s.length; i++) {
2492                 m.values[i][i] = s[i];
2493             }
2494             for (i = 0; i < e.length; i++) {
2495                 m.values[i][i+1] = e[i];
2496             }
2497         }
2498
2499         g = e[k];
2500         f = s[k+1];
2501
2502         for (i = k; i < u.nCol-2; i++){
2503             r = compute_rot(f, g, sinl, cosl);
2504             g = -e[i+1] * sinl[0];
2505             f = s[i+2];
2506             s[i+1] = r;
2507             e[i+1] = e[i+1] * cosl[0];
2508             update_u_split(k, i + 1, u, cosl, sinl, t, m);
2509         }
2510
2511         s[i+1] = compute_rot(f, g, sinl, cosl);
2512         update_u_split(k, i + 1, u, cosl, sinl, t, m);
2513     }
2514
2515     private static void update_v_split(int topr, int bottomr, GMatrix v,
2516                                        double[] cosr, double[] sinr,
2517                                        GMatrix t, GMatrix m) {
2518         int j;
2519         double vtemp;
2520
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];
2525         }
2526
2527         if (debug) {
2528             t.setIdentity();
2529             for (j = 0; j < v.nRow; j++) {
2530                 vtemp = t.values[j][topr];
2531                 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];
2535             }
2536         }
2537
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 =");
2543         checkMatrix(m);
2544         System.out.println("\nv =");
2545         checkMatrix(t);
2546         m.mul(m,t);
2547         System.out.println("\nt*m =");
2548         checkMatrix(m);
2549     }
2550
2551     private static void update_u_split(int topr, int bottomr, GMatrix u,
2552                                        double[] cosl, double[] sinl,
2553                                        GMatrix t, GMatrix m) {
2554         int j;
2555         double utemp;
2556
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];
2561         }
2562
2563         if(debug) {
2564             t.setIdentity();
2565             for (j = 0;j < u.nCol; j++) {
2566                 utemp = t.values[topr][j];
2567                 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];
2571             }
2572         }
2573         System.out.println("\nm=");
2574         checkMatrix(m);
2575         System.out.println("\nu=");
2576         checkMatrix(t);
2577         m.mul(t,m);
2578         System.out.println("\nt*m=");
2579         checkMatrix(m);
2580     }
2581
2582     private static void update_u(int index, GMatrix u,
2583                                  double[] cosl, double[] sinl) {
2584         int j;
2585         double utemp;
2586
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];
2593         }
2594     }
2595
2596     private static void print_m(GMatrix m, GMatrix u, GMatrix v) {
2597         GMatrix mtmp = new GMatrix(m.nCol, m.nRow);
2598
2599         mtmp.mul(u, mtmp);
2600         mtmp.mul(mtmp, v);
2601         System.out.println("\n m = \n" + mtmp.toString(mtmp));
2602
2603     }
2604
2605     private static String toString(GMatrix m)
2606     {
2607         StringBuffer buffer = new StringBuffer(m.nRow * m.nCol * 8);
2608         int i, j;
2609
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 ");
2614                 } else {
2615                     buffer.append(m.values[i][j]).append(" ");
2616                 }
2617             }
2618             buffer.append("\n");
2619         }
2620         return buffer.toString();
2621     }
2622
2623     private static void print_svd(double[] s, double[] e,
2624                                   GMatrix u, GMatrix v) {
2625         int i;
2626         GMatrix mtmp = new GMatrix(u.nCol, v.nRow);
2627
2628         System.out.println(" \ns = ");
2629         for (i = 0; i < s.length; i++) {
2630             System.out.println(" " + s[i]);
2631         }
2632
2633         System.out.println(" \ne = ");
2634         for (i = 0; i < e.length; i++) {
2635             System.out.println(" " + e[i]);
2636         }
2637
2638         System.out.println(" \nu  = \n" + u.toString());
2639         System.out.println(" \nv  = \n" + v.toString());
2640
2641         mtmp.setIdentity();
2642         for (i = 0; i < s.length; i++) {
2643             mtmp.values[i][i] = s[i];
2644         }
2645         for (i = 0; i < e.length; i++) {
2646             mtmp.values[i][i+1] = e[i];
2647         }
2648         System.out.println(" \nm  = \n"+mtmp.toString());
2649
2650         mtmp.mulTransposeLeft(u, mtmp);
2651         mtmp.mulTransposeRight(mtmp, v);
2652
2653         System.out.println(" \n u.transpose*m*v.transpose  = \n" +
2654                            mtmp.toString());
2655     }
2656
2657     static double max(double a, double b) {
2658         if (a > b)
2659             return a;
2660         else
2661             return b;
2662     }
2663
2664     static double min(double a, double b) {
2665         if (a < b)
2666             return a;
2667         else
2668             return b;
2669     }
2670
2671     static double compute_shift(double f, double g, double h) {
2672         double d__1, d__2;
2673         double fhmn, fhmx, c, fa, ga, ha, as, at, au;
2674         double ssmin;
2675
2676         fa = Math.abs(f);
2677         ga = Math.abs(g);
2678         ha = Math.abs(h);
2679         fhmn = min(fa,ha);
2680         fhmx = max(fa,ha);
2681
2682         if (fhmn == 0.0) {
2683             ssmin = 0.0;
2684             if (fhmx == 0.0) {
2685             } else {
2686                 d__1 = min(fhmx,ga) / max(fhmx,ga);
2687             }
2688         } else {
2689             if (ga < fhmx) {
2690                 as = fhmn / fhmx + 1.0;
2691                 at = (fhmx - fhmn) / fhmx;
2692                 d__1 = ga / fhmx;
2693                 au = d__1 * d__1;
2694                 c = 2.0 / (Math.sqrt(as * as + au) + Math.sqrt(at * at + au));
2695                 ssmin = fhmn * c;
2696             } else {
2697                 au = fhmx / ga;
2698                 if (au == 0.0) {
2699                     ssmin = fhmn * fhmx / ga;
2700                 } else {
2701                     as = fhmn / fhmx + 1.0;
2702                     at = (fhmx - fhmn) / fhmx;
2703                     d__1 = as * au;
2704                     d__2 = at * au;
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;
2708                     ssmin += ssmin;
2709                 }
2710             }
2711         }
2712
2713         return ssmin;
2714     }
2715
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) {
2719
2720         double c_b3 = 2.0;
2721         double c_b4 = 1.0;
2722
2723         double d__1;
2724         int pmax;
2725         double temp;
2726         boolean swap;
2727         double a, d, l, m, r, s, t, tsign, fa, ga, ha;
2728         double ft, gt, ht, mm;
2729         boolean gasmal;
2730         double tt, clt, crt, slt, srt;
2731         double ssmin,ssmax;
2732
2733         ssmax = single_values[0];
2734         ssmin = single_values[1];
2735         clt = 0.0;
2736         crt = 0.0;
2737         slt = 0.0;
2738         srt = 0.0;
2739         tsign = 0.0;
2740
2741         ft = f;
2742         fa = Math.abs(ft);
2743         ht = h;
2744         ha = Math.abs(h);
2745
2746         pmax = 1;
2747         if (ha > fa)
2748             swap = true;
2749         else
2750             swap = false;
2751
2752         if (swap) {
2753             pmax = 3;
2754             temp = ft;
2755             ft = ht;
2756             ht = temp;
2757             temp = fa;
2758             fa = ha;
2759             ha = temp;
2760
2761         }
2762
2763         gt = g;
2764         ga = Math.abs(gt);
2765         if (ga == 0.0) {
2766             single_values[1] = ha;
2767             single_values[0] = fa;
2768             clt = 1.0;
2769             crt = 1.0;
2770             slt = 0.0;
2771             srt = 0.0;
2772         } else {
2773             gasmal = true;
2774             if (ga > fa) {
2775                 pmax = 2;
2776                 if (fa / ga < EPS) {
2777                     gasmal = false;
2778                     ssmax = ga;
2779
2780                     if (ha > 1.0) {
2781                         ssmin = fa / (ga / ha);
2782                     } else {
2783                         ssmin = fa / ga * ha;
2784                     }
2785                     clt = 1.0;
2786                     slt = ht / gt;
2787                     srt = 1.0;
2788                     crt = ft / gt;
2789                 }
2790             }
2791             if (gasmal) {
2792                 d = fa - ha;
2793                 if (d == fa) {
2794
2795                     l = 1.0;
2796                 } else {
2797                     l = d / fa;
2798                 }
2799
2800                 m = gt / ft;
2801                 t = 2.0 - l;
2802                 mm = m * m;
2803                 tt = t * t;
2804                 s = Math.sqrt(tt + mm);
2805
2806                 if (l == 0.0) {
2807                     r = Math.abs(m);
2808                 } else {
2809                     r = Math.sqrt(l * l + mm);
2810                 }
2811
2812                 a = (s + r) * 0.5;
2813                 if (ga > fa) {
2814                     pmax = 2;
2815                     if (fa / ga < EPS) {
2816                         gasmal = false;
2817                         ssmax = ga;
2818                         if (ha > 1.0) {
2819                             ssmin = fa / (ga / ha);
2820                         } else {
2821                             ssmin = fa / ga * ha;
2822                         }
2823                         clt = 1.0;
2824                         slt = ht / gt;
2825                         srt = 1.0;
2826                         crt = ft / gt;
2827                     }
2828                 }
2829                 if (gasmal) {
2830                     d = fa - ha;
2831                     if (d == fa) {
2832                         l = 1.0;
2833                     } else {
2834                         l = d / fa;
2835                     }
2836
2837                     m = gt / ft;
2838                     t = 2.0 - l;
2839
2840                     mm = m * m;
2841                     tt = t * t;
2842                     s = Math.sqrt(tt + mm);
2843
2844                     if (l == 0.) {
2845                         r = Math.abs(m);
2846                     } else {
2847                         r = Math.sqrt(l * l + mm);
2848                     }
2849
2850                     a = (s + r) * 0.5;
2851                     ssmin = ha / a;
2852                     ssmax = fa * a;
2853
2854                     if (mm == 0.0) {
2855                         if (l == 0.0) {
2856                             t = d_sign(c_b3, ft) * d_sign(c_b4, gt);
2857                         } else {
2858                             t = gt / d_sign(d, ft) + m / t;
2859                         }
2860                     } else {
2861                         t = (m / (s + t) + m / (r + l)) * (a + 1.0);
2862                     }
2863
2864                     l = Math.sqrt(t * t + 4.0);
2865                     crt = 2.0 / l;
2866                     srt = t / l;
2867                     clt = (crt + srt * m) / a;
2868                     slt = ht / ft * srt / a;
2869                 }
2870             }
2871             if (swap) {
2872                 csl[0] = srt;
2873                 snl[0] = crt;
2874                 csr[0] = slt;
2875                 snr[0] = clt;
2876             } else {
2877                 csl[0] = clt;
2878                 snl[0] = slt;
2879                 csr[0] = crt;
2880                 snr[0] = srt;
2881             }
2882
2883             if (pmax == 1) {
2884                 tsign = d_sign(c_b4, csr[0]) *
2885                     d_sign(c_b4, csl[0]) * d_sign(c_b4, f);
2886             }
2887             if (pmax == 2) {
2888                 tsign = d_sign(c_b4, snr[0]) *
2889                     d_sign(c_b4, csl[0]) * d_sign(c_b4, g);
2890             }
2891             if (pmax == 3) {
2892                 tsign = d_sign(c_b4, snr[0]) *
2893                     d_sign(c_b4, snl[0]) * d_sign(c_b4, h);
2894             }
2895
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);
2899         }
2900
2901         return 0;
2902     }
2903
2904     static double compute_rot(double f, double g, double[] sin, double[] cos) {
2905         int i__1;
2906         double d__1, d__2;
2907         double cs, sn;
2908         int i;
2909         double scale;
2910         int count;
2911         double f1, g1;
2912         double r;
2913         final double safmn2 = 2.002083095183101E-146;
2914         final double safmx2 = 4.994797680505588E+145;
2915
2916         if (g == 0.0) {
2917             cs = 1.0;
2918             sn = 0.0;
2919             r = f;
2920         } else if (f == 0.0) {
2921             cs = 0.0;
2922             sn = 1.0;
2923             r = g;
2924         } else {
2925             f1 = f;
2926             g1 = g;
2927             scale = max(Math.abs(f1),Math.abs(g1));
2928             if (scale >= safmx2) {
2929                 count = 0;
2930                 while(scale >= safmx2) {
2931                     ++count;
2932                     f1 *= safmn2;
2933                     g1 *= safmn2;
2934                     scale = max(Math.abs(f1), Math.abs(g1));
2935                 }
2936                 r = Math.sqrt(f1*f1 + g1*g1);
2937                 cs = f1 / r;
2938                 sn = g1 / r;
2939                 i__1 = count;
2940                 for (i = 1; i <= count; ++i) {
2941                     r *= safmx2;
2942                 }
2943             } else if (scale <= safmn2) {
2944                 count = 0;
2945                 while(scale <= safmn2) {
2946                     ++count;
2947                     f1 *= safmx2;
2948                     g1 *= safmx2;
2949                     scale = max(Math.abs(f1), Math.abs(g1));
2950                 }
2951                 r = Math.sqrt(f1*f1 + g1*g1);
2952                 cs = f1 / r;
2953                 sn = g1 / r;
2954                 i__1 = count;
2955                 for (i = 1; i <= count; ++i) {
2956                     r *= safmn2;
2957                 }
2958             } else {
2959                 r = Math.sqrt(f1*f1 + g1*g1);
2960                 cs = f1 / r;
2961                 sn = g1 / r;
2962             }
2963             if (Math.abs(f) > Math.abs(g) && cs < 0.0) {
2964                 cs = -cs;
2965                 sn = -sn;
2966                 r = -r;
2967             }
2968         }
2969         sin[0] = sn;
2970         cos[0] = cs;
2971         return r;
2972     }
2973
2974     static double d_sign(double a, double b) {
2975         double x;
2976         x = (a >= 0 ? a : - a);
2977         return (b >= 0 ? x : -x);
2978     }
2979
2980     /**
2981      * Creates a new object of the same class as this object.
2982      *
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
2987      */
2988     public Object clone() {
2989         GMatrix m1 = null;
2990         try {
2991             m1 = (GMatrix)super.clone();
2992         } catch (CloneNotSupportedException e) {
2993             // this shouldn't happen, since we are Cloneable
2994             throw new InternalError();
2995         }
2996
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];
3002            }
3003         }
3004
3005         return m1;
3006     }
3007
3008 }