source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/math/Matrix.java

Last change on this file was 127, checked in by Wouter Pasman, 6 years ago

#41 ROLL BACK of rev.126 . So this version is equal to rev. 125

File size: 99.2 KB
Line 
1/* Class Matrix
2*
3* Defines a matrix and includes the methods needed
4* for standard matrix manipulations, e.g. multiplation,
5* and related procedures, e.g. solution of linear
6* simultaneous equations
7*
8* See class ComplexMatrix and PhasorMatrix for complex matrix arithmetic
9*
10* WRITTEN BY: Dr Michael Thomas Flanagan
11*
12* DATE: June 2002
13* UPDATES: 21 April 2004, 16 February 2006, 31 March 2006, 22 April 2006, 1 July 2007, 17 July 2007
14* 18 August 2007, 7 October 2007, 27 February 2008, 7 April 2008, 5 July 2008, 6-15 September 2008
15* 7-14 October 2008, 16 February 2009, 16 June 2009, 15 October 2009, 4-5 November 2009
16* 12 January 2010, 19 February 2010, 14 November 2010, 12 January 2011, 20 January 2011
17*
18*
19* DOCUMENTATION:
20* See Michael Thomas Flanagan's Java library on-line web page:
21* http://www.ee.ucl.ac.uk/~mflanaga/java/Matrix.html
22* http://www.ee.ucl.ac.uk/~mflanaga/java/
23*
24* Copyright (c) 2002 - 2011 Michael Thomas Flanagan
25* PERMISSION TO COPY:
26*
27* Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
28* provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
29* and associated documentation or publications.
30*
31* Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice, this list of conditions
32* and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
33*
34* Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
35* the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission from the Michael Thomas Flanagan:
36*
37* Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
38* Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
39* or its derivatives.
40*
41***************************************************************************************/
42
43package agents.anac.y2015.agentBuyogV2.flanagan.math;
44
45import java.util.ArrayList;
46import java.util.Vector;
47
48import agents.anac.y2015.agentBuyogV2.flanagan.analysis.Regression;
49import agents.anac.y2015.agentBuyogV2.flanagan.analysis.RegressionFunction;
50import agents.anac.y2015.agentBuyogV2.flanagan.analysis.Stat;
51import agents.anac.y2015.agentBuyogV2.flanagan.math.ArrayMaths;
52import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
53
54import java.math.BigDecimal;
55import java.math.BigInteger;
56
57public class Matrix{
58
59 private int numberOfRows = 0; // number of rows
60 private int numberOfColumns = 0; // number of columns
61 private double matrix[][] = null; // 2-D Matrix
62 private double hessenberg[][] = null; // 2-D Hessenberg equivalent
63 private boolean hessenbergDone = false; // = true when Hessenberg matrix calculated
64 private int permutationIndex[] = null; // row permutation index
65 private double rowSwapIndex = 1.0D; // row swap index
66 private double[] eigenValues = null; // eigen values of the matrix
67 private double[][] eigenVector = null; // eigen vectors of the matrix
68 private double[] sortedEigenValues = null; // eigen values of the matrix sorted into descending order
69 private double[][] sortedEigenVector = null; // eigen vectors of the matrix sorted to matching descending eigen value order
70 private int numberOfRotations = 0; // number of rotations in Jacobi transformation
71 private int[] eigenIndices = null; // indices of the eigen values before sorting into descending order
72 private int maximumJacobiIterations = 100; // maximum number of Jacobi iterations
73 private boolean eigenDone = false; // = true when eigen values and vectors calculated
74 private boolean matrixCheck = true; // check on matrix status
75 // true - no problems encountered in LU decomposition
76 // false - attempted a LU decomposition on a singular matrix
77
78 private boolean supressErrorMessage = false; // true - LU decompostion failure message supressed
79
80 private double tiny = 1.0e-100; // small number replacing zero in LU decomposition
81
82 // CONSTRUCTORS
83 // Construct a numberOfRows x numberOfColumns matrix of variables all equal to zero
84 public Matrix(int numberOfRows, int numberOfColumns){
85 this.numberOfRows = numberOfRows;
86 this.numberOfColumns = numberOfColumns;
87 this.matrix = new double[numberOfRows][numberOfColumns];
88 this.permutationIndex = new int[numberOfRows];
89 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
90 }
91
92 // Construct a numberOfRows x numberOfColumns matrix of variables all equal to the number const
93 public Matrix(int numberOfRows, int numberOfColumns, double constant){
94 this.numberOfRows = numberOfRows;
95 this.numberOfColumns = numberOfColumns;
96 this.matrix = new double[numberOfRows][numberOfColumns];
97 for(int i=0;i<numberOfRows;i++){
98 for(int j=0;j<numberOfColumns;j++)this.matrix[i][j]=constant;
99 }
100 this.permutationIndex = new int[numberOfRows];
101 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
102 }
103
104 // Construct matrix with a copy of an existing numberOfRows x numberOfColumns 2-D array of variables
105 public Matrix(double[][] twoD){
106 this.numberOfRows = twoD.length;
107 this.numberOfColumns = twoD[0].length;
108 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
109 for(int i=0; i<numberOfRows; i++){
110 if(twoD[i].length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
111 for(int j=0; j<numberOfColumns; j++){
112 this.matrix[i][j]=twoD[i][j];
113 }
114 }
115 this.permutationIndex = new int[numberOfRows];
116 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
117 }
118
119 // Construct matrix with a copy of an existing numberOfRows x numberOfColumns 2-D array of floats
120 public Matrix(float[][] twoD){
121 this.numberOfRows = twoD.length;
122 this.numberOfColumns = twoD[0].length;
123 for(int i=1; i<numberOfRows; i++){
124 if(twoD[i].length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
125 }
126 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
127 for(int i=0; i<numberOfRows; i++){
128 for(int j=0; j<numberOfColumns; j++){
129 this.matrix[i][j] = (double)twoD[i][j];
130 }
131 }
132 this.permutationIndex = new int[numberOfRows];
133 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
134 }
135
136 // Construct matrix with a copy of an existing numberOfRows x numberOfColumns 2-D array of longs
137 public Matrix(long[][] twoD){
138 this.numberOfRows = twoD.length;
139 this.numberOfColumns = twoD[0].length;
140 for(int i=1; i<numberOfRows; i++){
141 if(twoD[i].length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
142 }
143 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
144 for(int i=0; i<numberOfRows; i++){
145 for(int j=0; j<numberOfColumns; j++){
146 this.matrix[i][j] = (double)twoD[i][j];
147 }
148 }
149 this.permutationIndex = new int[numberOfRows];
150 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
151 }
152
153
154 // Construct matrix with a copy of an existing numberOfRows x numberOfColumns 2-D array of ints
155 public Matrix(int[][] twoD){
156 this.numberOfRows = twoD.length;
157 this.numberOfColumns = twoD[0].length;
158 for(int i=1; i<numberOfRows; i++){
159 if(twoD[i].length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
160 }
161 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
162 for(int i=0; i<numberOfRows; i++){
163 for(int j=0; j<numberOfColumns; j++){
164 this.matrix[i][j] = (double)twoD[i][j];
165 }
166 }
167 this.permutationIndex = new int[numberOfRows];
168 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
169 }
170
171 // Construct matrix with a copy of an existing numberOfRows 1-D array of ArrayMaths
172 public Matrix(ArrayMaths[] twoD){
173 this.numberOfRows = twoD.length;
174 this.numberOfColumns = twoD[0].length();
175 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
176 for(int i=0; i<numberOfRows; i++){
177 double[] arrayh = (twoD[i].copy()).array();
178 if(arrayh.length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
179 this.matrix[i] = arrayh;
180 }
181 this.permutationIndex = new int[numberOfRows];
182 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
183 }
184
185 // Construct matrix with a copy of an existing numberOfRows 1-D array of ArrayLists<Object>
186 public Matrix(ArrayList<Object>[] twoDal){
187 this.numberOfRows = twoDal.length;
188 ArrayMaths[] twoD = new ArrayMaths[this.numberOfRows];
189 for(int i=0; i<this.numberOfRows; i++){
190 twoD[i] = new ArrayMaths(twoDal[i]);
191 }
192
193 this.numberOfColumns = twoD[0].length();
194 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
195 for(int i=0; i<numberOfRows; i++){
196 double[] arrayh = (twoD[i].copy()).array();
197 if(arrayh.length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
198 this.matrix[i] = arrayh;
199 }
200 this.permutationIndex = new int[numberOfRows];
201 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
202 }
203
204 // Construct matrix with a copy of an existing numberOfRows 1-D array of Vector<Object>
205 public Matrix(Vector<Object>[] twoDv){
206 this.numberOfRows = twoDv.length;
207 ArrayMaths[] twoD = new ArrayMaths[this.numberOfRows];
208 for(int i=0; i<this.numberOfRows; i++){
209 twoD[i] = new ArrayMaths(twoDv[i]);
210 }
211
212 this.numberOfColumns = twoD[0].length();
213 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
214 for(int i=0; i<numberOfRows; i++){
215 double[] arrayh = (twoD[i].copy()).array();
216 if(arrayh.length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
217 this.matrix[i] = arrayh;
218 }
219 this.permutationIndex = new int[numberOfRows];
220 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
221 }
222
223 // Construct matrix with a copy of an existing numberOfRows x numberOfColumns 2-D array of BigDecimals
224 public Matrix(BigDecimal[][] twoD){
225 this.numberOfRows = twoD.length;
226 this.numberOfColumns = twoD[0].length;
227 for(int i=1; i<numberOfRows; i++){
228 if(twoD[i].length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
229 }
230 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
231 for(int i=0; i<numberOfRows; i++){
232 for(int j=0; j<numberOfColumns; j++){
233 this.matrix[i][j] = twoD[i][j].doubleValue();
234 }
235 }
236 this.permutationIndex = new int[numberOfRows];
237 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
238 }
239
240 // Construct matrix with a copy of an existing numberOfRows x numberOfColumns 2-D array of BigIntegers
241 public Matrix(BigInteger[][] twoD){
242 this.numberOfRows = twoD.length;
243 this.numberOfColumns = twoD[0].length;
244 for(int i=1; i<numberOfRows; i++){
245 if(twoD[i].length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
246 }
247 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
248 for(int i=0; i<numberOfRows; i++){
249 for(int j=0; j<numberOfColumns; j++){
250 this.matrix[i][j] = twoD[i][j].doubleValue();
251 }
252 }
253 this.permutationIndex = new int[numberOfRows];
254 for(int i=0;i<numberOfRows;i++)this.permutationIndex[i]=i;
255 }
256
257
258
259 // Construct matrix with a copy of the 2D matrix and permutation index of an existing Matrix bb.
260 public Matrix(Matrix bb){
261 this.numberOfRows = bb.numberOfRows;
262 this.numberOfColumns = bb.numberOfColumns;
263 this.matrix = new double[this.numberOfRows][this.numberOfColumns];
264 for(int i=0; i<numberOfRows; i++){
265 for(int j=0; j<numberOfColumns; j++){
266 this.matrix[i][j] = bb.matrix[i][j];
267 }
268 }
269 this.permutationIndex = Conv.copy(bb.permutationIndex);
270 this.rowSwapIndex = bb.rowSwapIndex;
271 }
272
273
274 // METHODS
275 // SET VALUES
276 // reset value of tiny used to replace zero in LU decompostions
277 // If not set: 1e-100 used
278 public void resetLUzero(double zeroValue){
279 this.tiny = zeroValue;
280 }
281
282 // Set the matrix with a copy of an existing numberOfRows x numberOfColumns 2-D matrix of variables
283 public void setTwoDarray(double[][] aarray){
284 if(this.numberOfRows != aarray.length)throw new IllegalArgumentException("row length of this Matrix differs from that of the 2D array argument");
285 if(this.numberOfColumns != aarray[0].length)throw new IllegalArgumentException("column length of this Matrix differs from that of the 2D array argument");
286 for(int i=0; i<numberOfRows; i++){
287 if(aarray[i].length!=numberOfColumns)throw new IllegalArgumentException("All rows must have the same length");
288 for(int j=0; j<numberOfColumns; j++){
289 this.matrix[i][j]=aarray[i][j];
290 }
291 }
292 }
293
294 // Set an individual array element
295 // i = row index
296 // j = column index
297 // aa = value of the element
298 public void setElement(int i, int j, double aa){
299 this.matrix[i][j]=aa;
300 }
301
302
303 // Set a sub-matrix starting with row index i, column index j
304 public void setSubMatrix(int i, int j, double[][] subMatrix){
305 int k = subMatrix.length;
306 int l = subMatrix[0].length;
307 if(i+k-1>=this.numberOfRows)throw new IllegalArgumentException("Sub-matrix position is outside the row bounds of this Matrix");
308 if(j+l-1>=this.numberOfColumns)throw new IllegalArgumentException("Sub-matrix position is outside the column bounds of this Matrix");
309
310 int m = 0;
311 int n = 0;
312 for(int p=0; p<k; p++){
313 n = 0;
314 for(int q=0; q<l; q++){
315 this.matrix[i+p][j+q] = subMatrix[m][n];
316 n++;
317 }
318 m++;
319 }
320 }
321
322 // Set a sub-matrix starting with row index i, column index j
323 // and ending with row index k, column index l
324 // See setSubMatrix above - this method has been retained for compatibility purposes
325 public void setSubMatrix(int i, int j, int k, int l, double[][] subMatrix){
326 this.setSubMatrix(i, j, subMatrix);
327 }
328
329 // Set a sub-matrix
330 // row = array of row indices
331 // col = array of column indices
332 public void setSubMatrix(int[] row, int[] col, double[][] subMatrix){
333 int n=row.length;
334 int m=col.length;
335 for(int p=0; p<n; p++){
336 for(int q=0; q<m; q++){
337 this.matrix[row[p]][col[q]] = subMatrix[p][q];
338 }
339 }
340 }
341
342 // Get the value of matrixCheck
343 public boolean getMatrixCheck(){
344 return this.matrixCheck;
345 }
346
347 // SPECIAL MATRICES
348 // Construct an identity matrix
349 public static Matrix identityMatrix(int numberOfRows){
350 Matrix special = new Matrix(numberOfRows, numberOfRows);
351 for(int i=0; i<numberOfRows; i++){
352 special.matrix[i][i]=1.0;
353 }
354 return special;
355 }
356
357 // Construct a square unit matrix
358 public static Matrix unitMatrix(int numberOfRows){
359 Matrix special = new Matrix(numberOfRows, numberOfRows);
360 for(int i=0; i<numberOfRows; i++){
361 for(int j=0; j<numberOfRows; j++){
362 special.matrix[i][j]=1.0;
363 }
364 }
365 return special;
366 }
367
368 // Construct a rectangular unit matrix
369 public static Matrix unitMatrix(int numberOfRows, int numberOfColumns){
370 Matrix special = new Matrix(numberOfRows, numberOfColumns);
371 for(int i=0; i<numberOfRows; i++){
372 for(int j=0; j<numberOfColumns; j++){
373 special.matrix[i][j]=1.0;
374 }
375 }
376 return special;
377 }
378
379 // Construct a square scalar matrix
380 public static Matrix scalarMatrix(int numberOfRows, double diagconst){
381 Matrix special = new Matrix(numberOfRows, numberOfRows);
382 double[][] specialArray = special.getArrayReference();
383 for(int i=0; i<numberOfRows; i++){
384 for(int j=i; j<numberOfRows; j++){
385 if(i==j){
386 specialArray[i][j]= diagconst;
387 }
388 }
389 }
390 return special;
391 }
392
393 // Construct a rectangular scalar matrix
394 public static Matrix scalarMatrix(int numberOfRows, int numberOfColumns, double diagconst){
395 Matrix special = new Matrix(numberOfRows, numberOfColumns);
396 double[][] specialArray = special.getArrayReference();
397 for(int i=0; i<numberOfRows; i++){
398 for(int j=i; j<numberOfColumns; j++){
399 if(i==j){
400 specialArray[i][j]= diagconst;
401 }
402 }
403 }
404 return special;
405 }
406
407 // Construct a square diagonal matrix
408 public static Matrix diagonalMatrix(int numberOfRows, double[] diag){
409 if(diag.length!=numberOfRows)throw new IllegalArgumentException("matrix dimension differs from diagonal array length");
410 Matrix special = new Matrix(numberOfRows, numberOfRows);
411 double[][] specialArray = special.getArrayReference();
412 for(int i=0; i<numberOfRows; i++){
413 specialArray[i][i]=diag[i];
414 }
415 return special;
416 }
417
418 // Construct a rectangular diagonal matrix
419 public static Matrix diagonalMatrix(int numberOfRows, int numberOfColumns, double[] diag){
420 if(diag.length!=numberOfRows)throw new IllegalArgumentException("matrix dimension differs from diagonal array length");
421 Matrix special = new Matrix(numberOfRows, numberOfColumns);
422 double[][] specialArray = special.getArrayReference();
423 for(int i=0; i<numberOfRows; i++){
424 for(int j=i; j<numberOfColumns; j++){
425 if(i==j){
426 specialArray[i][j]= diag[i];
427 }
428 }
429 }
430 return special;
431 }
432
433 // GET VALUES
434 // Return the number of rows
435 public int getNumberOfRows(){
436 return this.numberOfRows;
437 }
438
439 // Return the number of rows
440 public int getNrow(){
441 return this.numberOfRows;
442 }
443
444 // Return the number of columns
445 public int getNumberOfColumns(){
446 return this.numberOfColumns;
447 }
448
449 // Return the number of columns
450 public int getNcol(){
451 return this.numberOfColumns;
452 }
453
454 // Return a reference to the internal 2-D array
455 public double[][] getArrayReference(){
456 return this.matrix;
457 }
458
459 // Return a reference to the internal 2-D array
460 // included for backward compatibility with incorrect earlier documentation
461 public double[][] getArrayPointer(){
462 return this.matrix;
463 }
464
465 // Return a copy of the internal 2-D array
466 public double[][] getArrayCopy(){
467 double[][] c = new double[this.numberOfRows][this.numberOfColumns];
468 for(int i=0; i<numberOfRows; i++){
469 for(int j=0; j<numberOfColumns; j++){
470 c[i][j]=this.matrix[i][j];
471 }
472 }
473 return c;
474 }
475
476 // Return a copy of a row
477 public double[] getRowCopy(int i){
478 if(i>=this.numberOfRows)throw new IllegalArgumentException("Row index, " + i + ", must be less than the number of rows, " + this.numberOfRows);
479 if(i<0)throw new IllegalArgumentException("Row index, " + i + ", must be zero or positive");
480 return Conv.copy(this.matrix[i]);
481 }
482
483 // Return a copy of a column
484 public double[] getColumnCopy(int ii){
485 if(ii>=this.numberOfColumns)throw new IllegalArgumentException("Column index, " + ii + ", must be less than the number of columns, " + this.numberOfColumns);
486 if(ii<0)throw new IllegalArgumentException("column index, " + ii + ", must be zero or positive");
487 double[] col = new double[this.numberOfRows];
488 for(int i=0; i<numberOfRows; i++){
489 col[i]=this.matrix[i][ii];
490 }
491 return col;
492 }
493
494
495 // Return a single element of the internal 2-D array
496 public double getElement(int i, int j){
497 return this.matrix[i][j];
498 }
499
500 // Return a single element of the internal 2-D array
501 // included for backward compatibility with incorrect earlier documentation
502 public double getElementCopy(int i, int j){
503 return this.matrix[i][j];
504 }
505
506 // Return a single element of the internal 2-D array
507 // included for backward compatibility with incorrect earlier documentation
508 public double getElementPointer(int i, int j){
509 return this.matrix[i][j];
510 }
511
512 // Return a sub-matrix starting with row index i, column index j
513 // and ending with row index k, column index l
514 public Matrix getSubMatrix(int i, int j, int k, int l){
515 if(i>k)throw new IllegalArgumentException("row indices inverted");
516 if(j>l)throw new IllegalArgumentException("column indices inverted");
517 if(k>=this.numberOfRows)throw new IllegalArgumentException("Sub-matrix position is outside the row bounds of this Matrix" );
518 if(l>=this.numberOfColumns)throw new IllegalArgumentException("Sub-matrix position is outside the column bounds of this Matrix" + i + " " +l);
519
520 int n=k-i+1, m=l-j+1;
521 Matrix subMatrix = new Matrix(n, m);
522 double[][] sarray = subMatrix.getArrayReference();
523 for(int p=0; p<n; p++){
524 for(int q=0; q<m; q++){
525 sarray[p][q]= this.matrix[i+p][j+q];
526 }
527 }
528 return subMatrix;
529 }
530
531 // Return a sub-matrix
532 // row = array of row indices
533 // col = array of column indices
534 public Matrix getSubMatrix(int[] row, int[] col){
535 int n = row.length;
536 int m = col.length;
537 Matrix subMatrix = new Matrix(n, m);
538 double[][] sarray = subMatrix.getArrayReference();
539 for(int i=0; i<n; i++){
540 for(int j=0; j<m; j++){
541 sarray[i][j]= this.matrix[row[i]][col[j]];
542 }
543 }
544 return subMatrix;
545 }
546
547 // Return a reference to the permutation index array
548 public int[] getIndexReference(){
549 return this.permutationIndex;
550 }
551
552 // Return a reference to the permutation index array
553 // included for backward compatibility with incorrect earlier documentation
554 public int[] getIndexPointer(){
555 return this.permutationIndex;
556 }
557
558 // Return a copy of the permutation index array
559 public int[] getIndexCopy(){
560 int[] indcopy = new int[this.numberOfRows];
561 for(int i=0; i<this.numberOfRows; i++){
562 indcopy[i]=this.permutationIndex[i];
563 }
564 return indcopy;
565 }
566
567 // Return the row swap index
568 public double getSwap(){
569 return this.rowSwapIndex;
570 }
571
572 // COPY
573 // Copy a Matrix [static method]
574 public static Matrix copy(Matrix a){
575 if(a==null){
576 return null;
577 }
578 else{
579 int nr = a.getNumberOfRows();
580 int nc = a.getNumberOfColumns();
581 double[][] aarray = a.getArrayReference();
582 Matrix b = new Matrix(nr,nc);
583 b.numberOfRows = nr;
584 b.numberOfColumns = nc;
585 double[][] barray = b.getArrayReference();
586 for(int i=0; i<nr; i++){
587 for(int j=0; j<nc; j++){
588 barray[i][j]=aarray[i][j];
589 }
590 }
591 for(int i=0; i<nr; i++)b.permutationIndex[i] = a.permutationIndex[i];
592 return b;
593 }
594 }
595
596 // Copy a Matrix [instance method]
597 public Matrix copy(){
598 if(this==null){
599 return null;
600 }
601 else{
602 int nr = this.numberOfRows;
603 int nc = this.numberOfColumns;
604 Matrix b = new Matrix(nr,nc);
605 double[][] barray = b.getArrayReference();
606 b.numberOfRows = nr;
607 b.numberOfColumns = nc;
608 for(int i=0; i<nr; i++){
609 for(int j=0; j<nc; j++){
610 barray[i][j]=this.matrix[i][j];
611 }
612 }
613 for(int i=0; i<nr; i++)b.permutationIndex[i] = this.permutationIndex[i];
614 return b;
615 }
616 }
617
618 // Clone a Matrix
619 public Object clone(){
620 if(this==null){
621 return null;
622 }
623 else{
624 int nr = this.numberOfRows;
625 int nc = this.numberOfColumns;
626 Matrix b = new Matrix(nr,nc);
627 double[][] barray = b.getArrayReference();
628 b.numberOfRows = nr;
629 b.numberOfColumns = nc;
630 for(int i=0; i<nr; i++){
631 for(int j=0; j<nc; j++){
632 barray[i][j]=this.matrix[i][j];
633 }
634 }
635 for(int i=0; i<nr; i++)b.permutationIndex[i] = this.permutationIndex[i];
636 return (Object) b;
637 }
638 }
639
640 // COLUMN MATRICES
641 // Converts a 1-D array of doubles to a column matrix
642 public static Matrix columnMatrix(double[] darray){
643 int nr = darray.length;
644 Matrix pp = new Matrix(nr, 1);
645 for(int i=0; i<nr; i++)pp.matrix[i][0] = darray[i];
646 return pp;
647 }
648
649 // ROW MATRICES
650 // Converts a 1-D array of doubles to a row matrix
651 public static Matrix rowMatrix(double[] darray){
652 int nc = darray.length;
653 Matrix pp = new Matrix(1, nc);
654 for(int i=0; i<nc; i++)pp.matrix[0][i] = darray[i];
655 return pp;
656 }
657
658 // ADDITION
659 // Add this matrix to matrix B. This matrix remains unaltered [instance method]
660 public Matrix plus(Matrix bmat){
661 if((this.numberOfRows!=bmat.numberOfRows)||(this.numberOfColumns!=bmat.numberOfColumns)){
662 throw new IllegalArgumentException("Array dimensions do not agree");
663 }
664 int nr=bmat.numberOfRows;
665 int nc=bmat.numberOfColumns;
666 Matrix cmat = new Matrix(nr,nc);
667 double[][] carray = cmat.getArrayReference();
668 for(int i=0; i<nr; i++){
669 for(int j=0; j<nc; j++){
670 carray[i][j]=this.matrix[i][j] + bmat.matrix[i][j];
671 }
672 }
673 return cmat;
674 }
675
676 // Add this matrix to 2-D array B. This matrix remains unaltered [instance method]
677 public Matrix plus(double[][] bmat){
678 int nr=bmat.length;
679 int nc=bmat[0].length;
680 if((this.numberOfRows!=nr)||(this.numberOfColumns!=nc)){
681 throw new IllegalArgumentException("Array dimensions do not agree");
682 }
683
684 Matrix cmat = new Matrix(nr,nc);
685 double[][] carray = cmat.getArrayReference();
686 for(int i=0; i<nr; i++){
687 for(int j=0; j<nc; j++){
688 carray[i][j]=this.matrix[i][j] + bmat[i][j];
689 }
690 }
691 return cmat;
692 }
693
694
695 // Add matrices A and B [static method]
696 public static Matrix plus(Matrix amat, Matrix bmat){
697 if((amat.numberOfRows!=bmat.numberOfRows)||(amat.numberOfColumns!=bmat.numberOfColumns)){
698 throw new IllegalArgumentException("Array dimensions do not agree");
699 }
700 int nr=amat.numberOfRows;
701 int nc=amat.numberOfColumns;
702 Matrix cmat = new Matrix(nr,nc);
703 double[][] carray = cmat.getArrayReference();
704 for(int i=0; i<nr; i++){
705 for(int j=0; j<nc; j++){
706 carray[i][j]=amat.matrix[i][j] + bmat.matrix[i][j];
707 }
708 }
709 return cmat;
710 }
711
712 // Add matrix B to this matrix [equivalence of +=]
713 public void plusEquals(Matrix bmat){
714 if((this.numberOfRows!=bmat.numberOfRows)||(this.numberOfColumns!=bmat.numberOfColumns)){
715 throw new IllegalArgumentException("Array dimensions do not agree");
716 }
717 int nr=bmat.numberOfRows;
718 int nc=bmat.numberOfColumns;
719
720 for(int i=0; i<nr; i++){
721 for(int j=0; j<nc; j++){
722 this.matrix[i][j] += bmat.matrix[i][j];
723 }
724 }
725 }
726
727 // SUBTRACTION
728 // Subtract matrix B from this matrix. This matrix remains unaltered [instance method]
729 public Matrix minus(Matrix bmat){
730 if((this.numberOfRows!=bmat.numberOfRows)||(this.numberOfColumns!=bmat.numberOfColumns)){
731 throw new IllegalArgumentException("Array dimensions do not agree");
732 }
733 int nr=this.numberOfRows;
734 int nc=this.numberOfColumns;
735 Matrix cmat = new Matrix(nr,nc);
736 double[][] carray = cmat.getArrayReference();
737 for(int i=0; i<nr; i++){
738 for(int j=0; j<nc; j++){
739 carray[i][j]=this.matrix[i][j] - bmat.matrix[i][j];
740 }
741 }
742 return cmat;
743 }
744
745 // Subtract a 2-D array from this matrix. This matrix remains unaltered [instance method]
746 public Matrix minus(double[][] bmat){
747 int nr=bmat.length;
748 int nc=bmat[0].length;
749 if((this.numberOfRows!=nr)||(this.numberOfColumns!=nc)){
750 throw new IllegalArgumentException("Array dimensions do not agree");
751 }
752
753 Matrix cmat = new Matrix(nr,nc);
754 double[][] carray = cmat.getArrayReference();
755 for(int i=0; i<nr; i++){
756 for(int j=0; j<nc; j++){
757 carray[i][j]=this.matrix[i][j] - bmat[i][j];
758 }
759 }
760 return cmat;
761 }
762
763
764 // Subtract matrix B from matrix A [static method]
765 public static Matrix minus(Matrix amat, Matrix bmat){
766 if((amat.numberOfRows!=bmat.numberOfRows)||(amat.numberOfColumns!=bmat.numberOfColumns)){
767 throw new IllegalArgumentException("Array dimensions do not agree");
768 }
769 int nr=amat.numberOfRows;
770 int nc=amat.numberOfColumns;
771 Matrix cmat = new Matrix(nr,nc);
772 double[][] carray = cmat.getArrayReference();
773 for(int i=0; i<nr; i++){
774 for(int j=0; j<nc; j++){
775 carray[i][j]=amat.matrix[i][j] - bmat.matrix[i][j];
776 }
777 }
778 return cmat;
779 }
780
781 // Subtract matrix B from this matrix [equivlance of -=]
782 public void minusEquals(Matrix bmat){
783 if((this.numberOfRows!=bmat.numberOfRows)||(this.numberOfColumns!=bmat.numberOfColumns)){
784 throw new IllegalArgumentException("Array dimensions do not agree");
785 }
786 int nr=bmat.numberOfRows;
787 int nc=bmat.numberOfColumns;
788
789 for(int i=0; i<nr; i++){
790 for(int j=0; j<nc; j++){
791 this.matrix[i][j] -= bmat.matrix[i][j];
792 }
793 }
794 }
795
796 // MULTIPLICATION
797 // Multiply this matrix by a matrix. [instance method]
798 // This matrix remains unaltered.
799 public Matrix times(Matrix bmat){
800 if(this.numberOfColumns!=bmat.numberOfRows)throw new IllegalArgumentException("Nonconformable matrices");
801
802 Matrix cmat = new Matrix(this.numberOfRows, bmat.numberOfColumns);
803 double [][] carray = cmat.getArrayReference();
804 double sum = 0.0D;
805
806 for(int i=0; i<this.numberOfRows; i++){
807 for(int j=0; j<bmat.numberOfColumns; j++){
808 sum=0.0D;
809 for(int k=0; k<this.numberOfColumns; k++){
810 sum += this.matrix[i][k]*bmat.matrix[k][j];
811 }
812 carray[i][j]=sum;
813 }
814 }
815 return cmat;
816 }
817
818 // Multiply this matrix by a 2-D array. [instance method]
819 // This matrix remains unaltered.
820 public Matrix times(double[][] bmat){
821 int nr=bmat.length;
822 int nc=bmat[0].length;
823
824 if(this.numberOfColumns!=nr)throw new IllegalArgumentException("Nonconformable matrices");
825
826 Matrix cmat = new Matrix(this.numberOfRows, nc);
827 double [][] carray = cmat.getArrayReference();
828 double sum = 0.0D;
829
830 for(int i=0; i<this.numberOfRows; i++){
831 for(int j=0; j<nc; j++){
832 sum=0.0D;
833 for(int k=0; k<this.numberOfColumns; k++){
834 sum += this.matrix[i][k]*bmat[k][j];
835 }
836 carray[i][j]=sum;
837 }
838 }
839 return cmat;
840 }
841
842 // Multiply this matrix by a constant [instance method]
843 // This matrix remains unaltered
844 public Matrix times(double constant){
845 Matrix cmat = new Matrix(this.numberOfRows, this.numberOfColumns);
846 double [][] carray = cmat.getArrayReference();
847
848 for(int i=0; i<this.numberOfRows; i++){
849 for(int j=0; j<this.numberOfColumns; j++){
850 carray[i][j] = this.matrix[i][j]*constant;
851 }
852 }
853 return cmat;
854 }
855
856 // Multiply two matrices {static method]
857 public static Matrix times(Matrix amat, Matrix bmat){
858 if(amat.numberOfColumns!=bmat.numberOfRows)throw new IllegalArgumentException("Nonconformable matrices");
859
860 Matrix cmat = new Matrix(amat.numberOfRows, bmat.numberOfColumns);
861 double [][] carray = cmat.getArrayReference();
862 double sum = 0.0D;
863
864 for(int i=0; i<amat.numberOfRows; i++){
865 for(int j=0; j<bmat.numberOfColumns; j++){
866 sum=0.0D;
867 for(int k=0; k<amat.numberOfColumns; k++){
868 sum += (amat.matrix[i][k]*bmat.matrix[k][j]);
869 }
870 carray[i][j]=sum;
871 }
872 }
873 return cmat;
874 }
875
876 // Multiply a Matrix by a 2-D array of doubles [static method]
877 public static Matrix times(Matrix amat, double[][] bmat){
878 if(amat.numberOfColumns!=bmat.length)throw new IllegalArgumentException("Nonconformable matrices");
879
880 Matrix cmat = new Matrix(amat.numberOfRows, bmat[0].length);
881 Matrix dmat = new Matrix(bmat);
882 double [][] carray = cmat.getArrayReference();
883 double sum = 0.0D;
884
885 for(int i=0; i<amat.numberOfRows; i++){
886 for(int j=0; j<dmat.numberOfColumns; j++){
887 sum=0.0D;
888 for(int k=0; k<amat.numberOfColumns; k++){
889 sum += (amat.matrix[i][k]*dmat.matrix[k][j]);
890 }
891 carray[i][j]=sum;
892 }
893 }
894 return cmat;
895 }
896
897 // Multiply a matrix by a constant [static method]
898 public static Matrix times(Matrix amat, double constant){
899 Matrix cmat = new Matrix(amat.numberOfRows, amat.numberOfColumns);
900 double [][] carray = cmat.getArrayReference();
901
902 for(int i=0; i<amat.numberOfRows; i++){
903 for(int j=0; j<amat.numberOfColumns; j++){
904 carray[i][j] = amat.matrix[i][j]*constant;
905 }
906 }
907 return cmat;
908 }
909
910 // Multiply this matrix by a matrix [equivalence of *=]
911 public void timesEquals(Matrix bmat){
912 if(this.numberOfColumns!=bmat.numberOfRows)throw new IllegalArgumentException("Nonconformable matrices");
913
914 Matrix cmat = new Matrix(this.numberOfRows, bmat.numberOfColumns);
915 double [][] carray = cmat.getArrayReference();
916 double sum = 0.0D;
917
918 for(int i=0; i<this.numberOfRows; i++){
919 for(int j=0; j<bmat.numberOfColumns; j++){
920 sum=0.0D;
921 for(int k=0; k<this.numberOfColumns; k++){
922 sum += this.matrix[i][k]*bmat.matrix[k][j];
923 }
924 carray[i][j]=sum;
925 }
926 }
927
928 this.numberOfRows = cmat.numberOfRows;
929 this.numberOfColumns = cmat.numberOfColumns;
930 for(int i=0; i<this.numberOfRows; i++){
931 for(int j=0; j<this.numberOfColumns; j++){
932 this.matrix[i][j] = cmat.matrix[i][j];
933 }
934 }
935 }
936
937 // Multiply this matrix by a constant [equivalence of *=]
938 public void timesEquals(double constant){
939
940 for(int i=0; i<this.numberOfRows; i++){
941 for(int j=0; j<this.numberOfColumns; j++){
942 this.matrix[i][j] *= constant;
943 }
944 }
945 }
946
947 // DIVISION
948 // Divide this Matrix by a Matrix - instance method
949 public Matrix over(Matrix bmat){
950 if((this.numberOfRows!=bmat.numberOfRows)||(this.numberOfColumns!=bmat.numberOfColumns)){
951 throw new IllegalArgumentException("Array dimensions do not agree");
952 }
953 return this.times(bmat.inverse());
954 }
955
956 // Divide a Matrix by a Matrix - static method.
957 public Matrix over(Matrix amat, Matrix bmat){
958 if((amat.numberOfRows!=bmat.numberOfRows)||(amat.numberOfColumns!=bmat.numberOfColumns)){
959 throw new IllegalArgumentException("Array dimensions do not agree");
960 }
961 return amat.times(bmat.inverse());
962 }
963
964
965 // Divide this Matrix by a 2-D array of doubles.
966 public Matrix over(double[][] bmat){
967 int nr=bmat.length;
968 int nc=bmat[0].length;
969 if((this.numberOfRows!=nr)||(this.numberOfColumns!=nc)){
970 throw new IllegalArgumentException("Array dimensions do not agree");
971 }
972
973 Matrix cmat = new Matrix(bmat);
974 return this.times(cmat.inverse());
975 }
976
977 // Divide a Matrix by a 2-D array of doubles - static method.
978 public Matrix over(Matrix amat, double[][] bmat){
979 int nr=bmat.length;
980 int nc=bmat[0].length;
981 if((amat.numberOfRows!=nr)||(amat.numberOfColumns!=nc)){
982 throw new IllegalArgumentException("Array dimensions do not agree");
983 }
984
985 Matrix cmat = new Matrix(bmat);
986 return amat.times(cmat.inverse());
987 }
988
989 // Divide a 2-D array of doubles by a Matrix - static method.
990 public Matrix over(double[][] amat, Matrix bmat){
991 int nr=amat.length;
992 int nc=amat[0].length;
993 if((bmat.numberOfRows!=nr)||(bmat.numberOfColumns!=nc)){
994 throw new IllegalArgumentException("Array dimensions do not agree");
995 }
996
997 Matrix cmat = new Matrix(amat);
998 return cmat.times(bmat.inverse());
999 }
1000
1001 // Divide a 2-D array of doubles by a 2-D array of doubles - static method.
1002 public Matrix over(double[][] amat, double[][] bmat){
1003 int nr=amat.length;
1004 int nc=amat[0].length;
1005 if((bmat.length!=nr)||(bmat[0].length!=nc)){
1006 throw new IllegalArgumentException("Array dimensions do not agree");
1007 }
1008
1009 Matrix cmat = new Matrix(amat);
1010 Matrix dmat = new Matrix(bmat);
1011 return cmat.times(dmat.inverse());
1012 }
1013
1014 // Divide a this matrix by a matrix[equivalence of /=]
1015 public void overEquals(Matrix bmat){
1016 if((this.numberOfRows!=bmat.numberOfRows)||(this.numberOfColumns!=bmat.numberOfColumns)){
1017 throw new IllegalArgumentException("Array dimensions do not agree");
1018 }
1019 Matrix cmat = new Matrix(bmat);
1020 this.timesEquals(cmat.inverse());
1021 }
1022
1023 // Divide this Matrix by a 2D array of doubles [equivalence of /=]
1024 public void overEquals(double[][] bmat){
1025 Matrix pmat = new Matrix(bmat);
1026 this.overEquals(pmat);
1027 }
1028
1029 // INVERSE
1030 // Inverse of a square matrix [instance method]
1031 public Matrix inverse(){
1032 int n = this.numberOfRows;
1033 if(n!=this.numberOfColumns)throw new IllegalArgumentException("Matrix is not square");
1034 Matrix invmat = new Matrix(n, n);
1035
1036 if(n==1){
1037 double[][] hold = this.getArrayCopy();
1038 if(hold[0][0]==0.0)throw new IllegalArgumentException("Matrix is singular");
1039 hold[0][0] = 1.0/hold[0][0];
1040 invmat = new Matrix(hold);
1041 }
1042 else{
1043 if(n==2){
1044 double[][] hold = this.getArrayCopy();
1045 double det = hold[0][0]*hold[1][1] - hold[0][1]*hold[1][0];
1046 if(det==0.0)throw new IllegalArgumentException("Matrix is singular");
1047 double[][] hold2 = new double[2][2];
1048 hold2[0][0] = hold[1][1]/det;
1049 hold2[1][1] = hold[0][0]/det;
1050 hold2[1][0] = -hold[1][0]/det;
1051 hold2[0][1] = -hold[0][1]/det;
1052 invmat = new Matrix(hold2);
1053 }
1054 else{
1055 double[] col = new double[n];
1056 double[] xvec = new double[n];
1057 double[][] invarray = invmat.getArrayReference();
1058 Matrix ludmat;
1059
1060 ludmat = this.luDecomp();
1061 for(int j=0; j<n; j++){
1062 for(int i=0; i<n; i++)col[i]=0.0D;
1063 col[j]=1.0;
1064 xvec=ludmat.luBackSub(col);
1065 for(int i=0; i<n; i++)invarray[i][j]=xvec[i];
1066 }
1067 }
1068 }
1069 return invmat;
1070 }
1071
1072 // Inverse of a square matrix [static method]
1073 public static Matrix inverse(Matrix amat){
1074 int n = amat.numberOfRows;
1075 if(n!=amat.numberOfColumns)throw new IllegalArgumentException("Matrix is not square");
1076 Matrix invmat = new Matrix(n, n);
1077
1078 if(n==1){
1079 double[][] hold = amat.getArrayCopy();
1080 if(hold[0][0]==0.0)throw new IllegalArgumentException("Matrix is singular");
1081 hold[0][0] = 1.0/hold[0][0];
1082 invmat = new Matrix(hold);
1083 }
1084 else{
1085 if(n==2){
1086 double[][] hold = amat.getArrayCopy();
1087 double det = hold[0][0]*hold[1][1] - hold[0][1]*hold[1][0];
1088 if(det==0.0)throw new IllegalArgumentException("Matrix is singular");
1089 double[][] hold2 = new double[2][2];
1090 hold2[0][0] = hold[1][1]/det;
1091 hold2[1][1] = hold[0][0]/det;
1092 hold2[1][0] = -hold[1][0]/det;
1093 hold2[0][1] = -hold[0][1]/det;
1094 invmat = new Matrix(hold2);
1095 }
1096 else{
1097 double[] col = new double[n];
1098 double[] xvec = new double[n];
1099 double[][] invarray = invmat.getArrayReference();
1100 Matrix ludmat;
1101
1102 ludmat = amat.luDecomp();
1103 for(int j=0; j<n; j++){
1104 for(int i=0; i<n; i++)col[i]=0.0D;
1105 col[j]=1.0;
1106 xvec=ludmat.luBackSub(col);
1107 for(int i=0; i<n; i++)invarray[i][j]=xvec[i];
1108 }
1109 }
1110 }
1111 return invmat;
1112 }
1113
1114 // TRANSPOSE
1115 // Transpose of a matrix [instance method]
1116 public Matrix transpose(){
1117 Matrix tmat = new Matrix(this.numberOfColumns, this.numberOfRows);
1118 double[][] tarray = tmat.getArrayReference();
1119 for(int i=0; i<this.numberOfColumns; i++){
1120 for(int j=0; j<this.numberOfRows; j++){
1121 tarray[i][j]=this.matrix[j][i];
1122 }
1123 }
1124 return tmat;
1125 }
1126
1127 // Transpose of a matrix [static method]
1128 public static Matrix transpose(Matrix amat){
1129 Matrix tmat = new Matrix(amat.numberOfColumns, amat.numberOfRows);
1130 double[][] tarray = tmat.getArrayReference();
1131 for(int i=0; i<amat.numberOfColumns; i++){
1132 for(int j=0; j<amat.numberOfRows; j++){
1133 tarray[i][j]=amat.matrix[j][i];
1134 }
1135 }
1136 return tmat;
1137 }
1138
1139 // OPPOSITE
1140 // Opposite of a matrix [instance method]
1141 public Matrix opposite(){
1142 Matrix opp = Matrix.copy(this);
1143 for(int i=0; i<this.numberOfRows; i++){
1144 for(int j=0; j<this.numberOfColumns; j++){
1145 opp.matrix[i][j]=-this.matrix[i][j];
1146 }
1147 }
1148 return opp;
1149 }
1150
1151 // Opposite of a matrix [static method]
1152 public static Matrix opposite(Matrix amat){
1153 Matrix opp = Matrix.copy(amat);
1154 for(int i=0; i<amat.numberOfRows; i++){
1155 for(int j=0; j<amat.numberOfColumns; j++){
1156 opp.matrix[i][j]=-amat.matrix[i][j];
1157 }
1158 }
1159 return opp;
1160 }
1161
1162 // TRACE
1163 // Trace of a matrix [instance method]
1164 public double trace(){
1165 double trac = 0.0D;
1166 for(int i=0; i<Math.min(this.numberOfColumns,this.numberOfColumns); i++){
1167 trac += this.matrix[i][i];
1168 }
1169 return trac;
1170 }
1171
1172 // Trace of a matrix [static method]
1173 public static double trace(Matrix amat){
1174 double trac = 0.0D;
1175 for(int i=0; i<Math.min(amat.numberOfColumns,amat.numberOfColumns); i++){
1176 trac += amat.matrix[i][i];
1177 }
1178 return trac;
1179 }
1180
1181 // DETERMINANT
1182 // Returns the determinant of a square matrix [instance method]
1183 public double determinant(){
1184 int n = this.numberOfRows;
1185 if(n!=this.numberOfColumns)throw new IllegalArgumentException("Matrix is not square");
1186 double det = 0.0D;
1187 if(n==2){
1188 det = this.matrix[0][0]*this.matrix[1][1] - this.matrix[0][1]*this.matrix[1][0];
1189 }
1190 else{
1191 Matrix ludmat = this.luDecomp();
1192 det = ludmat.rowSwapIndex;
1193 for(int j=0; j<n; j++){
1194 det *= ludmat.matrix[j][j];
1195 }
1196 }
1197 return det;
1198 }
1199
1200 // Returns the determinant of a square matrix [static method] - Matrix input
1201 public static double determinant(Matrix amat){
1202 int n = amat.numberOfRows;
1203 if(n!=amat.numberOfColumns)throw new IllegalArgumentException("Matrix is not square");
1204 double det = 0.0D;
1205
1206 if(n==2){
1207 double[][] hold = amat.getArrayCopy();
1208 det = hold[0][0]*hold[1][1] - hold[0][1]*hold[1][0];
1209 }
1210 else{
1211 Matrix ludmat = amat.luDecomp();
1212 det = ludmat.rowSwapIndex;
1213 for(int j=0; j<n; j++){
1214 det *= (ludmat.matrix[j][j]);
1215 }
1216 }
1217 return det;
1218 }
1219
1220 // Returns the determinant of a square matrix [static method] - [][] array input
1221 public static double determinant(double[][]mat){
1222 int n = mat.length;
1223 for(int i=0; i<n; i++)if(n!=mat[i].length)throw new IllegalArgumentException("Matrix is not square");
1224 double det = 0.0D;
1225
1226 if(n==2){
1227 det = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
1228 }
1229 else{
1230 Matrix amat = new Matrix(mat);
1231 Matrix ludmat = amat.luDecomp();
1232 det = ludmat.rowSwapIndex;
1233 for(int j=0; j<n; j++){
1234 det *= (ludmat.matrix[j][j]);
1235 }
1236 }
1237 return det;
1238 }
1239
1240 // Returns the log(determinant) of a square matrix [instance method].
1241 // Useful if determinant() underflows or overflows.
1242 public double logDeterminant(){
1243 int n = this.numberOfRows;
1244 if(n!=this.numberOfColumns)throw new IllegalArgumentException("Matrix is not square");
1245 double det = 0.0D;
1246 Matrix ludmat = this.luDecomp();
1247
1248 det = ludmat.rowSwapIndex;
1249 det=Math.log(det);
1250 for(int j=0; j<n; j++){
1251 det += Math.log(ludmat.matrix[j][j]);
1252 }
1253 return det;
1254 }
1255
1256 // Returns the log(determinant) of a square matrix [static method] - matrix input.
1257 // Useful if determinant() underflows or overflows.
1258 public static double logDeterminant(Matrix amat){
1259 int n = amat.numberOfRows;
1260 if(n!=amat.numberOfColumns)throw new IllegalArgumentException("Matrix is not square");
1261 double det = 0.0D;
1262 Matrix ludmat = amat.luDecomp();
1263
1264 det = ludmat.rowSwapIndex;
1265 det=Math.log(det);
1266 for(int j=0; j<n; j++){
1267 det += Math.log(ludmat.matrix[j][j]);
1268 }
1269 return det;
1270 }
1271
1272 // Returns the log(determinant) of a square matrix [static method] double[][] input.
1273 // Useful if determinant() underflows or overflows.
1274 public static double logDeterminant(double[][] mat){
1275 int n = mat.length;
1276 for(int i=0; i<n; i++)if(n!=mat[i].length)throw new IllegalArgumentException("Matrix is not square");
1277 Matrix amat = new Matrix(mat);
1278 return amat.determinant();
1279 }
1280
1281 // REDUCED ROW ECHELON FORM
1282 public Matrix reducedRowEchelonForm(){
1283
1284 double[][] mat = new double[this.numberOfRows][this.numberOfColumns];
1285 for(int i=0; i<this.numberOfRows; i++){
1286 for(int j=0; j<this.numberOfColumns; j++){
1287 mat[i][j] = this.matrix[i][j];
1288 }
1289 }
1290
1291 int leadingCoeff = 0;
1292 int rowPointer = 0;
1293
1294 boolean testOuter = true;
1295 while(testOuter){
1296 int counter = rowPointer;
1297 boolean testInner = true;
1298 while(testInner && mat[counter][leadingCoeff] == 0) {
1299 counter++;
1300 if(counter == this.numberOfRows){
1301 counter = rowPointer;
1302 leadingCoeff++;
1303 if(leadingCoeff == this.numberOfColumns)testInner=false;
1304 }
1305 }
1306 if(testInner){
1307 double[] temp = mat[rowPointer];
1308 mat[rowPointer] = mat[counter];
1309 mat[counter] = temp;
1310
1311 double pivot = mat[rowPointer][leadingCoeff];
1312 for(int j=0; j<this.numberOfColumns; j++)mat[rowPointer][j] /= pivot;
1313
1314 for(int i=0; i<this.numberOfRows; i++){
1315 if (i!=rowPointer) {
1316 pivot = mat[i][leadingCoeff];
1317 for (int j=0; j<this.numberOfColumns; j++)mat[i][j] -= pivot * mat[rowPointer][j];
1318 }
1319 }
1320 leadingCoeff++;
1321 if(leadingCoeff>=this.numberOfColumns)testOuter = false;
1322 }
1323 rowPointer++;
1324 if(rowPointer>=this.numberOfRows || !testInner)testOuter = false;
1325 }
1326
1327 for(int i=0; i<this.numberOfRows; i++){
1328 for (int j=0; j<this.numberOfColumns; j++){
1329 if(mat[i][j]==-0.0)mat[i][j] = 0.0;
1330 }
1331 }
1332
1333 return new Matrix(mat);
1334 }
1335
1336 // FROBENIUS NORM of a matrix
1337 // Sometimes referred to as the EUCLIDEAN NORM
1338 public double frobeniusNorm(){
1339 double norm=0.0D;
1340 for(int i=0; i<this.numberOfRows; i++){
1341 for(int j=0; j<this.numberOfColumns; j++){
1342 norm=hypot(norm, Math.abs(matrix[i][j]));
1343 }
1344 }
1345 return norm;
1346 }
1347
1348
1349 // ONE NORM of a matrix
1350 public double oneNorm(){
1351 double norm = 0.0D;
1352 double sum = 0.0D;
1353 for(int i=0; i<this.numberOfRows; i++){
1354 sum=0.0D;
1355 for(int j=0; j<this.numberOfColumns; j++){
1356 sum+=Math.abs(this.matrix[i][j]);
1357 }
1358 norm=Math.max(norm,sum);
1359 }
1360 return norm;
1361 }
1362
1363 // INFINITY NORM of a matrix
1364 public double infinityNorm(){
1365 double norm = 0.0D;
1366 double sum = 0.0D;
1367 for(int i=0; i<this.numberOfRows; i++){
1368 sum=0.0D;
1369 for(int j=0; j<this.numberOfColumns; j++){
1370 sum+=Math.abs(this.matrix[i][j]);
1371 }
1372 norm=Math.max(norm,sum);
1373 }
1374 return norm;
1375 }
1376
1377 // SUM OF THE ELEMENTS
1378 // Returns sum of all elements
1379 public double sum(){
1380 double sum = 0.0;
1381 for(int i=0; i<this.numberOfRows; i++){
1382 for(int j=0; j<this.numberOfColumns; j++){
1383 sum += this.matrix[i][j];
1384 }
1385 }
1386 return sum;
1387 }
1388
1389 // Returns sums of the rows
1390 public double[] rowSums(){
1391 double[] sums = new double[this.numberOfRows];
1392 for(int i=0; i<this.numberOfRows; i++){
1393 sums[i] = 0.0;
1394 for(int j=0; j<this.numberOfColumns; j++){
1395 sums[i] += this.matrix[i][j];
1396 }
1397 }
1398 return sums;
1399 }
1400
1401 // Returns sums of the columns
1402 public double[] columnSums(){
1403 double[] sums = new double[this.numberOfColumns];
1404 for(int i=0; i<this.numberOfColumns; i++){
1405 sums[i] = 0.0;
1406 for(int j=0; j<this.numberOfRows; j++){
1407 sums[i] += this.matrix[j][i];
1408 }
1409 }
1410 return sums;
1411 }
1412
1413
1414
1415 // MEAN OF THE ELEMENTS
1416 // Returns mean of all elements
1417 public double mean(){
1418 double mean = 0.0;
1419 for(int i=0; i<this.numberOfRows; i++){
1420 for(int j=0; j<this.numberOfColumns; j++){
1421 mean += this.matrix[i][j];
1422 }
1423 }
1424 mean /= this.numberOfRows*this.numberOfColumns;
1425 return mean;
1426 }
1427
1428 // Returns means of the rows
1429 public double[] rowMeans(){
1430 double[] means = new double[this.numberOfRows];
1431 for(int i=0; i<this.numberOfRows; i++){
1432 means[i] = 0.0;
1433 for(int j=0; j<this.numberOfColumns; j++){
1434 means[i] += this.matrix[i][j];
1435 }
1436 means[i] /= this.numberOfColumns;
1437 }
1438 return means;
1439 }
1440
1441 // Returns means of the columns
1442 public double[] columnMeans(){
1443 double[] means = new double[this.numberOfColumns];
1444 for(int i=0; i<this.numberOfColumns; i++){
1445 means[i] = 0.0;
1446 for(int j=0; j<this.numberOfRows; j++){
1447 means[i] += this.matrix[j][i];
1448 }
1449 means[i] /= this.numberOfRows;
1450 }
1451 return means;
1452 }
1453
1454 // SUBTRACT THE MEAN OF THE ELEMENTS
1455 // Returns a matrix whose elements are the original elements minus the mean of all elements
1456 public Matrix subtractMean(){
1457 Matrix mat = new Matrix(this.numberOfRows, this.numberOfColumns);
1458
1459 double mean = 0.0;
1460 for(int i=0; i<this.numberOfRows; i++){
1461 for(int j=0; j<this.numberOfColumns; j++){
1462 mean += this.matrix[i][j];
1463 }
1464 }
1465 mean /= this.numberOfRows*this.numberOfColumns;
1466 for(int i=0; i<this.numberOfRows; i++){
1467 for(int j=0; j<this.numberOfColumns; j++){
1468 mat.matrix[i][j] = this.matrix[i][j] - mean;
1469 }
1470 }
1471 return mat;
1472 }
1473
1474 // Returns a matrix whose rows are the elements are the original row minus the mean of the elements of that row
1475 public Matrix subtractRowMeans(){
1476 Matrix mat = new Matrix(this.numberOfRows, this.numberOfColumns);
1477
1478 for(int i=0; i<this.numberOfRows; i++){
1479 double mean = 0.0;
1480 for(int j=0; j<this.numberOfColumns; j++){
1481 mean += this.matrix[i][j];
1482 }
1483 mean /= this.numberOfColumns;
1484 for(int j=0; j<this.numberOfColumns; j++){
1485 mat.matrix[i][j] = this.matrix[i][j] - mean;
1486 }
1487 }
1488 return mat;
1489 }
1490
1491 // Returns matrix whose columns are the elements are the original column minus the mean of the elements of that olumnc
1492 public Matrix subtractColumnMeans(){
1493 Matrix mat = new Matrix(this.numberOfRows, this.numberOfColumns);
1494
1495 for(int i=0; i<this.numberOfColumns; i++){
1496 double mean = 0.0;
1497 for(int j=0; j<this.numberOfRows; j++){
1498 mean += this.matrix[j][i];
1499 }
1500 mean /= this.numberOfRows;
1501 for(int j=0; j<this.numberOfRows; j++){
1502 mat.matrix[j][i] = this.matrix[j][i] - mean;
1503 }
1504 }
1505 return mat;
1506 }
1507
1508
1509
1510 // MEDIAN OF THE ELEMENTS
1511 // Returns median of all elements
1512 public double median(){
1513 Stat st = new Stat(this.matrix[0]);
1514
1515 for(int i=1; i<this.numberOfRows; i++){
1516 st.concatenate(this.matrix[i]);
1517 }
1518
1519 return st.median();
1520 }
1521
1522 // Returns medians of the rows
1523 public double[] rowMedians(){
1524 double[] medians = new double[this.numberOfRows];
1525 for(int i=0; i<this.numberOfRows; i++){
1526 Stat st = new Stat(this.matrix[i]);
1527 medians[i] = st.median();
1528 }
1529
1530 return medians;
1531 }
1532
1533 // Returns medians of the columns
1534 public double[] columnMedians(){
1535 double[] medians = new double[this.numberOfRows];
1536 for(int i=0; i<this.numberOfColumns; i++){
1537 double[] hold = new double[this.numberOfRows];
1538 for(int j=0; j<this.numberOfRows; j++){
1539 hold[i] = this.matrix[j][i];
1540 }
1541 Stat st = new Stat(hold);
1542 medians[i] = st.median();
1543 }
1544
1545 return medians;
1546 }
1547
1548 // SET THE DENOMINATOR OF THE VARIANCES AND STANDARD DEVIATIONS TO NUMBER OF ELEMENTS, n
1549 // Default value = n-1
1550 public void setDenominatorToN(){
1551 Stat.setStaticDenominatorToN();
1552 }
1553
1554
1555 // VARIANCE OF THE ELEMENTS
1556 // Returns variance of all elements
1557 public double variance(){
1558 Stat st = new Stat(this.matrix[0]);
1559
1560 for(int i=1; i<this.numberOfRows; i++){
1561 st.concatenate(this.matrix[i]);
1562 }
1563
1564 return st.variance();
1565 }
1566
1567 // Returns variances of the rows
1568 public double[] rowVariances(){
1569 double[] variances = new double[this.numberOfRows];
1570 for(int i=0; i<this.numberOfRows; i++){
1571 Stat st = new Stat(this.matrix[i]);
1572 variances[i] = st.variance();
1573 }
1574
1575 return variances;
1576 }
1577
1578 // Returns variances of the columns
1579 public double[] columnVariances(){
1580 double[] variances = new double[this.numberOfColumns];
1581 for(int i=0; i<this.numberOfColumns; i++){
1582 double[] hold = new double[this.numberOfRows];
1583 for(int j=0; j<this.numberOfRows; j++){
1584 hold[i] = this.matrix[j][i];
1585 }
1586 Stat st = new Stat(hold);
1587 variances[i] = st.variance();
1588 }
1589
1590 return variances;
1591 }
1592
1593
1594
1595 // STANDARD DEVIATION OF THE ELEMENTS
1596 // Returns standard deviation of all elements
1597 public double standardDeviation(){
1598 Stat st = new Stat(this.matrix[0]);
1599
1600 for(int i=1; i<this.numberOfRows; i++){
1601 st.concatenate(this.matrix[i]);
1602 }
1603
1604 return st.standardDeviation();
1605 }
1606
1607 // Returns standard deviations of the rows
1608 public double[] rowStandardDeviations(){
1609 double[] standardDeviations = new double[this.numberOfRows];
1610 for(int i=0; i<this.numberOfRows; i++){
1611 Stat st = new Stat(this.matrix[i]);
1612 standardDeviations[i] = st.standardDeviation();
1613 }
1614
1615 return standardDeviations;
1616 }
1617
1618 // Returns standard deviations of the columns
1619 public double[] columnStandardDeviations(){
1620 double[] standardDeviations = new double[this.numberOfColumns];
1621 for(int i=0; i<this.numberOfColumns; i++){
1622 double[] hold = new double[this.numberOfRows];
1623 for(int j=0; j<this.numberOfRows; j++){
1624 hold[i] = this.matrix[j][i];
1625 }
1626 Stat st = new Stat(hold);
1627 standardDeviations[i] = st.standardDeviation();
1628 }
1629
1630 return standardDeviations;
1631 }
1632
1633
1634 // STANDARD ERROR OF THE ELEMENTS
1635 // Returns standard error of all elements
1636 public double stanadardError(){
1637 Stat st = new Stat(this.matrix[0]);
1638
1639 for(int i=1; i<this.numberOfRows; i++){
1640 st.concatenate(this.matrix[i]);
1641 }
1642
1643 return st.standardError();
1644 }
1645
1646 // Returns standard errors of the rows
1647 public double[] rowStandardErrors(){
1648 double[] standardErrors = new double[this.numberOfRows];
1649 for(int i=0; i<this.numberOfRows; i++){
1650 Stat st = new Stat(this.matrix[i]);
1651 standardErrors[i] = st.standardError();
1652 }
1653
1654 return standardErrors;
1655 }
1656
1657 // Returns standard errors of the columns
1658 public double[] columnStandardErrors(){
1659 double[] standardErrors = new double[this.numberOfRows];
1660 for(int i=0; i<this.numberOfColumns; i++){
1661 double[] hold = new double[this.numberOfRows];
1662 for(int j=0; j<this.numberOfRows; j++){
1663 hold[i] = this.matrix[j][i];
1664 }
1665 Stat st = new Stat(hold);
1666 standardErrors[i] = st.standardError();
1667 }
1668
1669 return standardErrors;
1670 }
1671
1672
1673
1674 // MAXIMUM ELEMENT
1675 // Returns the value, row index and column index of the maximum element
1676 public double[] maximumElement(){
1677 double[] ret = new double[3];
1678 double[] holdD = new double[this.numberOfRows];
1679 ArrayMaths am = null;
1680 int[] holdI = new int [this.numberOfRows];
1681 for(int i=0; i<this.numberOfRows; i++){
1682 am = new ArrayMaths(this.matrix[i]);
1683 holdD[i] = am.maximum();
1684 holdI[i] = am.maximumIndex();
1685 }
1686 am = new ArrayMaths(holdD);
1687 ret[0] = am.maximum();
1688 int maxI = am.maximumIndex();
1689 ret[1] = (double)maxI;
1690 ret[2] = (double)holdI[maxI];
1691
1692 return ret;
1693 }
1694
1695 // Returns maxima of the rows
1696 public double[] rowMaxima(){
1697 double[] maxima = new double[this.numberOfRows];
1698 for(int i=0; i<this.numberOfRows; i++){
1699 Stat st = new Stat(this.matrix[i]);
1700 maxima[i] = st.maximum();
1701 }
1702
1703 return maxima;
1704 }
1705
1706 // Returns maxima of the columns
1707 public double[] columnMaxima(){
1708 double[] maxima = new double[this.numberOfRows];
1709 for(int i=0; i<this.numberOfColumns; i++){
1710 double[] hold = new double[this.numberOfRows];
1711 for(int j=0; j<this.numberOfRows; j++){
1712 hold[i] = this.matrix[j][i];
1713 }
1714 Stat st = new Stat(hold);
1715 maxima[i] = st.maximum();
1716 }
1717
1718 return maxima;
1719 }
1720
1721 // MINIMUM ELEMENT
1722 // Returns the value, row index and column index of the minimum element
1723 public double[] minimumElement(){
1724 double[] ret = new double[3];
1725 double[] holdD = new double[this.numberOfRows];
1726 ArrayMaths am = null;
1727 int[] holdI = new int [this.numberOfRows];
1728 for(int i=0; i<this.numberOfRows; i++){
1729 am = new ArrayMaths(this.matrix[i]);
1730 holdD[i] = am.minimum();
1731 holdI[i] = am.minimumIndex();
1732 }
1733 am = new ArrayMaths(holdD);
1734 ret[0] = am.minimum();
1735 int minI = am.minimumIndex();
1736 ret[1] = (double)minI;
1737 ret[2] = (double)holdI[minI];
1738
1739 return ret;
1740 }
1741
1742 // Returns minima of the rows
1743 public double[] rowMinima(){
1744 double[] minima = new double[this.numberOfRows];
1745 for(int i=0; i<this.numberOfRows; i++){
1746 Stat st = new Stat(this.matrix[i]);
1747 minima[i] = st.minimum();
1748 }
1749
1750 return minima;
1751 }
1752
1753 // Returns minima of the columns
1754 public double[] columnMinima(){
1755 double[] minima = new double[this.numberOfRows];
1756 for(int i=0; i<this.numberOfColumns; i++){
1757 double[] hold = new double[this.numberOfRows];
1758 for(int j=0; j<this.numberOfRows; j++){
1759 hold[i] = this.matrix[j][i];
1760 }
1761 Stat st = new Stat(hold);
1762 minima[i] = st.minimum();
1763 }
1764
1765 return minima;
1766 }
1767
1768 // RANGE
1769 // Returns the range of all the elements
1770 public double range(){
1771 return this.maximumElement()[0] - this.minimumElement()[0];
1772 }
1773
1774 // Returns ranges of the rows
1775 public double[] rowRanges(){
1776 double[] ranges = new double[this.numberOfRows];
1777 for(int i=0; i<this.numberOfRows; i++){
1778 Stat st = new Stat(this.matrix[i]);
1779 ranges[i] = st.maximum() - st.minimum();
1780 }
1781
1782 return ranges;
1783 }
1784
1785 // Returns ranges of the columns
1786 public double[] columnRanges(){
1787 double[] ranges = new double[this.numberOfRows];
1788 for(int i=0; i<this.numberOfColumns; i++){
1789 double[] hold = new double[this.numberOfRows];
1790 for(int j=0; j<this.numberOfRows; j++){
1791 hold[i] = this.matrix[j][i];
1792 }
1793 Stat st = new Stat(hold);
1794 ranges[i] = st.maximum() - st.minimum();
1795 }
1796
1797 return ranges;
1798 }
1799
1800 // PIVOT
1801 // Swaps rows and columns to place absolute maximum element in positiom matrix[0][0]
1802 public int[] pivot(){
1803 double[] max = this.maximumElement();
1804 int maxI = (int)max[1];
1805 int maxJ = (int)max[2];
1806 double[] min = this.minimumElement();
1807 int minI = (int)min[1];
1808 int minJ = (int)min[2];
1809 if(Math.abs(min[0])>Math.abs(max[0])){
1810 maxI = minI;
1811 maxJ = minJ;
1812 }
1813 int[] ret = {maxI, maxJ};
1814
1815 double[] hold1 = this.matrix[0];
1816 this.matrix[0] = this.matrix[maxI];
1817 this.matrix[maxI] = hold1;
1818 double hold2 = 0.0;
1819 for(int i=0; i<this.numberOfRows; i++){
1820 hold2 = this.matrix[i][0];
1821 this.matrix[i][0] = this.matrix[i][maxJ];
1822 this.matrix[i][maxJ] = hold2;
1823 }
1824
1825 return ret;
1826 }
1827
1828 // MATRIX TESTS
1829
1830 // Check if a matrix is square
1831 public boolean isSquare(){
1832 boolean test = false;
1833 if(this.numberOfRows==this.numberOfColumns)test = true;
1834 return test;
1835 }
1836
1837 // Check if a matrix is symmetric
1838 public boolean isSymmetric(){
1839 boolean test = true;
1840 if(this.numberOfRows==this.numberOfColumns){
1841 for(int i=0; i<this.numberOfRows; i++){
1842 for(int j=i+1; j<this.numberOfColumns; j++){
1843 if(this.matrix[i][j]!=this.matrix[j][i])test = false;
1844 }
1845 }
1846 }
1847 else{
1848 test = false;
1849 }
1850 return test;
1851 }
1852
1853 // Check if a matrix is zero
1854 public boolean isZero(){
1855 boolean test = true;
1856 for(int i=0; i<this.numberOfRows; i++){
1857 for(int j=0; j<this.numberOfColumns; j++){
1858 if(this.matrix[i][j]!=0.0D)test = false;
1859 }
1860 }
1861 return test;
1862 }
1863
1864 // Check if a matrix is unit
1865 public boolean isUnit(){
1866 boolean test = true;
1867 for(int i=0; i<this.numberOfRows; i++){
1868 for(int j=0; j<this.numberOfColumns; j++){
1869 if(this.matrix[i][j]!=1.0D)test = false;
1870 }
1871 }
1872 return test;
1873 }
1874
1875 // Check if a matrix is diagonal
1876 public boolean isDiagonal(){
1877 boolean test = true;
1878 for(int i=0; i<this.numberOfRows; i++){
1879 for(int j=0; j<this.numberOfColumns; j++){
1880 if(i!=j && this.matrix[i][j]!=0.0D)test = false;
1881 }
1882 }
1883 return test;
1884 }
1885
1886 // Check if a matrix is upper triagonal
1887 public boolean isUpperTriagonal(){
1888 boolean test = true;
1889 for(int i=0; i<this.numberOfRows; i++){
1890 for(int j=0; j<this.numberOfColumns; j++){
1891 if(j<i && this.matrix[i][j]!=0.0D)test = false;
1892 }
1893 }
1894 return test;
1895 }
1896
1897 // Check if a matrix is lower triagonal
1898 public boolean isLowerTriagonal(){
1899 boolean test = true;
1900 for(int i=0; i<this.numberOfRows; i++){
1901 for(int j=0; j<this.numberOfColumns; j++){
1902 if(i>j && this.matrix[i][j]!=0.0D)test = false;
1903 }
1904 }
1905 return test;
1906 }
1907
1908 // Check if a matrix is tridiagonal
1909 public boolean isTridiagonal(){
1910 boolean test = true;
1911 for(int i=0; i<this.numberOfRows; i++){
1912 for(int j=0; j<this.numberOfColumns; j++){
1913 if(i<(j+1) && this.matrix[i][j]!=0.0D)test = false;
1914 if(j>(i+1) && this.matrix[i][j]!=0.0D)test = false;
1915 }
1916 }
1917 return test;
1918 }
1919
1920 // Check if a matrix is upper Hessenberg
1921 public boolean isUpperHessenberg(){
1922 boolean test = true;
1923 for(int i=0; i<this.numberOfRows; i++){
1924 for(int j=0; j<this.numberOfColumns; j++){
1925 if(j<(i+1) && this.matrix[i][j]!=0.0D)test = false;
1926 }
1927 }
1928 return test;
1929 }
1930
1931 // Check if a matrix is lower Hessenberg
1932 public boolean isLowerHessenberg(){
1933 boolean test = true;
1934 for(int i=0; i<this.numberOfRows; i++){
1935 for(int j=0; j<this.numberOfColumns; j++){
1936 if(i>(j+1) && this.matrix[i][j]!=0.0D)test = false;
1937 }
1938 }
1939 return test;
1940 }
1941
1942 // Check if a matrix is a identity matrix
1943 public boolean isIdentity(){
1944 boolean test = true;
1945 if(this.numberOfRows==this.numberOfColumns){
1946 for(int i=0; i<this.numberOfRows; i++){
1947 if(this.matrix[i][i]!=1.0D)test = false;
1948 for(int j=i+1; j<this.numberOfColumns; j++){
1949 if(this.matrix[i][j]!=0.0D)test = false;
1950 if(this.matrix[j][i]!=0.0D)test = false;
1951 }
1952 }
1953 }
1954 else{
1955 test = false;
1956 }
1957 return test;
1958 }
1959
1960 // Check if a matrix is symmetric within a given tolerance
1961 public boolean isNearlySymmetric(double tolerance){
1962 boolean test = true;
1963 if(this.numberOfRows==this.numberOfColumns){
1964 for(int i=0; i<this.numberOfRows; i++){
1965 for(int j=i+1; j<this.numberOfColumns; j++){
1966 if(Math.abs(this.matrix[i][j]-this.matrix[j][i])>Math.abs(tolerance))test = false;
1967 }
1968 }
1969 }
1970 else{
1971 test = false;
1972 }
1973 return test;
1974 }
1975
1976 // Check if a matrix is zero within a given tolerance
1977 public boolean isNearlyZero(double tolerance){
1978 boolean test = true;
1979 for(int i=0; i<this.numberOfRows; i++){
1980 for(int j=0; j<this.numberOfColumns; j++){
1981 if(Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
1982 }
1983 }
1984 return test;
1985 }
1986
1987 // Check if a matrix is unit within a given tolerance
1988 public boolean isNearlyUnit(double tolerance){
1989 boolean test = true;
1990 for(int i=0; i<this.numberOfRows; i++){
1991 for(int j=0; j<this.numberOfColumns; j++){
1992 if(Math.abs(this.matrix[i][j] - 1.0D)>Math.abs(tolerance))test = false;
1993 }
1994 }
1995 return test;
1996 }
1997
1998
1999 // Check if a matrix is upper triagonal within a given tolerance
2000 public boolean isNearlyUpperTriagonal(double tolerance){
2001 boolean test = true;
2002 for(int i=0; i<this.numberOfRows; i++){
2003 for(int j=0; j<this.numberOfColumns; j++){
2004 if(j<i && Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2005 }
2006 }
2007 return test;
2008 }
2009
2010 // Check if a matrix is lower triagonal within a given tolerance
2011 public boolean isNearlyLowerTriagonal(double tolerance){
2012 boolean test = true;
2013 for(int i=0; i<this.numberOfRows; i++){
2014 for(int j=0; j<this.numberOfColumns; j++){
2015 if(i>j && Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2016 }
2017 }
2018 return test;
2019 }
2020
2021
2022
2023 // Check if a matrix is an identy matrix within a given tolerance
2024 public boolean isNearlyIdenty(double tolerance){
2025 boolean test = true;
2026 if(this.numberOfRows==this.numberOfColumns){
2027 for(int i=0; i<this.numberOfRows; i++){
2028 if(Math.abs(this.matrix[i][i]-1.0D)>Math.abs(tolerance))test = false;
2029 for(int j=i+1; j<this.numberOfColumns; j++){
2030 if(Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2031 if(Math.abs(this.matrix[j][i])>Math.abs(tolerance))test = false;
2032 }
2033 }
2034 }
2035 else{
2036 test = false;
2037 }
2038 return test;
2039 }
2040
2041 // Check if a matrix is tridiagonal within a given tolerance
2042 public boolean isTridiagonal(double tolerance){
2043 boolean test = true;
2044 for(int i=0; i<this.numberOfRows; i++){
2045 for(int j=0; j<this.numberOfColumns; j++){
2046 if(i<(j+1) && Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2047 if(j>(i+1) && Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2048 }
2049 }
2050 return test;
2051 }
2052
2053 // Check if a matrix is tridiagonal within a given tolerance
2054 public boolean isNearlyTridiagonal(double tolerance){
2055 boolean test = true;
2056 for(int i=0; i<this.numberOfRows; i++){
2057 for(int j=0; j<this.numberOfColumns; j++){
2058 if(i<(j+1) && Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2059 if(j>(i+1) && Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2060 }
2061 }
2062 return test;
2063 }
2064
2065 // Check if a matrix is upper Hessenberg within a given tolerance
2066 public boolean isNearlyUpperHessenberg(double tolerance){
2067 boolean test = true;
2068 for(int i=0; i<this.numberOfRows; i++){
2069 for(int j=0; j<this.numberOfColumns; j++){
2070 if(j<(i+1) && Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2071 }
2072 }
2073 return test;
2074 }
2075
2076 // Check if a matrix is lower Hessenberg within a given tolerance
2077 public boolean isNearlyLowerHessenberg(double tolerance){
2078 boolean test = true;
2079 for(int i=0; i<this.numberOfRows; i++){
2080 for(int j=0; j<this.numberOfColumns; j++){
2081 if(i>(j+1) && Math.abs(this.matrix[i][j])>Math.abs(tolerance))test = false;
2082 }
2083 }
2084 return test;
2085 }
2086
2087 // Check if a matrix is singular
2088 public boolean isSingular(){
2089 boolean test = false;
2090 double det = this.determinant();
2091 if(det==0.0)test = true;
2092 return test;
2093 }
2094
2095 // Check if a matrix is singular within a given tolerance
2096 public boolean isNearlySingular(double tolerance){
2097 boolean test = false;
2098 double det = this.determinant();
2099 if(Math.abs(det)<=Math.abs(tolerance))test = true;
2100 return test;
2101 }
2102
2103
2104 // Check for identical rows
2105 // Returns the number of pairs of identical rows followed by the row indices of the identical row pairs
2106 public ArrayList<Integer> identicalRows(){
2107 ArrayList<Integer> ret = new ArrayList<Integer>();
2108 int nIdentical = 0;
2109 for(int i=0; i<this.numberOfRows-1; i++){
2110 for(int j=i+1; j<this.numberOfRows; j++){
2111 int m = 0;
2112 for(int k=0; k<this.numberOfColumns; k++){
2113 if(this.matrix[i][k]==this.matrix[j][k])m++;
2114 }
2115 if(m==this.numberOfColumns){
2116 nIdentical++;
2117 ret.add(new Integer(i));
2118 ret.add(new Integer(j));
2119 }
2120 }
2121 }
2122 ret.add(0,new Integer(nIdentical));
2123 return ret;
2124 }
2125
2126 // Check for identical columnss
2127 // Returns the number of pairs of identical columns followed by the column indices of the identical column pairs
2128 public ArrayList<Integer> identicalColumns(){
2129 ArrayList<Integer> ret = new ArrayList<Integer>();
2130 int nIdentical = 0;
2131 for(int i=0; i<this.numberOfColumns; i++){
2132 for(int j=i+1; j<this.numberOfColumns-1; j++){
2133 int m = 0;
2134 for(int k=0; k<this.numberOfRows; k++){
2135 if(this.matrix[k][i]==this.matrix[k][j])m++;
2136 }
2137 if(m==this.numberOfRows){
2138 nIdentical++;
2139 ret.add(new Integer(i));
2140 ret.add(new Integer(j));
2141 }
2142 }
2143 }
2144 ret.add(0,new Integer(nIdentical));
2145 return ret;
2146 }
2147
2148 // Check for zero rows
2149 // Returns the number of columns of all zeros followed by the column indices
2150 public ArrayList<Integer> zeroRows(){
2151 ArrayList<Integer> ret = new ArrayList<Integer>();
2152 int nZero = 0;
2153 for(int i=0; i<this.numberOfRows; i++){
2154 int m = 0;
2155 for(int k=0; k<this.numberOfColumns; k++){
2156 if(this.matrix[i][k]==0.0)m++;
2157 }
2158 if(m==this.numberOfColumns){
2159 nZero++;
2160 ret.add(new Integer(i));
2161 }
2162 }
2163 ret.add(0,new Integer(nZero));
2164 return ret;
2165 }
2166
2167 // Check for zero columns
2168 // Returns the number of columns of all zeros followed by the column indices
2169 public ArrayList<Integer> zeroColumns(){
2170 ArrayList<Integer> ret = new ArrayList<Integer>();
2171 int nZero = 0;
2172 for(int i=0; i<this.numberOfColumns; i++){
2173 int m = 0;
2174 for(int k=0; k<this.numberOfRows; k++){
2175 if(this.matrix[k][i]==0.0)m++;
2176 }
2177 if(m==this.numberOfRows){
2178 nZero++;
2179 ret.add(new Integer(i));
2180 }
2181 }
2182 ret.add(0,new Integer(nZero));
2183 return ret;
2184 }
2185
2186
2187 // LU DECOMPOSITION OF MATRIX A
2188 // For details of LU decomposition
2189 // See Numerical Recipes, The Art of Scientific Computing
2190 // by W H Press, S A Teukolsky, W T Vetterling & B P Flannery
2191 // Cambridge University Press, http://www.nr.com/
2192 // This method has followed their approach but modified to an object oriented language
2193 // Matrix ludmat is the returned LU decompostion
2194 // int[] index is the vector of row permutations
2195 // rowSwapIndex returns +1.0 for even number of row interchanges
2196 // returns -1.0 for odd number of row interchanges
2197 public Matrix luDecomp(){
2198 if(this.numberOfRows!=this.numberOfColumns)throw new IllegalArgumentException("A matrix is not square");
2199 int n = this.numberOfRows;
2200 int imax = 0;
2201 double dum = 0.0D, temp = 0.0D, big = 0.0D;
2202 double[] vv = new double[n];
2203 double sum = 0.0D;
2204 double dumm = 0.0D;
2205
2206 this.matrixCheck = true;
2207
2208 Matrix ludmat = Matrix.copy(this);
2209 double[][] ludarray = ludmat.getArrayReference();
2210
2211 ludmat.rowSwapIndex=1.0D;
2212 for (int i=0;i<n;i++) {
2213 big=0.0D;
2214 for (int j=0;j<n;j++)if ((temp=Math.abs(ludarray[i][j])) > big) big=temp;
2215 if (big == 0.0D){
2216 if(!this.supressErrorMessage){
2217 System.out.println("Attempted LU Decomposition of a singular matrix in Matrix.luDecomp()");
2218 System.out.println("NaN matrix returned and matrixCheck set to false");
2219 }
2220 this.matrixCheck=false;
2221 for(int k=0;k<n;k++)for(int j=0;j<n;j++)ludarray[k][j]=Double.NaN;
2222 return ludmat;
2223 }
2224 vv[i]=1.0/big;
2225 }
2226 for (int j=0;j<n;j++) {
2227 for (int i=0;i<j;i++) {
2228 sum=ludarray[i][j];
2229 for (int k=0;k<i;k++) sum -= ludarray[i][k]*ludarray[k][j];
2230 ludarray[i][j]=sum;
2231 }
2232 big=0.0D;
2233 for (int i=j;i<n;i++) {
2234 sum=ludarray[i][j];
2235 for (int k=0;k<j;k++)
2236 sum -= ludarray[i][k]*ludarray[k][j];
2237 ludarray[i][j]=sum;
2238 if ((dum=vv[i]*Math.abs(sum)) >= big) {
2239 big=dum;
2240 imax=i;
2241 }
2242 }
2243 if (j != imax) {
2244 for (int k=0;k<n;k++) {
2245 dumm=ludarray[imax][k];
2246 ludarray[imax][k]=ludarray[j][k];
2247 ludarray[j][k]=dumm;
2248 }
2249 ludmat.rowSwapIndex = -ludmat.rowSwapIndex;
2250 vv[imax]=vv[j];
2251 }
2252 ludmat.permutationIndex[j]=imax;
2253
2254 if(ludarray[j][j]==0.0D){
2255 ludarray[j][j]=this.tiny;
2256 }
2257 if(j != n-1) {
2258 dumm=1.0/ludarray[j][j];
2259 for (int i=j+1;i<n;i++){
2260 ludarray[i][j]*=dumm;
2261 }
2262 }
2263 }
2264 return ludmat;
2265 }
2266
2267 // Solves the set of n linear equations A.X=B using not A but its LU decomposition
2268 // bvec is the vector B (input)
2269 // xvec is the vector X (output)
2270 // index is the permutation vector produced by luDecomp()
2271 public double[] luBackSub(double[] bvec){
2272 int ii = 0,ip = 0;
2273 int n=bvec.length;
2274 if(n!=this.numberOfColumns)throw new IllegalArgumentException("vector length is not equal to matrix dimension");
2275 if(this.numberOfColumns!=this.numberOfRows)throw new IllegalArgumentException("matrix is not square");
2276 double sum= 0.0D;
2277 double[] xvec=new double[n];
2278 for(int i=0; i<n; i++){
2279 xvec[i]=bvec[i];
2280 }
2281 for (int i=0;i<n;i++) {
2282 ip=this.permutationIndex[i];
2283 sum=xvec[ip];
2284 xvec[ip]=xvec[i];
2285 if (ii==0){
2286 for (int j=ii;j<=i-1;j++){
2287 sum -= this.matrix[i][j]*xvec[j];
2288 }
2289 }
2290 else{
2291 if(sum==0.0) ii=i;
2292 }
2293 xvec[i]=sum;
2294 }
2295 for(int i=n-1;i>=0;i--) {
2296 sum=xvec[i];
2297 for (int j=i+1;j<n;j++){
2298 sum -= this.matrix[i][j]*xvec[j];
2299 }
2300 xvec[i]= sum/matrix[i][i];
2301 }
2302 return xvec;
2303 }
2304
2305 // Solves the set of n linear equations A.X=B
2306 // bvec is the vector B (input)
2307 // xvec is the vector X (output)
2308 public double[] solveLinearSet(double[] bvec){
2309 double[] xvec = null;
2310 if(this.numberOfRows==this.numberOfColumns){
2311 // square matrix - LU decomposition used
2312 Matrix ludmat = this.luDecomp();
2313 xvec = ludmat.luBackSub(bvec);
2314 }
2315 else{
2316 if(this.numberOfRows>this.numberOfColumns){
2317 // overdetermined equations - least squares used - must be used with care
2318 int n = bvec.length;
2319 if(this.numberOfRows!=n)throw new IllegalArgumentException("Overdetermined equation solution - vector length is not equal to matrix column length");
2320 Matrix avecT = this.transpose();
2321 double[][] avec = avecT.getArrayCopy();
2322 Regression reg = new Regression(avec, bvec);
2323 reg.linearGeneral();
2324 xvec = reg.getCoeff();
2325 }
2326 else{
2327 throw new IllegalArgumentException("This class does not handle underdetermined equations");
2328 }
2329 }
2330 return xvec;
2331 }
2332
2333 //Supress printing of LU decompostion failure message
2334 public void supressErrorMessage(){
2335 this.supressErrorMessage = true;
2336 }
2337
2338
2339 // HESSENBERG MARTIX
2340
2341 // Calculates the Hessenberg equivalant of this matrix
2342 public void hessenbergMatrix(){
2343
2344 this.hessenberg = this.getArrayCopy();
2345 double pivot = 0.0D;
2346 int pivotIndex = 0;
2347 double hold = 0.0D;
2348
2349 for(int i = 1; i<this.numberOfRows-1; i++){
2350 // identify pivot
2351 pivot = 0.0D;
2352 pivotIndex = i;
2353 for(int j=i; j<this.numberOfRows; j++){
2354 if(Math.abs(this.hessenberg[j][i-1])> Math.abs(pivot)){
2355 pivot = this.hessenberg[j][i-1];
2356 pivotIndex = j;
2357 }
2358 }
2359
2360 // row and column interchange
2361 if(pivotIndex != i){
2362 for(int j = i-1; j<this.numberOfRows; j++){
2363 hold = this.hessenberg[pivotIndex][j];
2364 this.hessenberg[pivotIndex][j] = this.hessenberg[i][j];
2365 this.hessenberg[i][j] = hold;
2366 }
2367 for(int j = 0; j<this.numberOfRows; j++){
2368 hold = this.hessenberg[j][pivotIndex];
2369 this.hessenberg[j][pivotIndex] = this.hessenberg[j][i];
2370 this.hessenberg[j][i] = hold;
2371 }
2372
2373 // elimination
2374 if(pivot!=0.0){
2375 for(int j=i+1; j<this.numberOfRows; j++){
2376 hold = this.hessenberg[j][i-1];
2377 if(hold!=0.0){
2378 hold /= pivot;
2379 this.hessenberg[j][i-1] = hold;
2380 for(int k=i; k<this.numberOfRows; k++){
2381 this.hessenberg[j][k] -= hold*this.hessenberg[i][k];
2382 }
2383 for(int k=0; k<this.numberOfRows; k++){
2384 this.hessenberg[k][i] += hold*this.hessenberg[k][j];
2385 }
2386 }
2387 }
2388 }
2389 }
2390 }
2391 for(int i = 2; i<this.numberOfRows; i++){
2392 for(int j = 0; j<i-1; j++){
2393 this.hessenberg[i][j] = 0.0;
2394 }
2395 }
2396 this.hessenbergDone = true;
2397 }
2398
2399 // return the Hessenberg equivalent
2400 public double[][] getHessenbergMatrix(){
2401 if(!hessenbergDone)this.hessenbergMatrix();
2402 return this.hessenberg;
2403 }
2404
2405
2406 // EIGEN VALUES AND EIGEN VECTORS
2407 // For a discussion of eigen systems see
2408 // Numerical Recipes, The Art of Scientific Computing
2409 // by W H Press, S A Teukolsky, W T Vetterling & B P Flannery
2410 // Cambridge University Press, http://www.nr.com/
2411 // These methods follow their approach but modified to an object oriented language
2412
2413 // Return eigen values as calculated
2414 public double[] getEigenValues(){
2415 if(!this.eigenDone)symmetricEigen();
2416 return this.eigenValues;
2417 }
2418
2419 // Return eigen values in descending order
2420 public double[] getSortedEigenValues(){
2421 if(!this.eigenDone)symmetricEigen();
2422 return this.sortedEigenValues;
2423 }
2424
2425 // Return eigen vectors as calculated as columns
2426 // Each vector as a column
2427 public double[][] getEigenVectorsAsColumns(){
2428 if(!this.eigenDone)symmetricEigen();
2429 return this.eigenVector;
2430 }
2431 // Return eigen vectors as calculated as columns
2432 // Each vector as a column
2433 public double[][] getEigenVector(){
2434 if(!this.eigenDone)symmetricEigen();
2435 return this.eigenVector;
2436 }
2437
2438 // Return eigen vectors as calculated as rows
2439 // Each vector as a row
2440 public double[][] getEigenVectorsAsRows(){
2441 if(!this.eigenDone)symmetricEigen();
2442 double[][] ret = new double[this.numberOfRows][this.numberOfRows];
2443 for(int i=0; i<this.numberOfRows;i++){
2444 for(int j=0; j<this.numberOfRows;j++){
2445 ret[i][j] = this.eigenVector[j][i];
2446 }
2447 }
2448 return ret;
2449 }
2450
2451 // Return eigen vectors reordered to match a descending order of eigen values
2452 // Each vector as a column
2453 public double[][] getSortedEigenVectorsAsColumns(){
2454 if(!this.eigenDone)symmetricEigen();
2455 return this.sortedEigenVector;
2456 }
2457
2458 // Return eigen vectors reordered to match a descending order of eigen values
2459 // Each vector as a column
2460 public double[][] getSortedEigenVector(){
2461 if(!this.eigenDone)symmetricEigen();
2462 return this.sortedEigenVector;
2463 }
2464
2465 // Return eigen vectors reordered to match a descending order of eigen values
2466 // Each vector as a row
2467 public double[][] getSortedEigenVectorsAsRows(){
2468 if(!this.eigenDone)symmetricEigen();
2469 double[][] ret = new double[this.numberOfRows][this.numberOfRows];
2470 for(int i=0; i<this.numberOfRows;i++){
2471 for(int j=0; j<this.numberOfRows;j++){
2472 ret[i][j] = this.sortedEigenVector[j][i];
2473 }
2474 }
2475 return ret;
2476 }
2477
2478 // Return the number of rotations used in the Jacobi procedure
2479 public int getNumberOfJacobiRotations(){
2480 return this.numberOfRotations;
2481 }
2482
2483 // Returns the eigen values and eigen vectors of a symmetric matrix
2484 // Follows the approach of Numerical methods but adapted to object oriented programming (see above)
2485 private void symmetricEigen(){
2486
2487 if(!this.isSymmetric())throw new IllegalArgumentException("matrix is not symmetric");
2488 double[][] amat = this.getArrayCopy();
2489 this.eigenVector = new double[this.numberOfRows][this.numberOfRows];
2490 this.eigenValues = new double[this.numberOfRows];
2491 double threshold = 0.0D;
2492 double cot2rotationAngle = 0.0D;
2493 double tanHalfRotationAngle = 0.0D;
2494 double offDiagonalSum = 0.0D;
2495 double scaledOffDiagonal = 0.0D;
2496 double sElement = 0.0D;
2497 double cElement = 0.0D;
2498 double sOverC = 0.0D;
2499 double vectorDifference = 0.0D;
2500 double[] holdingVector1 = new double[this.numberOfRows];
2501 double[] holdingVector2 = new double[this.numberOfRows];
2502
2503 for(int p=0;p<this.numberOfRows;p++){
2504 for(int q=0;q<this.numberOfRows;q++) this.eigenVector[p][q] = 0.0;
2505 this.eigenVector[p][p] = 1.0;
2506 }
2507 for(int p=0;p<this.numberOfRows;p++){
2508 holdingVector1[p] = amat[p][p];
2509 this.eigenValues[p] = amat[p][p];
2510 holdingVector2[p] = 0.0;
2511 }
2512 this.numberOfRotations = 0;
2513 for(int i=1;i<=this.maximumJacobiIterations;i++){
2514 offDiagonalSum = 0.0;
2515 for(int p=0;p<this.numberOfRows-1;p++){
2516 for(int q=p+1;q<this.numberOfRows;q++){
2517 offDiagonalSum += Math.abs(amat[p][q]);
2518 }
2519 }
2520 if(offDiagonalSum==0.0){
2521 this.eigenDone = true;
2522 this.eigenSort();
2523 return;
2524 }
2525 if (i < 4){
2526 threshold = 0.2*offDiagonalSum/(this.numberOfRows*this.numberOfRows);
2527 }
2528 else{
2529 threshold = 0.0;
2530 }
2531 for(int p=0;p<this.numberOfRows-1;p++){
2532 for(int q=p+1;q<this.numberOfRows;q++){
2533 scaledOffDiagonal = 100.0*Math.abs(amat[p][q]);
2534 if (i > 4 && (Math.abs(this.eigenValues[p]) + scaledOffDiagonal) == Math.abs(this.eigenValues[p]) && (Math.abs(this.eigenValues[q]) + scaledOffDiagonal) == Math.abs(this.eigenValues[q])){
2535 amat[p][q] = 0.0;
2536 }
2537 else if(Math.abs(amat[p][q]) > threshold){
2538 vectorDifference = this.eigenValues[q] - this.eigenValues[p];
2539 if ((Math.abs(vectorDifference) + scaledOffDiagonal) == Math.abs(vectorDifference))
2540 sOverC = amat[p][q]/vectorDifference;
2541 else{
2542 cot2rotationAngle = 0.5*vectorDifference/amat[p][q];
2543 sOverC = 1.0/(Math.abs(cot2rotationAngle) + Math.sqrt(1.0 + cot2rotationAngle*cot2rotationAngle));
2544 if (cot2rotationAngle < 0.0) sOverC = -sOverC;
2545 }
2546 cElement = 1.0/Math.sqrt(1.0 + sOverC*sOverC);
2547 sElement = sOverC*cElement;
2548 tanHalfRotationAngle = sElement/(1.0 + cElement);
2549 vectorDifference = sOverC*amat[p][q];
2550 holdingVector2[p] -= vectorDifference;
2551 holdingVector2[q] += vectorDifference;
2552 this.eigenValues[p] -= vectorDifference;
2553 this.eigenValues[q] += vectorDifference;
2554 amat[p][q] = 0.0;
2555 for(int j=0;j<=p-1;j++) rotation(amat, tanHalfRotationAngle, sElement, j, p, j, q);
2556 for(int j=p+1;j<=q-1;j++) rotation(amat, tanHalfRotationAngle, sElement, p, j, j, q);
2557 for(int j=q+1;j<this.numberOfRows;j++) rotation(amat, tanHalfRotationAngle, sElement,p, j, q, j);
2558 for(int j=0;j<this.numberOfRows;j++) rotation(this.eigenVector, tanHalfRotationAngle, sElement, j, p, j, q);
2559 ++this.numberOfRotations;
2560 }
2561 }
2562 }
2563 for(int p=0;p<this.numberOfRows;p++){
2564 holdingVector1[p] += holdingVector2[p];
2565 this.eigenValues[p] = holdingVector1[p];
2566 holdingVector2[p] = 0.0;
2567 }
2568 }
2569 System.out.println("Maximum iterations, " + this.maximumJacobiIterations + ", reached - values at this point returned");
2570 this.eigenDone = true;
2571 this.eigenSort();
2572 }
2573
2574 // matrix rotaion required by symmetricEigen
2575 private void rotation(double[][] a, double tau, double sElement, int i, int j, int k, int l){
2576 double aHold1 = a[i][j];
2577 double aHold2 = a[k][l];
2578 a[i][j] = aHold1 - sElement*(aHold2 + aHold1*tau);
2579 a[k][l] = aHold2 + sElement*(aHold1 - aHold2*tau);
2580 }
2581
2582 // Sorts eigen values into descending order and rearranges eigen vecors to match
2583 // follows Numerical Recipes (see above)
2584 private void eigenSort(){
2585 int k = 0;
2586 double holdingElement;
2587 this.sortedEigenValues = Conv.copy(this.eigenValues);
2588 this.sortedEigenVector = Conv.copy(this.eigenVector);
2589 this.eigenIndices = new int[this.numberOfRows];
2590
2591 for(int i=0; i<this.numberOfRows-1; i++){
2592 holdingElement = this.sortedEigenValues[i];
2593 k = i;
2594 for(int j=i+1; j<this.numberOfRows; j++){
2595 if (this.sortedEigenValues[j] >= holdingElement){
2596 holdingElement = this.sortedEigenValues[j];
2597 k = j;
2598 }
2599 }
2600 if (k != i){
2601 this.sortedEigenValues[k] = this.sortedEigenValues[i];
2602 this.sortedEigenValues[i] = holdingElement;
2603
2604 for(int j=0; j<this.numberOfRows; j++){
2605 holdingElement = this.sortedEigenVector[j][i];
2606 this.sortedEigenVector[j][i] = this.sortedEigenVector[j][k];
2607 this.sortedEigenVector[j][k] = holdingElement;
2608 }
2609 }
2610 }
2611 this.eigenIndices = new int[this.numberOfRows];
2612 for(int i=0; i<this.numberOfRows; i++){
2613 boolean test = true;
2614 int j = 0;
2615 while(test){
2616 if(this.sortedEigenValues[i]==this.eigenValues[j]){
2617 this.eigenIndices[i] = j;
2618 test = false;
2619 }
2620 else{
2621 j++;
2622 }
2623 }
2624 }
2625 }
2626
2627 // Return indices of the eigen values before sorting into descending order
2628 public int[] eigenValueIndices(){
2629 if(!this.eigenDone)symmetricEigen();
2630 return this.eigenIndices;
2631 }
2632
2633
2634 // Method not in java.lang.maths required in this Class
2635 // See Fmath.class for public versions of this method
2636 private static double hypot(double aa, double bb){
2637 double cc = 0.0D, ratio = 0.0D;
2638 double amod=Math.abs(aa);
2639 double bmod=Math.abs(bb);
2640
2641 if(amod==0.0D){
2642 cc=bmod;
2643 }
2644 else{
2645 if(bmod==0.0D){
2646 cc=amod;
2647 }
2648 else{
2649 if(amod<=bmod){
2650 ratio=amod/bmod;
2651 cc=bmod*Math.sqrt(1.0D+ratio*ratio);
2652 }
2653 else{
2654 ratio=bmod/amod;
2655 cc=amod*Math.sqrt(1.0D+ratio*ratio);
2656 }
2657 }
2658 }
2659 return cc;
2660 }
2661
2662}
2663
2664
2665
2666
Note: See TracBrowser for help on using the repository browser.