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