1
2
3
4
5
6
7
8
9
10
11
12
13 package no.uib.cipr.matrix;
14
15 import java.io.IOException;
16 import java.util.Arrays;
17 import no.uib.cipr.matrix.io.MatrixVectorReader;
18
19 import org.astrogrid.matrix.Matrix;
20 import org.netlib.lapack.LAPACK;
21 import org.netlib.util.intW;
22
23 public class AGDenseMatrix extends DenseMatrix implements Matrix {
24
25
26
27
28
29 public AGDenseMatrix(final MatrixVectorReader r) throws IOException {
30 super(r);
31
32 }
33
34
35
36
37 public AGDenseMatrix(final Vector x) {
38 super(x);
39
40 }
41
42
43
44
45 public AGDenseMatrix(final double[][] values) {
46 super(values);
47 }
48
49 public AGDenseMatrix(){
50 this(new double[][]{{}});
51 }
52
53
54
55
56
57 public AGDenseMatrix(final int numRows, final int numColumns) {
58 super(numRows, numColumns);
59 }
60
61
62
63
64
65 public AGDenseMatrix(final Matrix A, final boolean deep) {
66 super(A, deep);
67 data = null;
68 }
69
70 public static AGDenseMatrix repeatColumn(final Vector v, final int n) {
71 if (v instanceof DenseVector) {
72 final DenseVector dv = (DenseVector) v;
73
74
75 final AGDenseMatrix retval = new AGDenseMatrix(v.size(),n);
76 for (int i = 0; i < n; i++) {
77 System.arraycopy(dv.getData(), 0, retval.data, i*retval.numRows, retval.numRows);
78 }
79 return retval;
80 } else {
81 throw new UnsupportedOperationException("this method only works with DenseVectors at the moment");
82
83 }
84 }
85
86 public static AGDenseMatrix repeatRow(final Vector v, final int n){
87 final AGDenseMatrix retval = new AGDenseMatrix(n, v.size());
88 for (MatrixEntry e : retval) {
89 e.set(v.get(e.column()));
90 }
91 return retval;
92 }
93
94
95
96
97 public AGDenseMatrix(final no.uib.cipr.matrix.Matrix A) {
98 super(A);
99 }
100
101
102
103
104
105
106
107 public AGDenseMatrix(final Vector v, final int ni, final int nj) {
108 super(v.size()*ni, nj);
109
110 int i = 0;
111 for (VectorEntry vectorEntry : v) {
112 this.data[i++] = vectorEntry.get();
113 }
114
115 final int size = v.size();
116 for (i = 1; i < ni * nj ; i++){
117 System.arraycopy(data, 0, data, i*size, size);
118 }
119 }
120
121
122
123
124
125
126
127 public AGDenseMatrix(final Matrix m, final int ni, final int nj) {
128 super(m.numRows()* ni, m.numColumns()*nj);
129
130 for (MatrixEntry me : m) {
131 this.set(me.row(), me.column(), me.get());
132 }
133
134
135 for(int j = 0; j < m.numColumns(); j++){
136 for (int i = 1; i < ni; i++) {
137 final int src = j*m.numRows()*ni;
138 System.arraycopy(data, src, data, src + i * m.numRows(), m.numRows());
139 }
140 }
141
142 final int blocksize = m.numRows()*ni*m.numColumns();
143 for (int j=1; j < nj; j++){
144 System.arraycopy(data, 0, data, j*blocksize, blocksize);
145 }
146 }
147
148
149
150
151
152
153
154 public AGDenseMatrix(int numRows, int numColumns, double[] indata) {
155 super(numRows, numColumns);
156 if (numColumns* numRows != indata.length) {
157 throw new IllegalArgumentException("datasize not equal to numrows*numcolumns");
158 }
159 this.data = indata;
160 }
161
162 public Matrix slice(final int rowstart, final int rowend, final int colstart, final int colend) {
163
164 throw new UnsupportedOperationException(
165 "Matrix.slice() not implemented");
166 }
167
168
169
170
171
172
173
174
175
176
177 public Vector slicev(final int idx, final boolean dorow) {
178 if (dorow) {
179 if (idx < numRows) {
180 final Vector retval = new DenseVector(numColumns);
181 for (int j = 0; j < numColumns; j++) {
182 retval.set(j, this.get(idx, j));
183 }
184 return retval;
185 } else {
186 throw new IllegalArgumentException("row index " + idx
187 + " greater than " + numRows);
188 }
189 } else {
190 if (idx < numColumns) {
191 final double[] vals = new double[numRows];
192 System.arraycopy(data, idx * numRows, vals, 0, numRows);
193 final Vector retval = new DenseVector(vals, false);
194 return retval;
195 } else {
196 throw new IllegalArgumentException("column index" + idx
197 + "greater than " + numColumns);
198 }
199
200 }
201 }
202
203 public Matrix slicem(final int idx, final boolean dorow) {
204 if (dorow) {
205 if (idx < numRows) {
206 final Matrix retval = new AGDenseMatrix(1,numColumns);
207 for (int j = 0; j < numColumns; j++) {
208 retval.set(0, j, this.get(idx, j));
209 }
210 return retval;
211 } else {
212 throw new IllegalArgumentException("row index " + idx
213 + " greater than " + numRows);
214 }
215 } else {
216 if (idx < numColumns) {
217 final double[][] vals = new double[numRows][1];
218 System.arraycopy(data, idx * numRows, vals, 0, numRows);
219 final Matrix retval = new AGDenseMatrix(vals);
220 return retval;
221 } else {
222 throw new IllegalArgumentException("column index" + idx
223 + "greater than " + numColumns);
224 }
225
226 }
227 }
228
229
230
231 public Vector sum(final int dim) {
232 Vector answ;
233 switch (dim) {
234 case 1:
235 final double[] colsum = new double[this.numColumns()];
236 for (MatrixEntry e : this) {
237 colsum[e.column()] += e.get();
238 }
239 answ = new DenseVector(colsum,false);
240 break;
241 case 2:
242 final double[] rowsum = new double[this.numRows()];
243 for (MatrixEntry e : this) {
244 rowsum[e.row()] += e.get();
245 }
246 answ = new DenseVector(rowsum, false);
247 break;
248 default:
249
250 throw new UnsupportedOperationException(
251 "Matrix.sum() not implemented for higher dimensions");
252
253 }
254
255 return answ;
256 }
257
258
259
260
261
262
263
264 public static Matrix identity(final int ncentres) {
265 final Matrix retval = new AGDenseMatrix(ncentres, ncentres);
266 for (int i = 0; i < ncentres; i++) {
267 retval.set(i, i, 1.0);
268 }
269 return retval;
270 }
271
272 public Matrix pow(final double i) {
273
274 final Matrix retval = new AGDenseMatrix(this.numRows, this.numColumns);
275 for (MatrixEntry e : this) {
276 retval.set(e.row(), e.column(), Math.pow(e.get(), i));
277 }
278 return retval;
279 }
280
281 public Matrix sliceCol(final int colstart, final int ncols) {
282 final AGDenseMatrix retval = new AGDenseMatrix(this.numRows, ncols);
283 System.arraycopy(this.data, colstart*this.numRows, retval.data, 0, ncols*this.numRows);
284
285 return retval;
286 }
287 public Vector sliceCol(final int colstart) {
288
289 final double vals[] = new double[this.numRows];
290 System.arraycopy(this.data, colstart*this.numRows, vals, 0, this.numRows);
291
292 return new DenseVector(vals);
293 }
294
295 public Matrix ones() {
296 Arrays.fill(data, 1.0);
297 return this;
298 }
299
300 public AGDenseMatrix reshape(final int irow, final int icol) {
301 if(this.numColumns !=0 && this.numRows !=0){
302 if(irow*icol != data.length){
303 throw new IllegalArgumentException("new total size of matrix is different");
304 }
305 else {
306 numColumns = icol;
307 numRows = irow;
308
309 }
310 }
311 return this;
312 }
313
314 public DenseVector asVector() {
315 return this.asVector(0,data.length - 1);
316 }
317
318 public DenseVector asVector(int istart, int iend) {
319
320 if( istart < 0){
321 istart = 0;
322 }
323 int ncopy;
324 if( iend <=0 || iend >= data.length){
325 iend = data.length -1;
326 }
327 ncopy = iend - istart + 1;
328 final double[] vecarray = new double[ncopy];
329 System.arraycopy(data, istart, vecarray, 0, ncopy);
330 return new DenseVector(vecarray, false);
331
332
333 }
334 public Vector asVector(int istart) {
335 return this.asVector(istart, data.length -1);
336 }
337
338
339
340
341
342
343
344 public double det() {
345 if(!isSquare()){
346 throw new IllegalArgumentException("matrix must be square for determinant");
347 }
348 double det;
349 final int[] piv = new int[numRows];
350 final double[] indata = new double[data.length];
351 System.arraycopy(data, 0, indata, 0, data.length);
352 final intW info = new intW(0);
353 LAPACK.getInstance().dgetrf(numRows, numColumns, indata, Matrices.ld(numColumns), piv, info);
354
355 if (info.val > 0)
356 throw new MatrixSingularException();
357 else if (info.val < 0)
358 throw new IllegalArgumentException();
359 det = 1;
360 for (int i = 0; i < numColumns; i++) {
361 if(i+1 == piv[i]){
362 det *= indata[i + i*numColumns];
363 } else
364 {
365 det *= -indata[i + i*numColumns];
366 }
367 }
368 return det;
369
370 }
371
372 public double trace() {
373 if (isSquare()) {
374 double retval = 0;
375 for(int i = 0 ; i < numColumns; i++){
376 retval += this.get(i, i);
377
378 }
379 return retval;
380 } else {
381
382
383 throw new UnsupportedOperationException("Matrix.trace() not only for square matrices");
384 }
385 }
386
387 public Matrix inv() {
388 if(!isSquare()){
389 throw new IllegalArgumentException("matrix must be square for inversion");
390 }
391 final int[] piv = new int[numRows];
392 final double[] indata = new double[data.length];
393 System.arraycopy(data, 0, indata, 0, data.length);
394 final intW info = new intW(0);
395 LAPACK.getInstance().dgetrf(numRows, numColumns, indata, Matrices.ld(numColumns), piv, info);
396
397 if (info.val > 0)
398 throw new MatrixSingularException();
399 else if (info.val < 0)
400 throw new IllegalArgumentException();
401
402 int nb = LAPACK.getInstance().ilaenv( 1, "DGETRI", " ", numRows, -1, -1, -1);
403 double[] work = new double[nb*numRows];
404 LAPACK.getInstance().dgetri(numRows, indata, numRows, piv, work, nb*numRows, info);
405 if (info.val > 0)
406 throw new MatrixSingularException();
407 else if (info.val < 0)
408 throw new IllegalArgumentException();
409
410
411 return new AGDenseMatrix(numRows, numColumns, indata);
412
413
414
415
416
417
418
419 }
420
421
422
423
424
425
426
427 public Matrix divide(final Matrix b, final Matrix c) {
428 checkSize(b);
429 checkSize(c);
430 for (MatrixEntry me : c) {
431 final int row = me.row(), col = me.column();
432 me.set(this.get(row, col) / b.get(row, col));
433 }
434 return c;
435 }
436
437
438
439
440
441
442 public Vector sliceRow(final int idx) {
443 return this.slicev(idx, true);
444 }
445
446 public Matrix add(double a) {
447 for (MatrixEntry me : this) {
448 me.set(me.get()+a);
449 }
450 return this;
451 }
452
453
454
455
456
457 public Matrix sliceRowM(int k) {
458 return slicem(k, true);
459 }
460
461
462 public int insert(DenseMatrix m){
463 return insert(m,0);
464 }
465
466
467 public int insert(DenseMatrix m, int position){
468 if(m.numRows * m.numColumns + position > this.numRows*this.numColumns)
469 {
470 throw new IllegalArgumentException("cannot insert m - too big");
471 }
472 else {
473 System.arraycopy(m.data, 0, this.data, position, m.numColumns*m.numRows);
474 return position + m.numRows*m.numColumns;
475 }
476 }
477 public int insert(DenseVector vector, int inpos) {
478 if (vector.size + inpos > this.numColumns*this.numRows) {
479 throw new IllegalArgumentException("cannot insert m - too big");
480 } else {
481 System.arraycopy(vector.getData(), 0, this.data, inpos, vector.size);
482 return inpos + vector.size;
483 }
484
485 }
486
487 public AGDenseMatrix append(Matrix mm) {
488
489 AGDenseMatrix m = (AGDenseMatrix) mm;
490 if(numColumns != 0 && numColumns != m.numColumns){
491 throw new IllegalArgumentException("appended matrix must have same number of columns") ;
492 }
493 else {
494 double arr[] = new double[numColumns*numRows + m.numColumns* m.numRows];
495 System.arraycopy(data, 0, arr, 0, numColumns*numRows);
496 System.arraycopy(m.data, 0, arr, numRows*numColumns, m.numColumns*m.numRows);
497 numRows += m.numRows;
498 data = arr;
499 }
500 return this;
501 }
502
503 public double asScalar() {
504
505 if(numColumns == 1 && numRows == 1)
506 {
507 return data[0];
508 }
509 else {
510 throw new IllegalStateException("cannot take scalar unless matix is 1x1");
511 }
512
513 }
514
515 public AGDenseMatrix append(Vector mm) {
516 DenseVector m = (DenseVector) mm;
517 if(!(numColumns == 1 || numColumns == 0) ){
518 throw new IllegalArgumentException("matrix must have 1 (or 0) column") ;
519 } else {
520 double arr[] = new double[numRows + m.size()];
521 System.arraycopy(data, 0, arr, 0, numRows);
522 System.arraycopy(m.getData(), 0, arr, numRows*numColumns, m.size());
523 numRows += m.size();
524 if(numColumns == 0) numColumns = 1;
525 data = arr;
526
527 }
528 return this;
529 }
530 public void appendCol(AGDenseMatrix r) {
531 if (r.numRows() != numRows) {
532 throw new IllegalArgumentException("matrices must have same number of rows") ;
533 } else {
534 double arr[] = new double[numRows*(numColumns+r.numColumns())];
535 System.arraycopy(data, 0, arr, 0, numColumns* numRows);
536 System.arraycopy(r.data, 0, arr, numRows*numColumns, r.numColumns()*r.numRows());
537 numColumns += r.numColumns();
538 data = arr;
539 }
540
541
542 }
543
544
545 public Matrix setColumn(int idx, Vector v) {
546 this.setColumn(idx, new AGDenseMatrix(v));
547 return this;
548 }
549
550 public Matrix setColumn(int idx, Matrix v) {
551 if(v.numColumns() != 1) {
552 throw new IllegalArgumentException("setting matrix must have 1 column");
553 } else {
554 if(v.numRows() != numRows){
555 throw new IllegalArgumentException("setting matrix must have same number of rows");
556 }
557 else {
558 for (int i = 0; i < v.numRows(); i++) {
559 this.set(i, idx, v.get(i, 0));
560 }
561 }
562 }
563 return this;
564 }
565
566 public AGDenseMatrix append(double d) {
567 if(!(numColumns == 1 || numColumns == 0)){
568 throw new IllegalArgumentException("matrix must have 1 column") ;
569 } else {
570 double arr[] = new double[numRows + 1];
571 System.arraycopy(data, 0, arr, 0, numRows);
572 arr[numRows] = d;
573 numRows += 1;
574 numColumns = 1;
575 data = arr;
576
577 }
578 return this;
579 }
580
581 public Matrix setRow(int k, Matrix m) {
582
583 if (m.numRows() == 1 && m.numColumns() == numColumns) {
584 for (int i = 0; i < m.numColumns(); i++) {
585 this.set(k, i, m.get(0, i));
586
587 }
588 } else {
589 throw new IllegalArgumentException("The setting matrix is not of the correct dimensions");
590 }
591 return this;
592 }
593
594 public Matrix setRow(int k, Vector v) {
595 if (v.size() == numColumns) {
596 for (int i = 0; i < v.size(); i++) {
597 this.set(k, i, v.get(i));
598 }
599 } else {
600 throw new IllegalArgumentException("The setting matrix is not of the correct dimensions");
601 }
602 return this;
603 }
604
605
606 public AGDenseMatrix selectCols(int cols[], int ncols){
607 AGDenseMatrix retval = null;
608 if (ncols <= cols.length && ncols <= numColumns) {
609 retval = new AGDenseMatrix(this.numRows, ncols);
610 for (int i = 0; i < ncols; i++) {
611 int icol = cols[i];
612 System.arraycopy(this.data, icol*this.numRows, retval.data, i*this.numRows, this.numRows);
613 }
614 } else {
615 throw new IllegalArgumentException("To many columns specified");
616 }
617 return retval;
618
619
620 }
621
622 public Matrix delCol(int k) {
623 if(k >= numColumns || k < 0){
624 throw new IllegalArgumentException("illegal column number");
625
626 }
627 if (k < numColumns -1){
628 System.arraycopy(data, (k+1)*numRows, data, k*numRows, (numColumns-k-1)*numRows);
629 }
630 numColumns--;
631 return this;
632
633 }
634
635 public Matrix delRow(int del) {
636
637 double[] newarr = new double[(numRows -1)*
638 numColumns];
639 int isrc, idest, i = 0, niter;
640 if(del == 0){
641 isrc = 1;
642 idest = 0;
643 niter = numColumns;
644 } else
645 {
646 System.arraycopy(data, 0, newarr, 0, del);
647 isrc = del +1;
648 idest = del;
649 niter = numColumns -1;
650
651 }
652 while ( i < niter)
653 {
654 System.arraycopy(data, isrc, newarr, idest, numRows -1);
655 isrc += numRows;
656 idest += numRows -1;
657 ++i;
658 }
659
660 if(del != 0){
661
662 System.arraycopy(data, isrc, newarr, idest, numRows -1 - del);
663 }
664 data = newarr;
665 numRows--;
666 return this;
667 }
668
669
670
671 }
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696