1 package org.catacomb.numeric.math;
2
3
4 public final class Matrix extends Object implements Cloneable {
5
6 public double a[][];
7 double ws[];
8 int n;
9 int perm[];
10 int sign;
11
12 int n1 = 0;
13 int n2 = 0;
14
15
16 public Matrix(int nnn) {
17 n1 = n2 = n = nnn;
18 a = new double[n][n];
19 perm = new int[n];
20 ws = new double[n];
21 }
22
23
24 public Matrix(double[][] a) {
25 this.a = a;
26 n1 = n2 = n = a.length;
27 perm = new int[n];
28 ws = new double[n];
29 }
30
31
32 public double[] flatten() {
33 double[] d = new double[n1 * n2];
34 for (int i = 0; i < n; i++) {
35 for (int j = 0; j < n; j++) {
36 d[n * i + j] = a[i][j];
37 }
38 }
39 return d;
40 }
41
42
43 public void Sp(String s) {
44 System.out.println(s);
45 }
46
47
48 public final int dim() {
49 return n;
50 }
51
52
53 public Matrix copy() {
54 Matrix m = new Matrix(n);
55
56 for (int i = 0; i < n; i++) {
57 for (int j = 0; j < n; j++) {
58 m.a[i][j] = a[i][j];
59 }
60 m.perm[i] = perm[i];
61 }
62
63 m.sign = sign;
64 m.setDims(n1, n2);
65 return m;
66 }
67
68
69 public void setDims(int d1, int d2) {
70 n1 = d1;
71 n2 = d2;
72 }
73
74
75 public void identise() {
76 for (int i = 0; i < n; i++) {
77 for (int j = 0; j < n; j++) {
78 a[i][j] = 0.;
79 }
80 a[i][i] = 1.;
81 }
82 }
83
84
85 public void randomise() {
86 for (int i = 0; i < n; i++) {
87 for (int j = 0; j < n; j++) {
88 a[i][j] = Math.random();
89 }
90 }
91 }
92
93
94 public void zero() {
95 for (int i = 0; i < n; i++) {
96 for (int j = 0; j < n; j++) {
97 a[i][j] = 0.;
98 }
99 }
100 }
101
102
103
104 public Matrix identity() {
105 Matrix m = copy();
106 m.identise();
107 return m;
108 }
109
110
111 public Matrix random() {
112 Matrix m = copy();
113 m.randomise();
114 return m;
115 }
116
117
118
119 public void add(double d) {
120 for (int i = 0; i < n; i++) {
121 for (int j = 0; j < n; j++) {
122 a[i][j] += d;
123 }
124 }
125 }
126
127
128
129 public static Matrix[] average(Matrix[] ma, Matrix[] mb, double f) {
130 int n = ma.length;
131 Matrix[] res = new Matrix[n];
132 for (int i = 0; i < n; i++) {
133 res[i] = average(ma[i], mb[i], f);
134 }
135 return res;
136 }
137
138
139 public static Matrix average(Matrix ma, Matrix mb, double f) {
140 double g = 1. - f;
141 int n = ma.n;
142 Matrix res = new Matrix(n);
143 for (int i = 0; i < n; i++) {
144 for (int j = 0; j < n; j++) {
145 res.a[i][j] = f * mb.a[i][j] + g * ma.a[i][j];
146 }
147 }
148 return res;
149 }
150
151
152
153 public void add(Matrix m) {
154 if (m.n != n) {
155 Sp("incompativle dims in Matrix.mplyBy " + n + " " + m.n);
156
157 } else {
158 for (int i = 0; i < n; i++) {
159 for (int j = 0; j < n; j++) {
160 a[i][j] += m.a[i][j];
161 }
162 }
163 }
164 }
165
166
167 public Matrix sum(Matrix m) {
168 Matrix mr = copy();
169 mr.zero();
170 if (m.n != n) {
171 Sp("incompativle dims in Matrix.mplyBy " + n + " " + m.n);
172
173 } else {
174 for (int i = 0; i < n; i++) {
175 for (int j = 0; j < n; j++) {
176 mr.a[i][j] = a[i][j] + m.a[i][j];
177 }
178 }
179 }
180 return mr;
181 }
182
183
184 public void mpyBy(double d) {
185 for (int i = 0; i < n; i++) {
186 for (int j = 0; j < n; j++) {
187 a[i][j] *= d;
188 }
189 }
190 }
191
192
193 public void mpyBy(Matrix m) {
194 a = (prod(m)).a;
195 }
196
197
198 public Matrix prod(Matrix m) {
199 Matrix mr = copy();
200 mr.zero();
201
202 if (m.n != n) {
203 Sp("incompativle dims in Matrix.mplyBy " + n + " " + m.n);
204
205 } else {
206 for (int i = 0; i < n; i++) {
207 for (int j = 0; j < n; j++) {
208 for (int k = 0; k < n; k++) {
209 mr.a[i][j] += a[i][k] * m.a[k][j];
210 }
211 }
212 }
213 }
214 return mr;
215 }
216
217
218 public double[] lvprod(double[] v) {
219 double[] r = new double[n];
220 if (v.length != n) {
221 Sp("incompatible dimensions in lvprod");
222 } else {
223 for (int i = 0; i < n; i++) {
224 for (int j = 0; j < n; j++) {
225 r[j] += v[i] * a[i][j];
226 }
227 }
228 }
229 return r;
230 }
231
232
233 public double[] rvprod(double[] v) {
234 double[] r = new double[n];
235 if (v.length != n) {
236 Sp("incompatible dimensions in lvprod");
237 } else {
238 for (int i = 0; i < n; i++) {
239 for (int j = 0; j < n; j++) {
240 r[i] += a[i][j] * v[j];
241 }
242 }
243 }
244 return r;
245 }
246
247
248
249 public void multiplyInto(double[] v) {
250 for (int i = 0; i < n; i++) {
251 ws[i] = 0.;
252 for (int j = 0; j < n; j++) {
253 ws[i] += a[i][j] * v[j];
254 }
255 }
256 for (int i = 0; i < n; i++) {
257 v[i] = ws[i];
258 }
259 }
260
261
262 public void rect2rvprod(double[] v, double[] r1, double[] r2) {
263 for (int i = 0; i < n1; i++) {
264 ws[i] = 0.0;
265 for (int j = 0; j < n1; j++) {
266 ws[i] += a[i][j] * v[j];
267 }
268 }
269
270 for (int i = 0; i < n2 - n1; i++) {
271 r2[i] = 0.0;
272 for (int j = 0; j < n1; j++) {
273 r2[i] += a[i + n1][j] * v[j];
274 }
275 }
276 for (int i = 0; i < n1; i++) {
277 r1[i] = ws[i];
278 }
279
280
281 }
282
283
284
285 public double rvprodOneElt(double[] v, int elt) {
286 double r = 0.0;
287 for (int j = 0; j < n; j++) {
288 r += a[elt][j] * v[j];
289 }
290 return r;
291 }
292
293
294 public Matrix transpose() {
295 Matrix m = copy();
296 for (int i = 0; i < n; i++) {
297 for (int j = 0; j < n; j++) {
298 m.a[i][j] = a[j][i];
299 }
300 }
301 return m;
302 }
303
304
305 public double det() {
306 Matrix t = copy();
307 t.LU();
308 double d = 1.0 * t.sign;
309 for (int i = 0; i < n; i++) {
310 d *= t.a[i][i];
311 }
312 return d;
313 }
314
315
316
317 public double[][] copyMat() {
318 double[][] ar = new double[n][n];
319 for (int i = 0; i < n; i++) {
320 for (int j = 0; j < n; j++) {
321 ar[i][j] = a[i][j];
322 }
323 }
324 return ar;
325 }
326
327
328
329 public void LU() {
330 int i, imax, j, k;
331 double big, dum, sum, temp;
332 double vv[] = new double[n];
333 double TINY = 1.0e-20;
334
335 sign = 1;
336
337 imax = -1;
338 for (i = 0; i < n; i++) {
339 big = 0.0;
340 for (j = 0; j < n; j++) {
341 if ((temp = Math.abs(a[i][j])) > big) {
342 big = temp;
343 }
344 }
345 if (big == 0.0) {
346 Sp("Singular Matrix in routine LUDCMP");
347 }
348 vv[i] = 1.0 / big;
349 }
350
351 for (j = 0; j < n; j++) {
352 for (i = 0; i < j; i++) {
353 sum = a[i][j];
354 for (k = 0; k < i; k++) {
355 sum -= a[i][k] * a[k][j];
356 }
357 a[i][j] = sum;
358 }
359 big = 0.0;
360 for (i = j; i < n; i++) {
361 sum = a[i][j];
362 for (k = 0; k < j; k++) {
363 sum -= a[i][k] * a[k][j];
364 }
365 a[i][j] = sum;
366 if ((dum = vv[i] * Math.abs(sum)) >= big) {
367 big = dum;
368 imax = i;
369 }
370 }
371 if (j != imax) {
372 for (k = 0; k < n; k++) {
373 dum = a[imax][k];
374 a[imax][k] = a[j][k];
375 a[j][k] = dum;
376 }
377 sign = -sign;
378 vv[imax] = vv[j];
379 }
380 perm[j] = imax;
381 if (a[j][j] == 0.0) {
382 a[j][j] = TINY;
383 }
384 if (j != n) {
385 dum = 1.0 / (a[j][j]);
386 for (i = j + 1; i < n; i++) {
387 a[i][j] *= dum;
388 }
389 }
390 }
391 }
392
393
394
395 public Matrix inverse() {
396 Matrix t, r;
397 t = copy();
398 r = copy();
399 t.LU();
400
401 double[] c = new double[n];
402 for (int j = 0; j < n; j++) {
403 for (int i = 0; i < n; i++) {
404 c[i] = 0.0;
405 }
406 c[j] = 1.0;
407 t.lubksb(c);
408 for (int i = 0; i < n; i++) {
409 r.a[i][j] = c[i];
410 }
411 }
412 return r;
413 }
414
415
416
417 public static double[] LUSolve(double[][] m, double[] R) {
418
419 Matrix M = new Matrix(m);
420 M.LU();
421 double[] W = M.lubksb(R);
422 return W;
423 }
424
425
426
427 public void invert() {
428 a = (inverse()).a;
429 }
430
431
432
433 public double[] lubksb(double[] b) {
434 int ip;
435 int ii = -1;
436 double sum;
437
438 for (int i = 0; i < n; i++) {
439 ip = perm[i];
440 sum = b[ip];
441 b[ip] = b[i];
442 if (ii >= 0) {
443 for (int j = ii; j < i; j++) {
444 sum -= a[i][j] * b[j];
445 }
446 } else if (sum != 0.0) {
447 ii = i;
448 }
449 b[i] = sum;
450 }
451
452 for (int i = n - 1; i >= 0; i--) {
453 sum = b[i];
454 for (int j = i + 1; j < n; j++) {
455 sum -= a[i][j] * b[j];
456 }
457 b[i] = sum / a[i][i];
458 }
459 return b;
460 }
461
462
463 public void round(double d) {
464 for (int i = 0; i < n; i++) {
465 for (int j = 0; j < n; j++) {
466 if (Math.abs(a[i][j]) < d) {
467 a[i][j] = 0.0;
468 }
469 }
470 }
471 }
472
473
474
475 public void round() {
476 round(1.0e-15);
477 }
478
479
480
481 public void print() {
482
483 String[] sa = new String[n];
484 for (int i = 0; i < n; i++) {
485 sa[i] = " ";
486 for (int j = 0; j < n; j++) {
487 sa[i] += (" " + a[i][j]);
488 }
489 }
490 Sp(" n1: " + n1 + " n2: " + n2);
491 for (int i = 0; i < sa.length; i++) {
492 Sp("" + i + " " + sa[i]);
493 }
494 }
495
496
497
498 public double maxAbsElt() {
499 double d = 0.0;
500 for (int i = 0; i < n; i++) {
501 for (int j = 0; j < n; j++) {
502 if (Math.abs(a[i][j]) > d) {
503 d = Math.abs(a[i][j]);
504 }
505 }
506 }
507 return d;
508 }
509
510
511
512 public Matrix power(int p) {
513 Matrix mr = identity();
514 int pg = 0;
515 int pl = 1;
516 Matrix mp = copy();
517 while (pg < p) {
518 if ((p & pl) > 0) {
519 pg += pl;
520 mr = mr.prod(mp);
521 }
522 pl *= 2;
523 mp = mp.prod(mp);
524 }
525 if (pg != p) {
526 Sp("got Matrix power wrong: " + p + " " + pg + " " + pl);
527 }
528 return mr;
529 }
530
531
532 public Matrix crudeExpOf(double t) {
533 Matrix m = copy();
534 m.mpyBy(t);
535 double eps = 1.0e-8;
536
537 double d = m.maxAbsElt();
538 int p = 0;
539 double f = 1;
540 for (p = 0; d * f > eps; f *= 0.5, p++) {
541 ;
542 }
543
544 m.mpyBy(f);
545 m.add(m.identity());
546
547 for (; p > 0; p--) {
548 m.mpyBy(m);
549 }
550 return m;
551 }
552
553
554
555 public Matrix expOf(double t) {
556 Matrix m = copy();
557 m.mpyBy(t);
558 double eps = 1.0e-12;
559
560 double d = m.maxAbsElt();
561 int p = 0;
562 double f = 1;
563 for (p = 0; d * f > eps; f *= 0.5, p++) {
564 ;
565 }
566
567 m.mpyBy(f);
568
569
570
571 for (; p > 0; p--) {
572 Matrix u = m.copy();
573 u.mpyBy(u);
574 m.add(m);
575 m.add(u);
576 }
577 m.add(m.identity());
578 return m;
579 }
580
581
582
583 public int randomIndexFromColumn(int c) {
584 return randomIndexFromColumn(c, Math.random());
585 }
586
587
588 public final int randomIndexFromColumn(int c, double rin) {
589 double r = rin;
590 int ir = 0;
591 while ((r -= a[ir][c]) > 0)
592 ir++;
593 return ir;
594 }
595
596
597
598
599 public final double[] getColumn(int ic) {
600 double[] c = new double[n2];
601 for (int i = 0; i < n2; i++) {
602 c[i] = a[i][ic];
603 }
604 return c;
605 }
606
607
608
609 public int randomIndexFromOffsetColumn(int c, int off) {
610 double r = Math.random();
611 int ir;
612 for (ir = off; (r -= a[ir][c]) > 0; ir++) {
613 ;
614 }
615 return ir;
616 }
617
618
619
620 public double[] ev1vec(int np) {
621
622
623
624
625
626 Matrix q = copy();
627 for (int i = 0; i < np; i++) {
628 q = q.prod(q);
629 }
630
631 double[] s = new double[n];
632 for (int i = 0; i < n; i++) {
633 s[i] = 1. / n;
634 }
635 s = q.rvprod(s);
636
637 double t = 0.0;
638 for (int i = 0; i < n; i++) {
639 t += s[i];
640 }
641 for (int i = 0; i < n; i++) {
642 s[i] /= t;
643 }
644 if (Math.abs(t - 1.) > 0.01) {
645 Sp("WARNING - ev1vec in class Matrix chnaged size " + t);
646 }
647 return s;
648 }
649
650
651 }