1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.linear;
19
20 import java.io.Serializable;
21 import org.apache.commons.math.util.MathUtils;
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52 public class RealMatrixImpl implements RealMatrix, Serializable {
53
54
55 private static final long serialVersionUID = 4237564493130426188L;
56
57
58 private double data[][] = null;
59
60
61
62
63 private double lu[][] = null;
64
65
66 private int[] permutation = null;
67
68
69 private int parity = 1;
70
71
72 protected static double TOO_SMALL = 10E-12;
73
74
75
76
77 public RealMatrixImpl() {
78 }
79
80
81
82
83
84
85
86
87
88 public RealMatrixImpl(int rowDimension, int columnDimension) {
89 if (rowDimension <= 0 || columnDimension <= 0) {
90 throw new IllegalArgumentException(
91 "row and column dimensions must be postive");
92 }
93 data = new double[rowDimension][columnDimension];
94 lu = null;
95 }
96
97
98
99
100
101
102
103
104
105
106
107
108 public RealMatrixImpl(double[][] d) {
109 this.copyIn(d);
110 lu = null;
111 }
112
113
114
115
116
117
118
119
120
121
122 public RealMatrixImpl(double[] v) {
123 int nRows = v.length;
124 data = new double[nRows][1];
125 for (int row = 0; row < nRows; row++) {
126 data[row][0] = v[row];
127 }
128 }
129
130
131
132
133
134
135 public RealMatrix copy() {
136 return new RealMatrixImpl(this.copyOut());
137 }
138
139
140
141
142
143
144
145
146 public RealMatrix add(RealMatrix m) throws IllegalArgumentException {
147 if (this.getColumnDimension() != m.getColumnDimension() ||
148 this.getRowDimension() != m.getRowDimension()) {
149 throw new IllegalArgumentException("matrix dimension mismatch");
150 }
151 int rowCount = this.getRowDimension();
152 int columnCount = this.getColumnDimension();
153 double[][] outData = new double[rowCount][columnCount];
154 for (int row = 0; row < rowCount; row++) {
155 for (int col = 0; col < columnCount; col++) {
156 outData[row][col] = data[row][col] + m.getEntry(row, col);
157 }
158 }
159 return new RealMatrixImpl(outData);
160 }
161
162
163
164
165
166
167
168
169 public RealMatrix subtract(RealMatrix m) throws IllegalArgumentException {
170 if (this.getColumnDimension() != m.getColumnDimension() ||
171 this.getRowDimension() != m.getRowDimension()) {
172 throw new IllegalArgumentException("matrix dimension mismatch");
173 }
174 int rowCount = this.getRowDimension();
175 int columnCount = this.getColumnDimension();
176 double[][] outData = new double[rowCount][columnCount];
177 for (int row = 0; row < rowCount; row++) {
178 for (int col = 0; col < columnCount; col++) {
179 outData[row][col] = data[row][col] - m.getEntry(row, col);
180 }
181 }
182 return new RealMatrixImpl(outData);
183 }
184
185
186
187
188
189
190
191 public RealMatrix scalarAdd(double d) {
192 int rowCount = this.getRowDimension();
193 int columnCount = this.getColumnDimension();
194 double[][] outData = new double[rowCount][columnCount];
195 for (int row = 0; row < rowCount; row++) {
196 for (int col = 0; col < columnCount; col++) {
197 outData[row][col] = data[row][col] + d;
198 }
199 }
200 return new RealMatrixImpl(outData);
201 }
202
203
204
205
206
207
208 public RealMatrix scalarMultiply(double d) {
209 int rowCount = this.getRowDimension();
210 int columnCount = this.getColumnDimension();
211 double[][] outData = new double[rowCount][columnCount];
212 for (int row = 0; row < rowCount; row++) {
213 for (int col = 0; col < columnCount; col++) {
214 outData[row][col] = data[row][col] * d;
215 }
216 }
217 return new RealMatrixImpl(outData);
218 }
219
220
221
222
223
224
225
226
227 public RealMatrix multiply(RealMatrix m) throws IllegalArgumentException {
228 if (this.getColumnDimension() != m.getRowDimension()) {
229 throw new IllegalArgumentException("Matrices are not multiplication compatible.");
230 }
231 int nRows = this.getRowDimension();
232 int nCols = m.getColumnDimension();
233 int nSum = this.getColumnDimension();
234 double[][] outData = new double[nRows][nCols];
235 double sum = 0;
236 for (int row = 0; row < nRows; row++) {
237 for (int col = 0; col < nCols; col++) {
238 sum = 0;
239 for (int i = 0; i < nSum; i++) {
240 sum += data[row][i] * m.getEntry(i, col);
241 }
242 outData[row][col] = sum;
243 }
244 }
245 return new RealMatrixImpl(outData);
246 }
247
248
249
250
251
252
253
254
255 public RealMatrix preMultiply(RealMatrix m) throws IllegalArgumentException {
256 return m.multiply(this);
257 }
258
259
260
261
262
263
264
265
266 public double[][] getData() {
267 return copyOut();
268 }
269
270
271
272
273
274
275
276
277 public double[][] getDataRef() {
278 return data;
279 }
280
281
282
283
284
285 public double getNorm() {
286 double maxColSum = 0;
287 for (int col = 0; col < this.getColumnDimension(); col++) {
288 double sum = 0;
289 for (int row = 0; row < this.getRowDimension(); row++) {
290 sum += Math.abs(data[row][col]);
291 }
292 maxColSum = Math.max(maxColSum, sum);
293 }
294 return maxColSum;
295 }
296
297
298
299
300
301
302
303
304
305
306
307
308
309 public RealMatrix getSubMatrix(int startRow, int endRow, int startColumn,
310 int endColumn) throws MatrixIndexException {
311 if (startRow < 0 || startRow > endRow || endRow > data.length ||
312 startColumn < 0 || startColumn > endColumn ||
313 endColumn > data[0].length ) {
314 throw new MatrixIndexException(
315 "invalid row or column index selection");
316 }
317 RealMatrixImpl subMatrix = new RealMatrixImpl(endRow - startRow+1,
318 endColumn - startColumn+1);
319 double[][] subMatrixData = subMatrix.getDataRef();
320 for (int i = startRow; i <= endRow; i++) {
321 for (int j = startColumn; j <= endColumn; j++) {
322 subMatrixData[i - startRow][j - startColumn] = data[i][j];
323 }
324 }
325 return subMatrix;
326 }
327
328
329
330
331
332
333
334
335
336
337
338
339 public RealMatrix getSubMatrix(int[] selectedRows, int[] selectedColumns)
340 throws MatrixIndexException {
341 if (selectedRows.length * selectedColumns.length == 0) {
342 throw new MatrixIndexException(
343 "selected row and column index arrays must be non-empty");
344 }
345 RealMatrixImpl subMatrix = new RealMatrixImpl(selectedRows.length,
346 selectedColumns.length);
347 double[][] subMatrixData = subMatrix.getDataRef();
348 try {
349 for (int i = 0; i < selectedRows.length; i++) {
350 for (int j = 0; j < selectedColumns.length; j++) {
351 subMatrixData[i][j] = data[selectedRows[i]][selectedColumns[j]];
352 }
353 }
354 }
355 catch (ArrayIndexOutOfBoundsException e) {
356 throw new MatrixIndexException("matrix dimension mismatch");
357 }
358 return subMatrix;
359 }
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388 public void setSubMatrix(double[][] subMatrix, int row, int column)
389 throws MatrixIndexException {
390 if ((row < 0) || (column < 0)){
391 throw new MatrixIndexException
392 ("invalid row or column index selection");
393 }
394 int nRows = subMatrix.length;
395 if (nRows == 0) {
396 throw new IllegalArgumentException(
397 "Matrix must have at least one row.");
398 }
399 int nCols = subMatrix[0].length;
400 if (nCols == 0) {
401 throw new IllegalArgumentException(
402 "Matrix must have at least one column.");
403 }
404 for (int r = 1; r < nRows; r++) {
405 if (subMatrix[r].length != nCols) {
406 throw new IllegalArgumentException(
407 "All input rows must have the same length.");
408 }
409 }
410 if (data == null) {
411 if ((row > 0)||(column > 0)) throw new MatrixIndexException
412 ("matrix must be initialized to perfom this method");
413 data = new double[nRows][nCols];
414 System.arraycopy(subMatrix, 0, data, 0, subMatrix.length);
415 }
416 if (((nRows + row) > this.getRowDimension()) ||
417 (nCols + column > this.getColumnDimension()))
418 throw new MatrixIndexException(
419 "invalid row or column index selection");
420 for (int i = 0; i < nRows; i++) {
421 System.arraycopy(subMatrix[i], 0, data[row + i], column, nCols);
422 }
423 lu = null;
424 }
425
426
427
428
429
430
431
432
433
434 public RealMatrix getRowMatrix(int row) throws MatrixIndexException {
435 if ( !isValidCoordinate( row, 0)) {
436 throw new MatrixIndexException("illegal row argument");
437 }
438 int ncols = this.getColumnDimension();
439 double[][] out = new double[1][ncols];
440 System.arraycopy(data[row], 0, out[0], 0, ncols);
441 return new RealMatrixImpl(out);
442 }
443
444
445
446
447
448
449
450
451
452 public RealMatrix getColumnMatrix(int column) throws MatrixIndexException {
453 if ( !isValidCoordinate( 0, column)) {
454 throw new MatrixIndexException("illegal column argument");
455 }
456 int nRows = this.getRowDimension();
457 double[][] out = new double[nRows][1];
458 for (int row = 0; row < nRows; row++) {
459 out[row][0] = data[row][column];
460 }
461 return new RealMatrixImpl(out);
462 }
463
464
465
466
467
468
469
470
471
472
473
474 public double[] getRow(int row) throws MatrixIndexException {
475 if ( !isValidCoordinate( row, 0 ) ) {
476 throw new MatrixIndexException("illegal row argument");
477 }
478 int ncols = this.getColumnDimension();
479 double[] out = new double[ncols];
480 System.arraycopy(data[row], 0, out, 0, ncols);
481 return out;
482 }
483
484
485
486
487
488
489
490
491
492
493
494 public double[] getColumn(int col) throws MatrixIndexException {
495 if ( !isValidCoordinate(0, col) ) {
496 throw new MatrixIndexException("illegal column argument");
497 }
498 int nRows = this.getRowDimension();
499 double[] out = new double[nRows];
500 for (int row = 0; row < nRows; row++) {
501 out[row] = data[row][col];
502 }
503 return out;
504 }
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521 public double getEntry(int row, int column)
522 throws MatrixIndexException {
523 if (!isValidCoordinate(row,column)) {
524 throw new MatrixIndexException("matrix entry does not exist");
525 }
526 return data[row][column];
527 }
528
529
530
531
532
533
534 public RealMatrix transpose() {
535 int nRows = this.getRowDimension();
536 int nCols = this.getColumnDimension();
537 RealMatrixImpl out = new RealMatrixImpl(nCols, nRows);
538 double[][] outData = out.getDataRef();
539 for (int row = 0; row < nRows; row++) {
540 for (int col = 0; col < nCols; col++) {
541 outData[col][row] = data[row][col];
542 }
543 }
544 return out;
545 }
546
547
548
549
550
551
552
553 public RealMatrix inverse() throws InvalidMatrixException {
554 return solve(MatrixUtils.createRealIdentityMatrix
555 (this.getRowDimension()));
556 }
557
558
559
560
561
562 public double getDeterminant() throws InvalidMatrixException {
563 if (!isSquare()) {
564 throw new InvalidMatrixException("matrix is not square");
565 }
566 if (isSingular()) {
567 return 0d;
568 } else {
569 double det = parity;
570 for (int i = 0; i < this.getRowDimension(); i++) {
571 det *= lu[i][i];
572 }
573 return det;
574 }
575 }
576
577
578
579
580 public boolean isSquare() {
581 return (this.getColumnDimension() == this.getRowDimension());
582 }
583
584
585
586
587 public boolean isSingular() {
588 if (lu == null) {
589 try {
590 luDecompose();
591 return false;
592 } catch (InvalidMatrixException ex) {
593 return true;
594 }
595 } else {
596 return false;
597 }
598 }
599
600
601
602
603 public int getRowDimension() {
604 return data.length;
605 }
606
607
608
609
610 public int getColumnDimension() {
611 return data[0].length;
612 }
613
614
615
616
617
618 public double getTrace() throws IllegalArgumentException {
619 if (!isSquare()) {
620 throw new IllegalArgumentException("matrix is not square");
621 }
622 double trace = data[0][0];
623 for (int i = 1; i < this.getRowDimension(); i++) {
624 trace += data[i][i];
625 }
626 return trace;
627 }
628
629
630
631
632
633
634 public double[] operate(double[] v) throws IllegalArgumentException {
635 if (v.length != this.getColumnDimension()) {
636 throw new IllegalArgumentException("vector has wrong length");
637 }
638 int nRows = this.getRowDimension();
639 int nCols = this.getColumnDimension();
640 double[] out = new double[v.length];
641 for (int row = 0; row < nRows; row++) {
642 double sum = 0;
643 for (int i = 0; i < nCols; i++) {
644 sum += data[row][i] * v[i];
645 }
646 out[row] = sum;
647 }
648 return out;
649 }
650
651
652
653
654
655
656 public double[] preMultiply(double[] v) throws IllegalArgumentException {
657 int nRows = this.getRowDimension();
658 if (v.length != nRows) {
659 throw new IllegalArgumentException("vector has wrong length");
660 }
661 int nCols = this.getColumnDimension();
662 double[] out = new double[nCols];
663 for (int col = 0; col < nCols; col++) {
664 double sum = 0;
665 for (int i = 0; i < nRows; i++) {
666 sum += data[i][col] * v[i];
667 }
668 out[col] = sum;
669 }
670 return out;
671 }
672
673
674
675
676
677
678
679
680
681
682
683
684 public double[] solve(double[] b) throws IllegalArgumentException, InvalidMatrixException {
685 int nRows = this.getRowDimension();
686 if (b.length != nRows) {
687 throw new IllegalArgumentException("constant vector has wrong length");
688 }
689 RealMatrix bMatrix = new RealMatrixImpl(b);
690 double[][] solution = ((RealMatrixImpl) (solve(bMatrix))).getDataRef();
691 double[] out = new double[nRows];
692 for (int row = 0; row < nRows; row++) {
693 out[row] = solution[row][0];
694 }
695 return out;
696 }
697
698
699
700
701
702
703
704
705
706
707
708
709 public RealMatrix solve(RealMatrix b) throws IllegalArgumentException, InvalidMatrixException {
710 if (b.getRowDimension() != this.getRowDimension()) {
711 throw new IllegalArgumentException("Incorrect row dimension");
712 }
713 if (!this.isSquare()) {
714 throw new InvalidMatrixException("coefficient matrix is not square");
715 }
716 if (this.isSingular()) {
717 throw new InvalidMatrixException("Matrix is singular.");
718 }
719
720 int nCol = this.getColumnDimension();
721 int nColB = b.getColumnDimension();
722 int nRowB = b.getRowDimension();
723
724
725 double[][] bp = new double[nRowB][nColB];
726 for (int row = 0; row < nRowB; row++) {
727 for (int col = 0; col < nColB; col++) {
728 bp[row][col] = b.getEntry(permutation[row], col);
729 }
730 }
731
732
733 for (int col = 0; col < nCol; col++) {
734 for (int i = col + 1; i < nCol; i++) {
735 for (int j = 0; j < nColB; j++) {
736 bp[i][j] -= bp[col][j] * lu[i][col];
737 }
738 }
739 }
740
741
742 for (int col = nCol - 1; col >= 0; col--) {
743 for (int j = 0; j < nColB; j++) {
744 bp[col][j] /= lu[col][col];
745 }
746 for (int i = 0; i < col; i++) {
747 for (int j = 0; j < nColB; j++) {
748 bp[i][j] -= bp[col][j] * lu[i][col];
749 }
750 }
751 }
752
753 RealMatrixImpl outMat = new RealMatrixImpl(bp);
754 return outMat;
755 }
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775 public void luDecompose() throws InvalidMatrixException {
776
777 int nRows = this.getRowDimension();
778 int nCols = this.getColumnDimension();
779 if (nRows != nCols) {
780 throw new InvalidMatrixException("LU decomposition requires that the matrix be square.");
781 }
782 lu = this.getData();
783
784
785 permutation = new int[nRows];
786 for (int row = 0; row < nRows; row++) {
787 permutation[row] = row;
788 }
789 parity = 1;
790
791
792 for (int col = 0; col < nCols; col++) {
793
794 double sum = 0;
795
796
797 for (int row = 0; row < col; row++) {
798 sum = lu[row][col];
799 for (int i = 0; i < row; i++) {
800 sum -= lu[row][i] * lu[i][col];
801 }
802 lu[row][col] = sum;
803 }
804
805
806 int max = col;
807 double largest = 0d;
808 for (int row = col; row < nRows; row++) {
809 sum = lu[row][col];
810 for (int i = 0; i < col; i++) {
811 sum -= lu[row][i] * lu[i][col];
812 }
813 lu[row][col] = sum;
814
815
816 if (Math.abs(sum) > largest) {
817 largest = Math.abs(sum);
818 max = row;
819 }
820 }
821
822
823 if (Math.abs(lu[max][col]) < TOO_SMALL) {
824 lu = null;
825 throw new InvalidMatrixException("matrix is singular");
826 }
827
828
829 if (max != col) {
830 double tmp = 0;
831 for (int i = 0; i < nCols; i++) {
832 tmp = lu[max][i];
833 lu[max][i] = lu[col][i];
834 lu[col][i] = tmp;
835 }
836 int temp = permutation[max];
837 permutation[max] = permutation[col];
838 permutation[col] = temp;
839 parity = -parity;
840 }
841
842
843 for (int row = col + 1; row < nRows; row++) {
844 lu[row][col] /= lu[col][col];
845 }
846 }
847 }
848
849
850
851
852
853 public String toString() {
854 StringBuffer res = new StringBuffer();
855 res.append("RealMatrixImpl{");
856 if (data != null) {
857 for (int i = 0; i < data.length; i++) {
858 if (i > 0)
859 res.append(",");
860 res.append("{");
861 for (int j = 0; j < data[0].length; j++) {
862 if (j > 0)
863 res.append(",");
864 res.append(data[i][j]);
865 }
866 res.append("}");
867 }
868 }
869 res.append("}");
870 return res.toString();
871 }
872
873
874
875
876
877
878
879
880
881
882 public boolean equals(Object object) {
883 if (object == this ) {
884 return true;
885 }
886 if (object instanceof RealMatrixImpl == false) {
887 return false;
888 }
889 RealMatrix m = (RealMatrix) object;
890 int nRows = getRowDimension();
891 int nCols = getColumnDimension();
892 if (m.getColumnDimension() != nCols || m.getRowDimension() != nRows) {
893 return false;
894 }
895 for (int row = 0; row < nRows; row++) {
896 for (int col = 0; col < nCols; col++) {
897 if (Double.doubleToLongBits(data[row][col]) !=
898 Double.doubleToLongBits(m.getEntry(row, col))) {
899 return false;
900 }
901 }
902 }
903 return true;
904 }
905
906
907
908
909
910
911 public int hashCode() {
912 int ret = 7;
913 int nRows = getRowDimension();
914 int nCols = getColumnDimension();
915 ret = ret * 31 + nRows;
916 ret = ret * 31 + nCols;
917 for (int row = 0; row < nRows; row++) {
918 for (int col = 0; col < nCols; col++) {
919 ret = ret * 31 + (11 * (row+1) + 17 * (col+1)) *
920 MathUtils.hash(data[row][col]);
921 }
922 }
923 return ret;
924 }
925
926
927
928
929
930
931
932
933
934
935
936 protected RealMatrix getIdentity(int dimension) {
937 return MatrixUtils.createRealIdentityMatrix(dimension);
938 }
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967 protected RealMatrix getLUMatrix() throws InvalidMatrixException {
968 if (lu == null) {
969 luDecompose();
970 }
971 return new RealMatrixImpl(lu);
972 }
973
974
975
976
977
978
979
980
981
982
983
984
985
986 protected int[] getPermutation() {
987 int[] out = new int[permutation.length];
988 System.arraycopy(permutation, 0, out, 0, permutation.length);
989 return out;
990 }
991
992
993
994
995
996
997
998
999 private double[][] copyOut() {
1000 int nRows = this.getRowDimension();
1001 double[][] out = new double[nRows][this.getColumnDimension()];
1002
1003 for (int i = 0; i < nRows; i++) {
1004 System.arraycopy(data[i], 0, out[i], 0, data[i].length);
1005 }
1006 return out;
1007 }
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019 private void copyIn(double[][] in) {
1020 setSubMatrix(in,0,0);
1021 }
1022
1023
1024
1025
1026
1027
1028
1029
1030 private boolean isValidCoordinate(int row, int col) {
1031 int nRows = this.getRowDimension();
1032 int nCols = this.getColumnDimension();
1033
1034 return !(row < 0 || row > nRows - 1 || col < 0 || col > nCols -1);
1035 }
1036
1037 }