1 package org.lcsim.spacegeom;
2 import java.io.*;
3
4
5
6
7
8 class Matrix
9 {
10
11 protected int _row;
12 protected int _col;
13 protected double _a[][];
14 protected Matrix _eigenvalues;
15 protected Eigensystem _eigensystem;
16
17 public Matrix()
18 {
19 _row = 0;
20 _col = 0;
21 createData();
22 }
23
24 public Matrix(int i, int j)
25 {
26 _row = i;
27 _col = j;
28 createData();
29 }
30
31 public Matrix(Matrix matrix)
32 {
33 _row = matrix._row;
34 _col = matrix._col;
35 _eigensystem = matrix._eigensystem;
36 createData();
37 for(int i = 0; i < _row; i++)
38 {
39 for(int j = 0; j < _col; j++)
40 _a[i][j] = matrix._a[i][j];
41
42 }
43
44 }
45
46 protected void QL(Matrix matrix, Matrix matrix1, Matrix matrix2, boolean flag)
47 {
48 int k1 = _row;
49 for(int i = 1; i < k1; i++)
50 matrix2.set(i - 1, 0, matrix2.at(i, 0));
51
52 matrix2.set(k1 - 1, 0, 0.0D);
53 for(int l = 0; l < k1; l++)
54 {
55 int j1 = 0;
56 int i1;
57 do
58 {
59 for(i1 = l; i1 < k1 - 1; i1++)
60 {
61 double d2 = Math.abs(matrix1.at(i1, 0)) + Math.abs(matrix1.at(i1 + 1, 0));
62 if(Math.abs(matrix2.at(i1, 0)) + d2 == d2)
63 break;
64 }
65
66 if(i1 != l)
67 {
68 if(j1++ == 30)
69 {
70 System.err.println("Matrix::QL -> Too many iterations.");
71 System.exit(0);
72 }
73 double d5 = (matrix1.at(l + 1, 0) - matrix1.at(l, 0)) / (2D * matrix2.at(l, 0));
74 double d7 = pythagoras(d5, 1.0D);
75 d5 = (matrix1.at(i1, 0) - matrix1.at(l, 0)) + matrix2.at(l, 0) / (d5 + SIGN(d7, d5));
76 double d1;
77 double d8 = d1 = 1.0D;
78 double d6 = 0.0D;
79 int j;
80 for(j = i1 - 1; j >= l; j--)
81 {
82 double d3 = d8 * matrix2.at(j, 0);
83 double d = d1 * matrix2.at(j, 0);
84 matrix2.set(j + 1, 0, d7 = pythagoras(d3, d5));
85 if(d7 == 0.0D)
86 {
87 matrix1.set(j + 1, 0, matrix1.at(j + 1, 0) - d6);
88 matrix2.set(i1, 0, 0.0D);
89 break;
90 }
91 d8 = d3 / d7;
92 d1 = d5 / d7;
93 d5 = matrix1.at(j + 1, 0) - d6;
94 d7 = (matrix1.at(j, 0) - d5) * d8 + 2D * d1 * d;
95 matrix1.set(j + 1, 0, d5 + (d6 = d8 * d7));
96 d5 = d1 * d7 - d;
97 if(flag)
98 {
99 for(int k = 0; k < k1; k++)
100 {
101 double d4 = matrix.at(k, j + 1);
102 matrix.set(k, j + 1, d8 * matrix.at(k, j) + d1 * d4);
103 matrix.set(k, j, d1 * matrix.at(k, j) - d8 * d4);
104 }
105
106 }
107 }
108
109 if(d7 != 0.0D || j < l)
110 {
111 matrix1.set(l, 0, matrix1.at(l, 0) - d6);
112 matrix2.set(l, 0, d5);
113 matrix2.set(i1, 0, 0.0D);
114 }
115 }
116 } while(i1 != l);
117 }
118
119 }
120
121 protected double SIGN(double d, double d1)
122 {
123 if(d1 >= 0.0D)
124 return Math.abs(d);
125 else
126 return -Math.abs(d);
127 }
128
129
130 public double at(int i, int j)
131 {
132 return _a[i][j];
133 }
134
135 public int cols()
136 {
137 return _col;
138 }
139
140 protected void createData()
141 {
142 int i = _row + 1;
143 int j = _col + 1;
144 _a = new double[i][j];
145 }
146
147 public Eigensystem eigensystem()
148 {
149 if(!isSymmetric())
150 {
151 System.err.println("Matrix::eigensystem -> Matrix is not symmetric.");
152 System.exit(0);
153 }
154 int i = _row;
155 Matrix matrix = new Matrix(this);
156 Matrix matrix1 = new Matrix(i, 1);
157 Matrix matrix2 = new Matrix(i, 1);
158 householder(matrix, matrix1, matrix2, true);
159 QL(matrix, matrix1, matrix2, true);
160 if(_eigensystem == null)
161 _eigensystem = new Eigensystem();
162 Matrix matrix3 = new Matrix(matrix);
163 _eigensystem.setEigenvectors(matrix3);
164 Matrix matrix4 = new Matrix(matrix1);
165 _eigensystem.setEigenvalues(matrix4);
166 return _eigensystem;
167 }
168
169 public Matrix eigenvalues()
170 {
171 if(!isSymmetric())
172 {
173 System.err.println("Matrix::eigenvalues -> Matrix is not symmetric.");
174 System.exit(0);
175 }
176 int i = _row;
177 Matrix matrix = new Matrix(this);
178 Matrix matrix1 = new Matrix(i, 1);
179 Matrix matrix2 = new Matrix(i, 1);
180 householder(matrix, matrix1, matrix2, false);
181 QL(matrix, matrix1, matrix2, false);
182 if(_eigenvalues == null)
183 _eigenvalues = new Matrix(i, 1);
184 _eigenvalues = new Matrix(matrix1);
185 return _eigenvalues;
186 }
187
188
189 public Eigensystem getEigensystem()
190 {
191 return _eigensystem;
192 }
193
194 public Matrix getEigenvalues()
195 {
196 return _eigenvalues;
197 }
198
199 protected void householder(Matrix matrix, Matrix matrix1, Matrix matrix2, boolean flag)
200 {
201 int j3 = 0;
202 int k3 = _row;
203 for(int i = k3 - 1; i > 0; i--)
204 {
205 j3 = i - 1;
206 double d;
207 double d7 = d = 0.0D;
208 if(j3 > 0)
209 {
210 for(int k1 = 0; k1 <= j3; k1++)
211 d += Math.abs(matrix.at(i, k1));
212
213 if(d == 0.0D)
214 {
215 matrix2.set(i, 0, matrix.at(i, i));
216 }
217 else
218 {
219 for(int l1 = 0; l1 <= j3; l1++)
220 {
221 matrix.set(i, l1, matrix.at(i, l1) / d);
222 d7 += matrix.at(i, l1) * matrix.at(i, l1);
223 }
224
225 double d1 = matrix.at(i, j3);
226 double d3 = d1 < 0.0D ? Math.sqrt(d7) : -Math.sqrt(d7);
227 matrix2.set(i, 0, d * d3);
228 d7 -= d1 * d3;
229 matrix.set(i, j3, d1 - d3);
230 d1 = 0.0D;
231 for(int k = 0; k <= j3; k++)
232 {
233 if(flag)
234 matrix.set(k, i, matrix.at(i, k) / d7);
235 double d4 = 0.0D;
236 for(int i2 = 0; i2 <= k; i2++)
237 d4 += matrix.at(k, i2) * matrix.at(i, i2);
238
239 for(int j2 = k + 1; j2 <= j3; j2++)
240 d4 += matrix.at(j2, k) * matrix.at(i, j2);
241
242 matrix2.set(k, 0, d4 / d7);
243 d1 += matrix2.at(k, 0) * matrix.at(i, k);
244 }
245
246 double d8 = d1 / (d7 + d7);
247 for(int l = 0; l <= j3; l++)
248 {
249 double d2 = matrix.at(i, l);
250 double d5 = matrix2.at(l, 0) - d8 * d2;
251 matrix2.set(l, 0, d5);
252 for(int k2 = 0; k2 <= l; k2++)
253 matrix.set(l, k2, matrix.at(l, k2) - (d2 * matrix2.at(k2, 0) + d5 * matrix.at(i, k2)));
254
255 }
256
257 }
258 }
259 else
260 {
261 matrix2.set(i, 0, matrix.at(i, j3));
262 }
263 matrix1.set(i, 0, d7);
264 }
265
266 if(flag)
267 matrix1.set(0, 0, 0.0D);
268 matrix2.set(0, 0, 0.0D);
269 for(int j = 0; j < k3; j++)
270 {
271 if(flag)
272 {
273 j3 = j - 1;
274 if(matrix1.at(j, 0) != 0.0D)
275 {
276 for(int i1 = 0; i1 <= j3; i1++)
277 {
278 double d6 = 0.0D;
279 for(int l2 = 0; l2 <= j3; l2++)
280 d6 += matrix.at(j, l2) * matrix.at(l2, i1);
281
282 for(int i3 = 0; i3 <= j3; i3++)
283 matrix.set(i3, i1, matrix.at(i3, i1) - d6 * matrix.at(i3, j));
284
285 }
286
287 }
288 }
289 matrix1.set(j, 0, matrix.at(j, j));
290 if(flag)
291 {
292 matrix.set(j, j, 1.0D);
293 for(int j1 = 0; j1 <= j3; j1++)
294 {
295 matrix.set(j, j1, 0.0D);
296 matrix.set(j1, j, matrix.at(j, j1));
297 }
298
299 }
300 }
301
302 }
303
304 public Matrix inverse()
305 {
306 if(!isSquare())
307 {
308 System.err.println("Matrix::inverse -> Matrix is not square.");
309 System.exit(0);
310 }
311 int i = _row;
312 Matrix matrix = new Matrix(i, 2 * i);
313 Matrix matrix1 = new Matrix(i, 1);
314 Matrix matrix2 = new Matrix(i, 1);
315 for(int j = 0; j < i; j++)
316 {
317 for(int i2 = 0; i2 < i; i2++)
318 matrix.set(j, i2, _a[j][i2]);
319
320 }
321
322 for(int k = 0; k < i; k++)
323 {
324 for(int j2 = i; j2 < 2 * i; j2++)
325 matrix.set(k, j2, 0.0D);
326
327 matrix.set(k, i + k, 1.0D);
328 }
329
330 for(int i4 = 0; i4 < i; i4++)
331 {
332 for(int l = i4; l < i; l++)
333 {
334 matrix2.set(l, 0, 0.0D);
335 for(int k2 = i4; k2 < i; k2++)
336 matrix2.set(l, 0, matrix2.at(l, 0) + matrix.at(l, k2));
337
338 matrix1.set(l, 0, Math.abs(matrix.at(l, i4)) / matrix2.at(l, 0));
339 }
340
341 int k4 = i4;
342 for(int i1 = i4 + 1; i1 < i; i1++)
343 if(matrix1.at(i1, 0) > matrix1.at(i1 - 1, 0))
344 k4 = i1;
345
346 if(matrix1.at(k4, 0) == 0.0D)
347 {
348 System.err.println("Matrix::inverse -> Matrix is singular.");
349 System.exit(0);
350 }
351 if(k4 != i4)
352 {
353 for(int l2 = 0; l2 < 2 * i; l2++)
354 {
355 double d = matrix.at(i4, l2);
356 matrix.set(i4, l2, matrix.at(k4, l2));
357 matrix.set(k4, l2, d);
358 }
359
360 }
361 for(int i3 = 2 * i - 1; i3 >= i4; i3--)
362 matrix.set(i4, i3, matrix.at(i4, i3) / matrix.at(i4, i4));
363
364 for(int j1 = i4 + 1; j1 < i; j1++)
365 {
366 for(int j3 = 2 * i - 1; j3 >= i4 + 1; j3--)
367 matrix.set(j1, j3, matrix.at(j1, j3) - matrix.at(i4, j3) * matrix.at(j1, i4));
368
369 }
370
371 }
372
373 for(int j4 = i - 1; j4 >= 0; j4--)
374 {
375 for(int k1 = j4 - 1; k1 >= 0; k1--)
376 {
377 for(int k3 = i; k3 < 2 * i; k3++)
378 matrix.set(k1, k3, matrix.at(k1, k3) - matrix.at(j4, k3) * matrix.at(k1, j4));
379
380 }
381
382 }
383
384 Matrix matrix3 = new Matrix(i, i);
385 for(int l1 = 0; l1 < i; l1++)
386 {
387 for(int l3 = 0; l3 < i; l3++)
388 matrix3.set(l1, l3, matrix.at(l1, l3 + i));
389
390 }
391
392 return matrix3;
393 }
394
395 public boolean isLowerDiagonal()
396 {
397 for(int i = 0; i < _row - 1; i++)
398 {
399 for(int j = i + 1; j < _col; j++)
400 if(_a[i][j] != 0.0D)
401 return false;
402
403 }
404
405 return true;
406 }
407
408 public boolean isSquare()
409 {
410 return _row == _col;
411 }
412
413 public boolean isSymmetric()
414 {
415 if(!isSquare())
416 return false;
417 for(int i = 0; i < _row; i++)
418 {
419 for(int j = i + 1; j < _col; j++)
420 if(_a[i][j] != _a[j][i])
421 return false;
422
423 }
424
425 return true;
426 }
427
428 public boolean isUpperDiagonal()
429 {
430 for(int i = 1; i < _row; i++)
431 {
432 for(int j = 0; j < i; j++)
433 if(_a[i][j] != 0.0D)
434 return false;
435
436 }
437
438 return true;
439 }
440
441 public void makeIdentity()
442 {
443 if(_row != _col)
444 {
445 System.err.println("Matrix::makeIdentity -> Matrix is not square.");
446 System.exit(0);
447 }
448 int i = _row;
449 for(int j = 0; j < i; j++)
450 {
451 for(int k = 0; k < i; k++)
452 _a[j][k] = 0.0D;
453
454 _a[j][j] = 1.0D;
455 }
456
457 }
458
459 public Matrix times(double d)
460 {
461 int i = _row;
462 int j = _col;
463 Matrix matrix = new Matrix(i, j);
464 for(int k = 0; k < i; k++)
465 {
466 for(int l = 0; l < j; l++)
467 matrix.set(k, l, d * _a[k][l]);
468
469 }
470
471 return matrix;
472 }
473
474 public Matrix times(Matrix matrix)
475 {
476 int i = _row;
477 int j = _col;
478 int k = matrix.cols();
479 if(_col != matrix.rows())
480 {
481 System.err.println("Matrix::operator* -> Invalid dimensions of matrices.");
482 System.exit(0);
483 }
484 Matrix matrix1 = new Matrix(i, k);
485 for(int l = 0; l < i; l++)
486 {
487 for(int i1 = 0; i1 < k; i1++)
488 {
489 double d = 0.0D;
490 for(int j1 = 0; j1 < j; j1++)
491 d += _a[l][j1] * matrix.at(j1, i1);
492
493 matrix1.set(l, i1, d);
494 }
495
496 }
497
498 return matrix1;
499 }
500
501 protected double pythagoras(double d, double d1)
502 {
503 double d2 = Math.abs(d);
504 double d3 = Math.abs(d1);
505 if(d2 > d3)
506 return d2 * Math.sqrt(1.0D + square(d3 / d2));
507 else
508 return d3 != 0.0D ? d3 * Math.sqrt(1.0D + square(d2 / d3)) : 0.0D;
509 }
510
511
512
513 public int rows()
514 {
515 return _row;
516 }
517
518 public void set(int i, int j, double d)
519 {
520 _a[i][j] = d;
521 }
522
523 public void plusequal(int i, int j, double d)
524 {
525 _a[i][j] += d;
526 }
527
528 public void plusEqual(Matrix b)
529 {
530 for(int i=0; i<rows(); ++i)
531 {
532 for(int j=0; j<cols(); ++j)
533 {
534 _a[i][j] += b._a[i][j];
535 }
536 }
537 }
538
539 protected double square(double d)
540 {
541 if(d == 0.0D)
542 return 0.0D;
543 else
544 return d * d;
545 }
546
547 public Matrix submatrix(int i, int j, int k, int l)
548 {
549 Matrix matrix = new Matrix(i, j);
550 for(int i1 = 0; i1 < i; i1++)
551 {
552 for(int j1 = 0; j1 < j; j1++)
553 matrix._a[i1][j1] = _a[i1 + k][j1 + l];
554
555 }
556
557 return matrix;
558 }
559
560 public Matrix subtract(Matrix matrix)
561 {
562 int i = _row;
563 int j = _col;
564 if(i != matrix.rows() || j != matrix.cols())
565 {
566 System.err.println("Matrix::operator- -> Invalid dimensions of matrices.");
567 System.exit(0);
568 }
569 Matrix matrix1 = new Matrix(i, j);
570 for(int k = 0; k < i; k++)
571 {
572 for(int l = 0; l < j; l++)
573 matrix1.set(k, l, _a[k][l] - matrix.at(k, l));
574
575 }
576
577 return matrix1;
578 }
579
580 public Matrix transposed()
581 {
582 int i = _row;
583 int j = _col;
584 Matrix matrix = new Matrix(j, i);
585 for(int k = 0; k < j; k++)
586 {
587 for(int l = 0; l < i; l++)
588 matrix.set(k, l, _a[l][k]);
589
590 }
591
592 return matrix;
593 }
594
595 public Matrix upperDiagonal()
596 {
597 int i2 = _row;
598 int j2 = _col;
599 if(i2 < j2)
600 {
601 System.err.println("Matrix::upperDiagonal -> More rows than columns required.");
602 System.exit(0);
603 }
604 Matrix matrix = new Matrix(this);
605 int l1 = j2;
606 if(i2 - 1 < j2)
607 l1 = i2 - 1;
608 for(int i1 = 0; i1 <= l1; i1++)
609 {
610 double d1 = 0.0D;
611 for(int i = i1; i < i2; i++)
612 d1 += matrix.at(i, i1) * matrix.at(i, i1);
613
614 double d4 = 0.0D;
615 if(matrix.at(i1, i1) < 0.0D)
616 d4 = -1D;
617 if(matrix.at(i1, i1) > 0.0D)
618 d4 = 1.0D;
619 double d2 = -d4 * Math.sqrt(d1);
620 double d = d1 - matrix.at(i1, i1) * d2;
621 matrix.set(i1, i1, matrix.at(i1, i1) - d2);
622 for(int k1 = i1 + 1; k1 < j2; k1++)
623 {
624 double d3 = 0.0D;
625 for(int j = i1; j < i2; j++)
626 d3 += matrix.at(j, i1) * matrix.at(j, k1);
627
628 double d5 = d3 / d;
629 for(int k = i1; k < i2; k++)
630 matrix.set(k, k1, matrix.at(k, k1) - matrix.at(k, i1) * d5);
631
632 }
633
634 matrix.set(i1, i1, d2);
635 }
636
637 for(int l = 1; l < i2; l++)
638 {
639 for(int j1 = 0; j1 < l; j1++)
640 matrix.set(l, j1, 0.0D);
641
642 }
643
644 if(matrix.at(j2 - 1, j2 - 1) == 0.0D)
645 {
646 System.err.println("Matrix::householder -> Rank of matrix < n");
647 System.exit(0);
648 }
649 int k2 = _row;
650 if(_col < _row)
651 k2 = _col;
652 return matrix.submatrix(k2, k2, 0, 0);
653 }
654
655 public void write(String s)
656 {
657 write(s, "");
658 }
659
660 public void write(String s, String s1)
661 {
662 System.out.println("Writing matrix to '" + s + "'.");
663 try
664 {
665 PrintWriter printwriter = new PrintWriter(new BufferedWriter(new FileWriter(s)));
666 String s2 = "# " + s1;
667 printwriter.println(s2);
668 String s3 = _row + " " + _col;
669 printwriter.println(s3);
670 String s4 = "";
671 for(int i = 0; i < _row; i++)
672 {
673 String s5 = "";
674 for(int j = 0; j < _col; j++)
675 s5 = s5 + _a[i][j] + " ";
676
677 printwriter.println(s5.substring(0, s5.length() - 1));
678 }
679
680 printwriter.close();
681 }
682 catch(IOException _ex)
683 {
684 System.err.println("Matrix::write -> Error writing to '" + s + "'.");
685 System.exit(0);
686 }
687 System.out.println("Matrix written.");
688 }
689
690 public String toString()
691 {
692 return "Matrix: ";
693 }
694
695 public Matrix column(int col)
696 {
697 int rows = _row;
698 int cols = 1;
699 Matrix matrix = new Matrix(rows, cols);
700 for(int i = 0; i < rows; ++i)
701 {
702 matrix.set(i, 0, _a[i][col]);
703 }
704
705 return matrix;
706 }
707
708 public Matrix row(int row)
709 {
710 int rows = 1;
711 int cols = _col;
712 Matrix matrix = new Matrix(rows, cols);
713 for(int i = 0; i < cols; ++i)
714 {
715 matrix.set(0, i, _a[row][i]);
716 }
717
718 return matrix;
719 }
720 }