001: /* AUTO-GENERATED */
002: package JSci.maths.matrices;
003:
004: import JSci.maths.ExtraMath;
005: import JSci.maths.Mapping;
006: import JSci.maths.DimensionException;
007: import JSci.maths.MaximumIterationsExceededException;
008: import JSci.maths.vectors.AbstractDoubleVector;
009: import JSci.maths.vectors.DoubleVector;
010: import JSci.maths.groups.AbelianGroup;
011: import JSci.maths.algebras.*;
012: import JSci.maths.fields.*;
013:
014: /**
015: * The DoubleDiagonalMatrix class provides an object for encapsulating double diagonal matrices.
016: * @version 2.2
017: * @author Mark Hale
018: */
019: public class DoubleDiagonalMatrix extends AbstractDoubleSquareMatrix
020: implements DiagonalMatrix {
021: /**
022: * Diagonal data.
023: */
024: protected final double diag[];
025:
026: /**
027: * Constructs an empty matrix.
028: * @param size the number of rows/columns
029: */
030: public DoubleDiagonalMatrix(final int size) {
031: this (new double[size]);
032: }
033:
034: /**
035: * Constructs a matrix from an array.
036: * Any non-diagonal elements in the array are ignored.
037: * @param array an assigned value
038: * @exception MatrixDimensionException If the array is not square.
039: */
040: public DoubleDiagonalMatrix(final double array[][]) {
041: this (array.length);
042: for (int i = 0; i < array.length; i++) {
043: if (array[i].length != array.length)
044: throw new MatrixDimensionException(
045: "Array is not square.");
046: diag[i] = array[i][i];
047: }
048: }
049:
050: /**
051: * Constructs a matrix by wrapping an array containing the diagonal elements.
052: * @param array an assigned value
053: */
054: public DoubleDiagonalMatrix(final double array[]) {
055: super (array.length);
056: diag = array;
057: }
058:
059: /**
060: * Creates an identity matrix.
061: * @param size the number of rows/columns
062: */
063: public static DoubleDiagonalMatrix identity(final int size) {
064: double array[] = new double[size];
065: for (int i = 0; i < size; i++)
066: array[i] = 1;
067: return new DoubleDiagonalMatrix(array);
068: }
069:
070: /**
071: * Compares two ${nativeTyp} matrices for equality.
072: * @param m a double matrix
073: */
074: public boolean equals(AbstractDoubleMatrix m, double tol) {
075: if (m instanceof DiagonalMatrix) {
076: if (numRows != m.rows() || numCols != m.columns())
077: return false;
078: double sumSqr = 0;
079: double delta = diag[0] - m.getElement(0, 0);
080: sumSqr += delta * delta;
081: for (int i = 1; i < numRows; i++) {
082: delta = diag[i] - m.getElement(i, i);
083: sumSqr += delta * delta;
084: }
085: return (sumSqr <= tol * tol);
086: } else {
087: return false;
088: }
089: }
090:
091: /**
092: * Returns a string representing this matrix.
093: */
094: public String toString() {
095: final StringBuffer buf = new StringBuffer(5 * numRows * numCols);
096: for (int i = 0; i < numRows; i++) {
097: for (int j = 0; j < numCols; j++) {
098: buf.append(getElement(i, j));
099: buf.append(' ');
100: }
101: buf.append('\n');
102: }
103: return buf.toString();
104: }
105:
106: /**
107: * Converts this matrix to an integer matrix.
108: * @return an integer matrix
109: */
110: public AbstractIntegerMatrix toIntegerMatrix() {
111: final int array[] = new int[numRows];
112: for (int i = 0; i < numRows; i++)
113: array[i] = Math.round((float) diag[i]);
114: return new IntegerDiagonalMatrix(array);
115: }
116:
117: /**
118: * Converts this matrix to a complex matrix.
119: * @return a complex matrix
120: */
121: public AbstractComplexMatrix toComplexMatrix() {
122: final double array[] = new double[numRows];
123: for (int i = 0; i < numRows; i++)
124: array[i] = diag[i];
125: return new ComplexDiagonalMatrix(array, new double[numRows]);
126: }
127:
128: /**
129: * Returns an element of the matrix.
130: * @param i row index of the element
131: * @param j column index of the element
132: * @exception MatrixDimensionException If attempting to access an invalid element.
133: */
134: public double getElement(int i, int j) {
135: if (i >= 0 && i < numRows && j >= 0 && j < numCols) {
136: if (i == j)
137: return diag[i];
138: else
139: return 0;
140: } else
141: throw new MatrixDimensionException(getInvalidElementMsg(i,
142: j));
143: }
144:
145: /**
146: * Sets the value of an element of the matrix.
147: * Should only be used to initialise this matrix.
148: * @param i row index of the element
149: * @param j column index of the element
150: * @param x a number
151: * @exception MatrixDimensionException If attempting to access an invalid element.
152: */
153: public void setElement(int i, int j, final double x) {
154: if (i >= 0 && i < numRows && j >= 0 && j < numCols) {
155: if (i == j)
156: diag[i] = x;
157: else
158: throw new MatrixDimensionException(
159: getInvalidElementMsg(i, j));
160: } else
161: throw new MatrixDimensionException(getInvalidElementMsg(i,
162: j));
163: }
164:
165: /**
166: * Returns true if this matrix is symmetric.
167: */
168: public boolean isSymmetric() {
169: return true;
170: }
171:
172: /**
173: * Returns the determinant.
174: */
175: public double det() {
176: double det = diag[0];
177: for (int i = 1; i < numRows; i++)
178: det *= diag[i];
179: return det;
180: }
181:
182: /**
183: * Returns the trace.
184: */
185: public double trace() {
186: double tr = diag[0];
187: for (int i = 1; i < numRows; i++)
188: tr += diag[i];
189: return tr;
190: }
191:
192: /**
193: * Returns the l<sup><img border=0 alt="infinity" src="doc-files/infinity.gif"></sup>-norm.
194: * @author Taber Smith
195: */
196: public double infNorm() {
197: double result = Math.abs(diag[0]);
198: double tmpResult;
199: for (int i = 1; i < numRows; i++) {
200: tmpResult = Math.abs(diag[i]);
201: if (tmpResult > result)
202: result = tmpResult;
203: }
204: return result;
205: }
206:
207: /**
208: * Returns the Frobenius (l<sup>2</sup>) norm.
209: * @author Taber Smith
210: */
211: public double frobeniusNorm() {
212: double result = diag[0];
213: for (int i = 1; i < numRows; i++)
214: result = ExtraMath.hypot(result, diag[i]);
215: return result;
216: }
217:
218: /**
219: * Returns the operator norm.
220: * @exception MaximumIterationsExceededException If it takes more than 50 iterations to determine an eigenvalue.
221: */
222: public double operatorNorm()
223: throws MaximumIterationsExceededException {
224: return infNorm();
225: }
226:
227: //============
228: // OPERATIONS
229: //============
230:
231: // ADDITION
232:
233: /**
234: * Returns the addition of this matrix and another.
235: * @param m a double matrix
236: * @exception MatrixDimensionException If the matrices are different sizes.
237: */
238: public AbstractDoubleSquareMatrix add(
239: final AbstractDoubleSquareMatrix m) {
240: if (m instanceof DoubleDiagonalMatrix)
241: return add((DoubleDiagonalMatrix) m);
242: if (m instanceof DiagonalMatrix)
243: return addDiagonal(m);
244: if (m instanceof DoubleTridiagonalMatrix)
245: return add((DoubleTridiagonalMatrix) m);
246: if (m instanceof TridiagonalMatrix)
247: return addTridiagonal(m);
248: if (m instanceof DoubleSquareMatrix)
249: return add((DoubleSquareMatrix) m);
250:
251: if (numRows == m.rows() && numCols == m.columns()) {
252: final double array[][] = new double[numRows][numCols];
253: for (int i = 0; i < numRows; i++) {
254: array[i][0] = m.getElement(i, 0);
255: for (int j = 1; j < numCols; j++)
256: array[i][j] = m.getElement(i, j);
257: }
258: for (int i = 0; i < numRows; i++)
259: array[i][i] += diag[i];
260: return new DoubleSquareMatrix(array);
261: } else {
262: throw new MatrixDimensionException(
263: "Matrices are different sizes.");
264: }
265: }
266:
267: public DoubleSquareMatrix add(final DoubleSquareMatrix m) {
268: if (numRows == m.numRows && numCols == m.numCols) {
269: final double array[][] = new double[numRows][numCols];
270: for (int i = 0; i < numRows; i++)
271: System.arraycopy(m.matrix[i], 0, array[i], 0, numRows);
272: for (int i = 0; i < numRows; i++)
273: array[i][i] += diag[i];
274: return new DoubleSquareMatrix(array);
275: } else
276: throw new MatrixDimensionException(
277: "Matrices are different sizes.");
278: }
279:
280: /**
281: * Returns the addition of this matrix and another.
282: * @param m a double tridiagonal matrix
283: * @exception MatrixDimensionException If the matrices are different sizes.
284: */
285: public DoubleTridiagonalMatrix add(final DoubleTridiagonalMatrix m) {
286: if (numRows == m.numRows) {
287: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
288: numRows);
289: System.arraycopy(m.ldiag, 0, ans.ldiag, 0, m.ldiag.length);
290: System.arraycopy(m.udiag, 0, ans.udiag, 0, m.udiag.length);
291: ans.diag[0] = diag[0] + m.diag[0];
292: for (int i = 1; i < numRows; i++)
293: ans.diag[i] = diag[i] + m.diag[i];
294: return ans;
295: } else
296: throw new MatrixDimensionException(
297: "Matrices are different sizes.");
298: }
299:
300: private DoubleTridiagonalMatrix addTridiagonal(
301: final AbstractDoubleSquareMatrix m) {
302: int mRow = numRows;
303: if (mRow == m.rows()) {
304: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
305: mRow);
306: ans.diag[0] = diag[0] + m.getElement(0, 0);
307: ans.udiag[0] = m.getElement(0, 1);
308: mRow--;
309: for (int i = 1; i < mRow; i++) {
310: ans.ldiag[i] = m.getElement(i, i - 1);
311: ans.diag[i] = diag[i] + m.getElement(i, i);
312: ans.udiag[i] = m.getElement(i, i + 1);
313: }
314: ans.ldiag[mRow] = m.getElement(mRow, mRow - 1);
315: ans.diag[mRow] = diag[mRow] + m.getElement(mRow, mRow);
316: return ans;
317: } else {
318: throw new MatrixDimensionException(
319: "Matrices are different sizes.");
320: }
321: }
322:
323: /**
324: * Returns the addition of this matrix and another.
325: * @param m a double diagonal matrix
326: * @exception MatrixDimensionException If the matrices are different sizes.
327: */
328: public DoubleDiagonalMatrix add(final DoubleDiagonalMatrix m) {
329: if (numRows == m.numRows) {
330: final double array[] = new double[numRows];
331: array[0] = diag[0] + m.diag[0];
332: for (int i = 1; i < numRows; i++)
333: array[i] = diag[i] + m.diag[i];
334: return new DoubleDiagonalMatrix(array);
335: } else
336: throw new MatrixDimensionException(
337: "Matrices are different sizes.");
338: }
339:
340: private DoubleDiagonalMatrix addDiagonal(
341: final AbstractDoubleSquareMatrix m) {
342: if (numRows == m.numRows) {
343: final double array[] = new double[numRows];
344: array[0] = diag[0] + m.getElement(0, 0);
345: for (int i = 1; i < numRows; i++)
346: array[i] = diag[i] + m.getElement(i, i);
347: return new DoubleDiagonalMatrix(array);
348: } else
349: throw new MatrixDimensionException(
350: "Matrices are different sizes.");
351: }
352:
353: // SUBTRACTION
354:
355: /**
356: * Returns the subtraction of this matrix by another.
357: * @param m a double matrix
358: * @exception MatrixDimensionException If the matrices are different sizes.
359: */
360: public AbstractDoubleSquareMatrix subtract(
361: final AbstractDoubleSquareMatrix m) {
362: if (m instanceof DoubleDiagonalMatrix)
363: return subtract((DoubleDiagonalMatrix) m);
364: if (m instanceof DiagonalMatrix)
365: return subtractDiagonal(m);
366: if (m instanceof DoubleTridiagonalMatrix)
367: return subtract((DoubleTridiagonalMatrix) m);
368: if (m instanceof TridiagonalMatrix)
369: return subtractTridiagonal(m);
370: if (m instanceof DoubleSquareMatrix)
371: return subtract((DoubleSquareMatrix) m);
372:
373: if (numRows == m.rows() && numCols == m.columns()) {
374: final double array[][] = new double[numRows][numCols];
375: for (int i = 0; i < numRows; i++) {
376: array[i][0] = -m.getElement(i, 0);
377: for (int j = 1; j < numCols; j++)
378: array[i][j] = -m.getElement(i, j);
379: }
380: for (int i = 0; i < numRows; i++)
381: array[i][i] += diag[i];
382: return new DoubleSquareMatrix(array);
383: } else {
384: throw new MatrixDimensionException(
385: "Matrices are different sizes.");
386: }
387: }
388:
389: public DoubleSquareMatrix subtract(final DoubleSquareMatrix m) {
390: if (numRows == m.numRows && numCols == m.numCols) {
391: final double array[][] = new double[numRows][numCols];
392: for (int i = 0; i < numRows; i++) {
393: array[i][0] = -m.matrix[i][0];
394: for (int j = 1; j < numCols; j++)
395: array[i][j] = -m.matrix[i][j];
396: }
397: for (int i = 0; i < numRows; i++)
398: array[i][i] += diag[i];
399: return new DoubleSquareMatrix(array);
400: } else
401: throw new MatrixDimensionException(
402: "Matrices are different sizes.");
403: }
404:
405: /**
406: * Returns the subtraction of this matrix and another.
407: * @param m a double tridiagonal matrix
408: * @exception MatrixDimensionException If the matrices are different sizes.
409: */
410: public DoubleTridiagonalMatrix subtract(
411: final DoubleTridiagonalMatrix m) {
412: int mRow = numRows;
413: if (mRow == m.numRows) {
414: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
415: mRow);
416: ans.diag[0] = diag[0] - m.diag[0];
417: ans.udiag[0] = -m.udiag[0];
418: mRow--;
419: for (int i = 1; i < mRow; i++) {
420: ans.ldiag[i] = -m.ldiag[i];
421: ans.diag[i] = diag[i] - m.diag[i];
422: ans.udiag[i] = -m.udiag[i];
423: }
424: ans.ldiag[mRow] = -m.ldiag[mRow];
425: ans.diag[mRow] = diag[mRow] - m.diag[mRow];
426: return ans;
427: } else
428: throw new MatrixDimensionException(
429: "Matrices are different sizes.");
430: }
431:
432: private DoubleTridiagonalMatrix subtractTridiagonal(
433: final AbstractDoubleSquareMatrix m) {
434: int mRow = numRows;
435: if (mRow == m.rows()) {
436: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
437: mRow);
438: ans.diag[0] = diag[0] - m.getElement(0, 0);
439: ans.udiag[0] = -m.getElement(0, 1);
440: mRow--;
441: for (int i = 1; i < mRow; i++) {
442: ans.ldiag[i] = -m.getElement(i, i - 1);
443: ans.diag[i] = diag[i] - m.getElement(i, i);
444: ans.udiag[i] = -m.getElement(i, i + 1);
445: }
446: ans.ldiag[mRow] = -m.getElement(mRow, mRow - 1);
447: ans.diag[mRow] = diag[mRow] - m.getElement(mRow, mRow);
448: return ans;
449: } else {
450: throw new MatrixDimensionException(
451: "Matrices are different sizes.");
452: }
453: }
454:
455: /**
456: * Returns the subtraction of this matrix and another.
457: * @param m a double diagonal matrix
458: * @exception MatrixDimensionException If the matrices are different sizes.
459: */
460: public DoubleDiagonalMatrix subtract(final DoubleDiagonalMatrix m) {
461: if (numRows == m.numRows) {
462: final double array[] = new double[numRows];
463: array[0] = diag[0] - m.diag[0];
464: for (int i = 1; i < numRows; i++)
465: array[i] = diag[i] - m.diag[i];
466: return new DoubleDiagonalMatrix(array);
467: } else
468: throw new MatrixDimensionException(
469: "Matrices are different sizes.");
470: }
471:
472: private DoubleDiagonalMatrix subtractDiagonal(
473: final AbstractDoubleSquareMatrix m) {
474: if (numRows == m.numRows) {
475: final double array[] = new double[numRows];
476: array[0] = diag[0] - m.getElement(0, 0);
477: for (int i = 1; i < numRows; i++)
478: array[i] = diag[i] - m.getElement(i, i);
479: return new DoubleDiagonalMatrix(array);
480: } else
481: throw new MatrixDimensionException(
482: "Matrices are different sizes.");
483: }
484:
485: // SCALAR MULTIPLICATION
486:
487: /**
488: * Returns the multiplication of this matrix by a scalar.
489: * @param x a double.
490: * @return a double diagonal matrix.
491: */
492: public AbstractDoubleMatrix scalarMultiply(final double x) {
493: final double array[] = new double[numRows];
494: array[0] = x * diag[0];
495: for (int i = 1; i < numRows; i++)
496: array[i] = x * diag[i];
497: return new DoubleDiagonalMatrix(array);
498: }
499:
500: // SCALAR DIVISON
501:
502: /**
503: * Returns the division of this matrix by a scalar.
504: * @param x a double.
505: * @return a double diagonal matrix.
506: */
507: public AbstractDoubleMatrix scalarDivide(final double x) {
508: final double array[] = new double[numRows];
509: array[0] = diag[0] / x;
510: for (int i = 1; i < numRows; i++)
511: array[i] = diag[i] / x;
512: return new DoubleDiagonalMatrix(array);
513: }
514:
515: // SCALAR PRODUCT
516:
517: /**
518: * Returns the scalar product of this matrix and another.
519: * @param m a double matrix.
520: * @exception MatrixDimensionException If the matrices are different sizes.
521: */
522: public double scalarProduct(final AbstractDoubleSquareMatrix m) {
523: if (m instanceof DoubleDiagonalMatrix)
524: return scalarProduct((DoubleDiagonalMatrix) m);
525: if (m instanceof DoubleTridiagonalMatrix)
526: return scalarProduct((DoubleTridiagonalMatrix) m);
527: if (m instanceof DoubleSquareMatrix)
528: return scalarProduct((DoubleSquareMatrix) m);
529:
530: if (numRows == m.rows() && numCols == m.columns()) {
531: double ans = diag[0] * m.getElement(0, 0);
532: for (int i = 1; i < numRows; i++)
533: ans += diag[i] * m.getElement(i, i);
534: return ans;
535: } else {
536: throw new MatrixDimensionException(
537: "Matrices are different sizes.");
538: }
539: }
540:
541: public double scalarProduct(final DoubleSquareMatrix m) {
542: if (numRows == m.numRows && numCols == m.numCols) {
543: double ans = diag[0] * m.matrix[0][0];
544: for (int i = 1; i < numRows; i++)
545: ans += diag[i] * m.matrix[i][i];
546: return ans;
547: } else
548: throw new MatrixDimensionException(
549: "Matrices are different sizes.");
550: }
551:
552: public double scalarProduct(final DoubleTridiagonalMatrix m) {
553: if (numRows == m.numRows) {
554: double ans = diag[0] * m.diag[0];
555: for (int i = 1; i < numRows; i++)
556: ans += diag[i] * m.diag[i];
557: return ans;
558: } else
559: throw new MatrixDimensionException(
560: "Matrices are different sizes.");
561: }
562:
563: public double scalarProduct(final DoubleDiagonalMatrix m) {
564: if (numRows == m.numRows) {
565: double ans = diag[0] * m.diag[0];
566: for (int i = 1; i < numRows; i++)
567: ans += diag[i] * m.diag[i];
568: return ans;
569: } else
570: throw new MatrixDimensionException(
571: "Matrices are different sizes.");
572: }
573:
574: // MATRIX MULTIPLICATION
575:
576: /**
577: * Returns the multiplication of a vector by this matrix.
578: * @param v a double vector.
579: * @exception DimensionException If the matrix and vector are incompatible.
580: */
581: public AbstractDoubleVector multiply(final AbstractDoubleVector v) {
582: if (numCols == v.dimension()) {
583: final double array[] = new double[numRows];
584: array[0] = diag[0] * v.getComponent(0);
585: for (int i = 1; i < numRows; i++)
586: array[i] = diag[i] * v.getComponent(i);
587: return new DoubleVector(array);
588: } else {
589: throw new DimensionException(
590: "Matrix and vector are incompatible.");
591: }
592: }
593:
594: /**
595: * Returns the multiplication of this matrix and another.
596: * @param m a double matrix
597: * @return a AbstractDoubleMatrix or a AbstractDoubleSquareMatrix as appropriate
598: * @exception MatrixDimensionException If the matrices are incompatible.
599: */
600: public AbstractDoubleSquareMatrix multiply(
601: final AbstractDoubleSquareMatrix m) {
602: if (m instanceof DoubleDiagonalMatrix)
603: return multiply((DoubleDiagonalMatrix) m);
604: if (m instanceof DiagonalMatrix)
605: return multiplyDiagonal(m);
606: if (m instanceof DoubleTridiagonalMatrix)
607: return multiply((DoubleTridiagonalMatrix) m);
608: if (m instanceof TridiagonalMatrix)
609: return multiplyTridiagonal(m);
610: if (m instanceof DoubleSquareMatrix)
611: return multiply((DoubleSquareMatrix) m);
612:
613: if (numCols == m.rows()) {
614: final int mColumns = m.columns();
615: final double array[][] = new double[numRows][mColumns];
616: for (int i = 0; i < numRows; i++) {
617: array[i][0] = diag[0] * m.getElement(i, 0);
618: for (int j = 1; j < mColumns; j++)
619: array[i][j] = diag[i] * m.getElement(i, j);
620: }
621: return new DoubleSquareMatrix(array);
622: } else {
623: throw new MatrixDimensionException("Incompatible matrices.");
624: }
625: }
626:
627: public DoubleSquareMatrix multiply(final DoubleSquareMatrix m) {
628: if (numCols == m.numRows) {
629: final double array[][] = new double[numRows][m.numCols];
630: for (int i = 0; i < numRows; i++) {
631: array[i][0] = diag[0] * m.matrix[i][0];
632: for (int j = 1; j < m.numCols; j++)
633: array[i][j] = diag[i] * m.matrix[i][j];
634: }
635: return new DoubleSquareMatrix(array);
636: } else
637: throw new MatrixDimensionException("Incompatible matrices.");
638: }
639:
640: public DoubleTridiagonalMatrix multiply(
641: final DoubleTridiagonalMatrix m) {
642: int mRow = numRows;
643: if (numCols == m.numRows) {
644: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
645: mRow);
646: ans.diag[0] = diag[0] * m.diag[0];
647: ans.udiag[0] = diag[0] * m.udiag[0];
648: mRow--;
649: for (int i = 1; i < mRow; i++) {
650: ans.ldiag[i] = diag[i] * m.ldiag[i];
651: ans.diag[i] = diag[i] * m.diag[i];
652: ans.udiag[i] = diag[i] * m.udiag[i];
653: }
654: ans.ldiag[mRow] = diag[mRow] * m.ldiag[mRow];
655: ans.diag[mRow] = diag[mRow] * m.diag[mRow];
656: return ans;
657: } else
658: throw new MatrixDimensionException("Incompatible matrices.");
659: }
660:
661: private DoubleTridiagonalMatrix multiplyTridiagonal(
662: final AbstractDoubleSquareMatrix m) {
663: int mRow = numRows;
664: if (numCols == m.rows()) {
665: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
666: mRow);
667: ans.diag[0] = diag[0] * m.getElement(0, 0);
668: ans.udiag[0] = diag[0] * m.getElement(0, 1);
669: mRow--;
670: for (int i = 1; i < mRow; i++) {
671: ans.ldiag[i] = diag[i] * m.getElement(i, i - 1);
672: ans.diag[i] = diag[i] * m.getElement(i, i);
673: ans.udiag[i] = diag[i] * m.getElement(i, i + 1);
674: }
675: ans.ldiag[mRow] = diag[mRow] * m.getElement(mRow, mRow - 1);
676: ans.diag[mRow] = diag[mRow] * m.getElement(mRow, mRow);
677: return ans;
678: } else {
679: throw new MatrixDimensionException("Incompatible matrices.");
680: }
681: }
682:
683: public DoubleDiagonalMatrix multiply(final DoubleDiagonalMatrix m) {
684: if (numCols == m.numRows) {
685: final double array[] = new double[numRows];
686: array[0] = diag[0] * m.diag[0];
687: for (int i = 1; i < numRows; i++) {
688: array[i] = diag[i] * m.diag[i];
689: }
690: return new DoubleDiagonalMatrix(array);
691: } else
692: throw new MatrixDimensionException("Incompatible matrices.");
693: }
694:
695: private DoubleDiagonalMatrix multiplyDiagonal(
696: final AbstractDoubleSquareMatrix m) {
697: if (numCols == m.rows()) {
698: final double array[] = new double[numRows];
699: array[0] = diag[0] * m.getElement(0, 0);
700: for (int i = 1; i < numRows; i++) {
701: array[i] = diag[i] * m.getElement(i, i);
702: }
703: return new DoubleDiagonalMatrix(array);
704: } else {
705: throw new MatrixDimensionException("Incompatible matrices.");
706: }
707: }
708:
709: // TRANSPOSE
710:
711: /**
712: * Returns the transpose of this matrix.
713: * @return a double matrix
714: */
715: public Matrix transpose() {
716: return this ;
717: }
718:
719: // INVERSE
720:
721: /**
722: * Returns the inverse of this matrix.
723: * @return a double diagonal matrix
724: */
725: public AbstractDoubleSquareMatrix inverse() {
726: final double array[] = new double[numRows];
727: array[0] = 1.0 / diag[0];
728: for (int i = 1; i < numRows; i++)
729: array[i] = 1.0 / diag[i];
730: return new DoubleDiagonalMatrix(array);
731: }
732:
733: // LU DECOMPOSITION
734:
735: /**
736: * Returns the LU decomposition of this matrix.
737: * @param pivot an empty array of length <code>rows()+1</code>
738: * to hold the pivot information (null if not interested).
739: * The last array element will contain the parity.
740: * @return an array with [0] containing the L-matrix
741: * and [1] containing the U-matrix.
742: */
743: public AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
744: if (LU != null) {
745: if (pivot != null)
746: System.arraycopy(LUpivot, 0, pivot, 0, pivot.length);
747: return LU;
748: }
749: if (pivot == null)
750: pivot = new int[numRows + 1];
751: for (int i = 0; i < numRows; i++)
752: pivot[i] = i;
753: pivot[numRows] = 1;
754: LU = new AbstractDoubleSquareMatrix[2];
755: LU[0] = DoubleDiagonalMatrix.identity(numRows);
756: LU[1] = this ;
757: LUpivot = new int[pivot.length];
758: System.arraycopy(pivot, 0, LUpivot, 0, pivot.length);
759: return LU;
760: }
761:
762: /**
763: * Returns the LU decomposition of this matrix.
764: * @return an array with [0] containing the L-matrix
765: * and [1] containing the U-matrix.
766: * @jsci.planetmath LUDecomposition
767: */
768: public AbstractDoubleSquareMatrix[] luDecompose() {
769: return luDecompose(null);
770: }
771:
772: // CHOLESKY DECOMPOSITION
773:
774: /**
775: * Returns the Cholesky decomposition of this matrix.
776: * Matrix must be symmetric and positive definite.
777: * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
778: */
779: public AbstractDoubleSquareMatrix[] choleskyDecompose() {
780: final AbstractDoubleSquareMatrix lu[] = new AbstractDoubleSquareMatrix[2];
781: final double array[] = new double[numRows];
782: array[0] = Math.sqrt(diag[0]);
783: for (int i = 1; i < numRows; i++)
784: array[i] = Math.sqrt(diag[i]);
785: lu[0] = new DoubleDiagonalMatrix(array);
786: lu[1] = lu[0];
787: return lu;
788: }
789:
790: // QR DECOMPOSITION
791:
792: /**
793: * Returns the QR decomposition of this matrix.
794: * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
795: * @jsci.planetmath QRDecomposition
796: */
797: public AbstractDoubleSquareMatrix[] qrDecompose() {
798: final AbstractDoubleSquareMatrix qr[] = new AbstractDoubleSquareMatrix[2];
799: qr[0] = DoubleDiagonalMatrix.identity(numRows);
800: qr[1] = this ;
801: return qr;
802: }
803:
804: // SINGULAR VALUE DECOMPOSITION
805:
806: /**
807: * Returns the singular value decomposition of this matrix.
808: * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
809: */
810: public AbstractDoubleSquareMatrix[] singularValueDecompose() {
811: final int N = numRows;
812: final int Nm1 = N - 1;
813: final double arrayU[] = new double[N];
814: final double arrayS[] = new double[N];
815: final double arrayV[] = new double[N];
816: for (int i = 0; i < Nm1; i++) {
817: arrayU[i] = -1.0;
818: arrayS[i] = Math.abs(diag[i]);
819: arrayV[i] = diag[i] < 0.0 ? 1.0 : -1.0;
820: }
821: arrayU[Nm1] = 1.0;
822: arrayS[Nm1] = Math.abs(diag[Nm1]);
823: arrayV[Nm1] = diag[Nm1] < 0.0 ? -1.0 : 1.0;
824: final AbstractDoubleSquareMatrix svd[] = new AbstractDoubleSquareMatrix[3];
825: svd[0] = new DoubleDiagonalMatrix(arrayU);
826: svd[1] = new DoubleDiagonalMatrix(arrayS);
827: svd[2] = new DoubleDiagonalMatrix(arrayV);
828: return svd;
829: }
830:
831: // MAP ELEMENTS
832:
833: /**
834: * Applies a function on all the matrix elements.
835: * @param f a user-defined function
836: * @return a double matrix
837: */
838: public AbstractDoubleMatrix mapElements(final Mapping f) {
839: double zeroValue = f.map(0.0);
840: if (Math.abs(zeroValue) <= JSci.GlobalSettings.ZERO_TOL)
841: return diagonalMap(f);
842: else
843: return generalMap(f, zeroValue);
844: }
845:
846: private AbstractDoubleMatrix diagonalMap(Mapping f) {
847: final double array[] = new double[numRows];
848: array[0] = f.map(diag[0]);
849: for (int i = 1; i < numRows; i++)
850: array[i] = f.map(diag[i]);
851: return new DoubleDiagonalMatrix(array);
852: }
853:
854: private AbstractDoubleMatrix generalMap(Mapping f, double zeroValue) {
855: final double array[][] = new double[numRows][numRows];
856: for (int i = 0; i < numRows; i++) {
857: for (int j = 0; j < numRows; j++) {
858: array[i][j] = zeroValue;
859: }
860: }
861: array[0][0] = f.map(diag[0]);
862: for (int i = 1; i < numRows; i++)
863: array[i][i] = f.map(diag[i]);
864: return new DoubleSquareMatrix(array);
865: }
866: }
|