001: package JSci.maths.matrices;
002:
003: import JSci.GlobalSettings;
004: import JSci.maths.Mapping;
005: import JSci.maths.DimensionException;
006: import JSci.maths.vectors.AbstractDoubleVector;
007: import JSci.maths.vectors.DoubleVector;
008:
009: /**
010: * The DoubleSparseSquareMatrix class provides an object for encapsulating sparse square matrices.
011: * Uses compressed row storage.
012: * @version 1.2
013: * @author Mark Hale
014: */
015: public final class DoubleSparseSquareMatrix extends
016: AbstractDoubleSquareMatrix {
017: /**
018: * Matrix elements.
019: */
020: private double elements[];
021: /**
022: * Sparse indexing data.
023: * Contains the column positions of each element,
024: * e.g. <code>colPos[n]</code> is the column position
025: * of the <code>n</code>th element.
026: */
027: private int colPos[];
028: /**
029: * Sparse indexing data.
030: * Contains the indices of the start of each row,
031: * e.g. <code>rows[i]</code> is the index
032: * where the <code>i</code>th row starts.
033: */
034: private int rows[];
035:
036: /**
037: * Constructs an empty matrix.
038: * @param size the number of rows/columns
039: */
040: public DoubleSparseSquareMatrix(final int size) {
041: super (size);
042: elements = new double[0];
043: colPos = new int[0];
044: rows = new int[numRows + 1];
045: }
046:
047: /**
048: * Constructs a matrix from an array.
049: * @param array an assigned value
050: * @exception MatrixDimensionException If the array is not square.
051: */
052: public DoubleSparseSquareMatrix(final double array[][]) {
053: super (array.length);
054: rows = new int[numRows + 1];
055: int n = 0;
056: for (int i = 0; i < numRows; i++) {
057: if (array[i].length != array.length)
058: throw new MatrixDimensionException(
059: "Array is not square.");
060: for (int j = 0; j < numCols; j++) {
061: if (Math.abs(array[i][j]) > GlobalSettings.ZERO_TOL)
062: n++;
063: }
064: }
065: elements = new double[n];
066: colPos = new int[n];
067: n = 0;
068: for (int i = 0; i < numRows; i++) {
069: rows[i] = n;
070: for (int j = 0; j < numCols; j++) {
071: if (Math.abs(array[i][j]) > GlobalSettings.ZERO_TOL) {
072: elements[n] = array[i][j];
073: colPos[n] = j;
074: n++;
075: }
076: }
077: }
078: rows[numRows] = n;
079: }
080:
081: /**
082: * Compares two double sparse square matrices for equality.
083: * @param m a double matrix
084: */
085: public boolean equals(AbstractDoubleSquareMatrix m, double tol) {
086: if (numRows == m.numRows && numCols == m.numCols) {
087: if (m instanceof DoubleSparseSquareMatrix) {
088: return this .equals((DoubleSparseSquareMatrix) m);
089: } else {
090: double sumSqr = 0;
091: for (int i = 0; i < numRows; i++) {
092: for (int j = 0; j < numCols; j++) {
093: double delta = getElement(i, j)
094: - m.getElement(i, j);
095: sumSqr += delta * delta;
096: }
097: }
098: return (sumSqr <= tol * tol);
099: }
100: } else
101: return false;
102: }
103:
104: public final boolean equals(DoubleSparseSquareMatrix m) {
105: return equals(m, GlobalSettings.ZERO_TOL);
106: }
107:
108: public boolean equals(DoubleSparseSquareMatrix m, double tol) {
109: if (numRows == m.numRows && numCols == m.numCols) {
110: if (colPos.length != m.colPos.length)
111: return false;
112: for (int i = 1; i < rows.length; i++) {
113: if (rows[i] != m.rows[i])
114: return false;
115: }
116: double sumSqr = 0.0;
117: for (int i = 1; i < colPos.length; i++) {
118: if (colPos[i] != m.colPos[i])
119: return false;
120: double delta = elements[i] - m.elements[i];
121: sumSqr += delta * delta;
122: }
123: return (sumSqr <= tol * tol);
124: } else
125: return false;
126: }
127:
128: /**
129: * Returns a string representing this matrix.
130: */
131: public String toString() {
132: final StringBuffer buf = new StringBuffer(numRows * numCols);
133: for (int i = 0; i < numRows; i++) {
134: for (int j = 0; j < numCols; j++) {
135: buf.append(getElement(i, j));
136: buf.append(' ');
137: }
138: buf.append('\n');
139: }
140: return buf.toString();
141: }
142:
143: /**
144: * Converts this matrix to an integer matrix.
145: * @return an integer square matrix
146: */
147: public AbstractIntegerMatrix toIntegerMatrix() {
148: final int ans[][] = new int[numRows][numCols];
149: for (int i = 0; i < numRows; i++) {
150: for (int j = 0; j < numCols; j++)
151: ans[i][j] = Math.round((float) getElement(i, j));
152: }
153: return new IntegerSquareMatrix(ans);
154: }
155:
156: /**
157: * Converts this matrix to a complex matrix.
158: * @return a complex square matrix
159: */
160: public AbstractComplexMatrix toComplexMatrix() {
161: final double arrayRe[][] = new double[numRows][numCols];
162: for (int i = 0; i < numRows; i++) {
163: for (int j = 0; j < numCols; j++)
164: arrayRe[i][j] = getElement(i, j);
165: }
166: return new ComplexSquareMatrix(arrayRe,
167: new double[numRows][numCols]);
168: }
169:
170: /**
171: * Returns an element of the matrix.
172: * @param i row index of the element
173: * @param j column index of the element
174: * @exception MatrixDimensionException If attempting to access an invalid element.
175: */
176: public double getElement(final int i, final int j) {
177: if (i >= 0 && i < numRows && j >= 0 && j < numCols) {
178: for (int k = rows[i]; k < rows[i + 1]; k++) {
179: if (colPos[k] == j)
180: return elements[k];
181: }
182: return 0.0;
183: } else
184: throw new MatrixDimensionException(getInvalidElementMsg(i,
185: j));
186: }
187:
188: /**
189: * Sets the value of an element of the matrix.
190: * @param i row index of the element
191: * @param j column index of the element
192: * @param x a number
193: * @exception MatrixDimensionException If attempting to access an invalid element.
194: */
195: public void setElement(final int i, final int j, final double x) {
196: if (i >= 0 && i < numRows && j >= 0 && j < numCols) {
197: if (Math.abs(x) <= GlobalSettings.ZERO_TOL)
198: return;
199: for (int k = rows[i]; k < rows[i + 1]; k++) {
200: if (colPos[k] == j) {
201: elements[k] = x;
202: return;
203: }
204: }
205: final double oldMatrix[] = elements;
206: final int oldColPos[] = colPos;
207: elements = new double[oldMatrix.length + 1];
208: colPos = new int[oldColPos.length + 1];
209: System.arraycopy(oldMatrix, 0, elements, 0, rows[i]);
210: System.arraycopy(oldColPos, 0, colPos, 0, rows[i]);
211: int k;
212: for (k = rows[i]; k < rows[i + 1] && oldColPos[k] < j; k++) {
213: elements[k] = oldMatrix[k];
214: colPos[k] = oldColPos[k];
215: }
216: elements[k] = x;
217: colPos[k] = j;
218: System.arraycopy(oldMatrix, k, elements, k + 1,
219: oldMatrix.length - k);
220: System.arraycopy(oldColPos, k, colPos, k + 1,
221: oldColPos.length - k);
222: for (k = i + 1; k < rows.length; k++)
223: rows[k]++;
224: } else
225: throw new MatrixDimensionException(getInvalidElementMsg(i,
226: j));
227: }
228:
229: /**
230: * Returns the determinant.
231: */
232: public double det() {
233: final AbstractDoubleSquareMatrix lu[] = this .luDecompose(null);
234: double det = lu[1].getElement(0, 0);
235: for (int i = 1; i < numRows; i++)
236: det *= lu[1].getElement(i, i);
237: return det * LUpivot[numRows];
238: }
239:
240: /**
241: * Returns the trace.
242: */
243: public double trace() {
244: double result = getElement(0, 0);
245: for (int i = 1; i < numRows; i++)
246: result += getElement(i, i);
247: return result;
248: }
249:
250: /**
251: * Returns the l<sup><img border=0 alt="infinity" src="doc-files/infinity.gif"></sup>-norm.
252: */
253: public double infNorm() {
254: double result = 0.0, tmpResult;
255: for (int j, i = 0; i < numRows; i++) {
256: tmpResult = 0.0;
257: for (j = rows[i]; j < rows[i + 1]; j++)
258: tmpResult += Math.abs(elements[j]);
259: if (tmpResult > result)
260: result = tmpResult;
261: }
262: return result;
263: }
264:
265: /**
266: * Returns the Frobenius (l<sup>2</sup>) norm.
267: */
268: public double frobeniusNorm() {
269: double result = 0.0;
270: for (int i = 0; i < colPos.length; i++)
271: result += elements[i] * elements[i];
272: return Math.sqrt(result);
273: }
274:
275: //============
276: // OPERATIONS
277: //============
278:
279: // ADDITION
280:
281: /**
282: * Returns the addition of this matrix and another.
283: * @param m a double matrix
284: * @exception MatrixDimensionException If the matrices are different sizes.
285: */
286: public AbstractDoubleSquareMatrix add(
287: final AbstractDoubleSquareMatrix m) {
288: if (m instanceof DoubleSparseSquareMatrix)
289: return add((DoubleSparseSquareMatrix) m);
290: if (m instanceof DoubleSquareMatrix)
291: return add((DoubleSquareMatrix) m);
292:
293: if (numRows == m.rows() && numCols == m.columns()) {
294: final double array[][] = new double[numRows][numCols];
295: for (int i = 0; i < numRows; i++) {
296: for (int j = rows[i]; j < rows[i + 1]; j++)
297: array[i][colPos[j]] = elements[j];
298: array[i][0] += m.getElement(i, 0);
299: for (int j = 1; j < numCols; j++)
300: array[i][j] += m.getElement(i, j);
301: }
302: return new DoubleSquareMatrix(array);
303: } else {
304: throw new MatrixDimensionException(
305: "Matrices are different sizes.");
306: }
307: }
308:
309: public DoubleSquareMatrix add(final DoubleSquareMatrix m) {
310: if (numRows == m.numRows && numCols == m.numCols) {
311: final double array[][] = new double[numRows][numCols];
312: for (int i = 0; i < numRows; i++) {
313: for (int j = rows[i]; j < rows[i + 1]; j++)
314: array[i][colPos[j]] = elements[j];
315: array[i][0] += m.matrix[i][0];
316: for (int j = 1; j < numCols; j++)
317: array[i][j] += m.matrix[i][j];
318: }
319: return new DoubleSquareMatrix(array);
320: } else
321: throw new MatrixDimensionException(
322: "Matrices are different sizes.");
323: }
324:
325: /**
326: * Returns the addition of this matrix and another.
327: * @param m a double sparse matrix
328: * @exception MatrixDimensionException If the matrices are different sizes.
329: */
330: public DoubleSparseSquareMatrix add(final DoubleSparseSquareMatrix m) {
331: if (numRows == m.numRows && numCols == m.numCols) {
332: DoubleSparseSquareMatrix ans = new DoubleSparseSquareMatrix(
333: numRows);
334: for (int i = 0; i < numRows; i++) {
335: for (int j = 0; j < numCols; j++)
336: ans.setElement(i, j, getElement(i, j)
337: + m.getElement(i, j));
338: }
339: return ans;
340: } else
341: throw new MatrixDimensionException(
342: "Matrices are different sizes.");
343: }
344:
345: // SUBTRACTION
346:
347: /**
348: * Returns the subtraction of this matrix and another.
349: * @param m a double matrix
350: * @exception MatrixDimensionException If the matrices are different sizes.
351: */
352: public AbstractDoubleSquareMatrix subtract(
353: final AbstractDoubleSquareMatrix m) {
354: if (m instanceof DoubleSparseSquareMatrix)
355: return subtract((DoubleSparseSquareMatrix) m);
356: if (m instanceof DoubleSquareMatrix)
357: return subtract((DoubleSquareMatrix) m);
358:
359: if (numRows == m.rows() && numCols == m.columns()) {
360: final double array[][] = new double[numRows][numCols];
361: for (int i = 0; i < numRows; i++) {
362: for (int j = rows[i]; j < rows[i + 1]; j++)
363: array[i][colPos[j]] = elements[j];
364: array[i][0] -= m.getElement(i, 0);
365: for (int j = 1; j < numCols; j++)
366: array[i][j] -= m.getElement(i, j);
367: }
368: return new DoubleSquareMatrix(array);
369: } else {
370: throw new MatrixDimensionException(
371: "Matrices are different sizes.");
372: }
373: }
374:
375: public DoubleSquareMatrix subtract(final DoubleSquareMatrix m) {
376: if (numRows == m.numRows && numCols == m.numCols) {
377: final double array[][] = new double[numRows][numCols];
378: for (int i = 0; i < numRows; i++) {
379: for (int j = rows[i]; j < rows[i + 1]; j++)
380: array[i][colPos[j]] = elements[j];
381: array[i][0] -= m.matrix[i][0];
382: for (int j = 1; j < numCols; j++)
383: array[i][j] -= m.matrix[i][j];
384: }
385: return new DoubleSquareMatrix(array);
386: } else
387: throw new MatrixDimensionException(
388: "Matrices are different sizes.");
389: }
390:
391: /**
392: * Returns the addition of this matrix and another.
393: * @param m a double sparse matrix
394: * @exception MatrixDimensionException If the matrices are different sizes.
395: */
396: public DoubleSparseSquareMatrix subtract(
397: final DoubleSparseSquareMatrix m) {
398: if (numRows == m.numRows && numCols == m.numCols) {
399: DoubleSparseSquareMatrix ans = new DoubleSparseSquareMatrix(
400: numRows);
401: for (int i = 0; i < numRows; i++) {
402: for (int j = 0; j < numCols; j++)
403: ans.setElement(i, j, getElement(i, j)
404: - m.getElement(i, j));
405: }
406: return ans;
407: } else
408: throw new MatrixDimensionException(
409: "Matrices are different sizes.");
410: }
411:
412: // SCALAR MULTIPLICATION
413:
414: /**
415: * Returns the multiplication of this matrix by a scalar.
416: * @param x a double
417: * @return a double sparse matrix
418: */
419: public AbstractDoubleMatrix scalarMultiply(final double x) {
420: final DoubleSparseSquareMatrix ans = new DoubleSparseSquareMatrix(
421: numRows);
422: ans.elements = new double[elements.length];
423: ans.colPos = new int[colPos.length];
424: System.arraycopy(colPos, 0, ans.colPos, 0, colPos.length);
425: System.arraycopy(rows, 0, ans.rows, 0, rows.length);
426: for (int i = 0; i < colPos.length; i++)
427: ans.elements[i] = x * elements[i];
428: return ans;
429: }
430:
431: public AbstractDoubleMatrix scalarDivide(final double x) {
432: final DoubleSparseSquareMatrix ans = new DoubleSparseSquareMatrix(
433: numRows);
434: ans.elements = new double[elements.length];
435: ans.colPos = new int[colPos.length];
436: System.arraycopy(colPos, 0, ans.colPos, 0, colPos.length);
437: System.arraycopy(rows, 0, ans.rows, 0, rows.length);
438: for (int i = 0; i < colPos.length; i++)
439: ans.elements[i] = elements[i] / x;
440: return ans;
441: }
442:
443: // SCALAR PRODUCT
444:
445: /**
446: * Returns the scalar product of this matrix and another.
447: * @param m a double matrix.
448: * @exception MatrixDimensionException If the matrices are different sizes.
449: */
450: public double scalarProduct(final AbstractDoubleSquareMatrix m) {
451: if (m instanceof DoubleSquareMatrix)
452: return scalarProduct((DoubleSquareMatrix) m);
453:
454: if (numRows == m.numRows && numCols == m.numCols) {
455: double ans = 0.0;
456: for (int i = 0; i < numRows; i++) {
457: for (int j = rows[i]; j < rows[i + 1]; j++)
458: ans += elements[j] * m.getElement(i, colPos[j]);
459: }
460: return ans;
461: } else
462: throw new MatrixDimensionException(
463: "Matrices are different sizes.");
464: }
465:
466: public double scalarProduct(final DoubleSquareMatrix m) {
467: if (numRows == m.numRows && numCols == m.numCols) {
468: double ans = 0.0;
469: for (int i = 0; i < numRows; i++) {
470: for (int j = rows[i]; j < rows[i + 1]; j++)
471: ans += elements[j] * m.matrix[i][colPos[j]];
472: }
473: return ans;
474: } else
475: throw new MatrixDimensionException(
476: "Matrices are different sizes.");
477: }
478:
479: // MATRIX MULTIPLICATION
480:
481: /**
482: * Returns the multiplication of a vector by this matrix.
483: * @param v a double vector
484: * @exception DimensionException If the matrix and vector are incompatible.
485: */
486: public AbstractDoubleVector multiply(final AbstractDoubleVector v) {
487: if (numCols == v.dimension()) {
488: final double array[] = new double[numRows];
489: for (int i = 0; i < numRows; i++) {
490: for (int j = rows[i]; j < rows[i + 1]; j++)
491: array[i] += elements[j] * v.getComponent(colPos[j]);
492: }
493: return new DoubleVector(array);
494: } else
495: throw new DimensionException(
496: "Matrix and vector are incompatible.");
497: }
498:
499: /**
500: * Returns the multiplication of this matrix and another.
501: * @param m a double matrix
502: * @exception MatrixDimensionException If the matrices are incompatible.
503: */
504: public AbstractDoubleSquareMatrix multiply(
505: final AbstractDoubleSquareMatrix m) {
506: if (m instanceof DoubleSparseSquareMatrix)
507: return multiply((DoubleSparseSquareMatrix) m);
508: if (m instanceof DoubleSquareMatrix)
509: return multiply((DoubleSquareMatrix) m);
510:
511: if (numCols == m.numRows) {
512: final double array[][] = new double[numRows][m.numCols];
513: for (int j = 0; j < numRows; j++) {
514: for (int k = 0; k < m.numCols; k++) {
515: for (int n = rows[j]; n < rows[j + 1]; n++)
516: array[j][k] += elements[n]
517: * m.getElement(colPos[n], k);
518: }
519: }
520: return new DoubleSquareMatrix(array);
521: } else
522: throw new MatrixDimensionException("Incompatible matrices.");
523: }
524:
525: public DoubleSquareMatrix multiply(final DoubleSquareMatrix m) {
526: if (numCols == m.numRows) {
527: final double array[][] = new double[numRows][m.numCols];
528: for (int j = 0; j < numRows; j++) {
529: for (int k = 0; k < m.numCols; k++) {
530: for (int n = rows[j]; n < rows[j + 1]; n++)
531: array[j][k] += elements[n]
532: * m.matrix[colPos[n]][k];
533: }
534: }
535: return new DoubleSquareMatrix(array);
536: } else
537: throw new MatrixDimensionException("Incompatible matrices.");
538: }
539:
540: /**
541: * Returns the multiplication of this matrix and another.
542: * @param m a double sparse matrix
543: * @exception MatrixDimensionException If the matrices are incompatible.
544: */
545: public DoubleSparseSquareMatrix multiply(
546: final DoubleSparseSquareMatrix m) {
547: if (numCols == m.numRows) {
548: DoubleSparseSquareMatrix ans = new DoubleSparseSquareMatrix(
549: numRows);
550: for (int j = 0; j < numRows; j++) {
551: for (int k = 0; k < numCols; k++) {
552: double tmp = 0.0;
553: for (int n = rows[j]; n < rows[j + 1]; n++)
554: tmp += elements[n] * m.getElement(colPos[n], k);
555: ans.setElement(j, k, tmp);
556: }
557: }
558: return ans;
559: } else
560: throw new MatrixDimensionException("Incompatible matrices.");
561: }
562:
563: // TRANSPOSE
564:
565: /**
566: * Returns the transpose of this matrix.
567: * @return a double sparse matrix
568: */
569: public Matrix transpose() {
570: final DoubleSparseSquareMatrix ans = new DoubleSparseSquareMatrix(
571: numRows);
572: for (int i = 0; i < numRows; i++) {
573: ans.setElement(0, i, getElement(i, 0));
574: for (int j = 1; j < numCols; j++)
575: ans.setElement(j, i, getElement(i, j));
576: }
577: return ans;
578: }
579:
580: // LU DECOMPOSITION
581:
582: /**
583: * Returns the LU decomposition of this matrix.
584: * @param pivot an empty array of length <code>rows()+1</code>
585: * to hold the pivot information (null if not interested)
586: * @return an array with [0] containing the L-matrix
587: * and [1] containing the U-matrix.
588: */
589: public AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
590: if (LU != null) {
591: if (pivot != null)
592: System.arraycopy(LUpivot, 0, pivot, 0, pivot.length);
593: return LU;
594: }
595: final double arrayL[][] = new double[numRows][numCols];
596: final double arrayU[][] = new double[numRows][numCols];
597: if (pivot == null)
598: pivot = new int[numRows + 1];
599: for (int i = 0; i < numRows; i++)
600: pivot[i] = i;
601: pivot[numRows] = 1;
602: // LU decomposition to arrayU
603: for (int j = 0; j < numCols; j++) {
604: for (int i = 0; i < j; i++) {
605: double tmp = getElement(pivot[i], j);
606: for (int k = 0; k < i; k++)
607: tmp -= arrayU[i][k] * arrayU[k][j];
608: arrayU[i][j] = tmp;
609: }
610: double max = 0.0;
611: int pivotrow = j;
612: for (int i = j; i < numRows; i++) {
613: double tmp = getElement(pivot[i], j);
614: for (int k = 0; k < j; k++)
615: tmp -= arrayU[i][k] * arrayU[k][j];
616: arrayU[i][j] = tmp;
617: // while we're here search for a pivot for arrayU[j][j]
618: tmp = Math.abs(tmp);
619: if (tmp > max) {
620: max = tmp;
621: pivotrow = i;
622: }
623: }
624: // swap row j with pivotrow
625: if (pivotrow != j) {
626: double[] tmprow = arrayU[j];
627: arrayU[j] = arrayU[pivotrow];
628: arrayU[pivotrow] = tmprow;
629: int k = pivot[j];
630: pivot[j] = pivot[pivotrow];
631: pivot[pivotrow] = k;
632: // update parity
633: pivot[numRows] = -pivot[numRows];
634: }
635: // divide by pivot
636: for (int i = j + 1; i < numRows; i++)
637: arrayU[i][j] /= arrayU[j][j];
638: }
639: // move lower triangular part to arrayL
640: for (int j = 0; j < numCols; j++) {
641: arrayL[j][j] = 1.0;
642: for (int i = j + 1; i < numRows; i++) {
643: arrayL[i][j] = arrayU[i][j];
644: arrayU[i][j] = 0.0;
645: }
646: }
647: LU = new AbstractDoubleSquareMatrix[2];
648: LU[0] = new DoubleSquareMatrix(arrayL);
649: LU[1] = new DoubleSquareMatrix(arrayU);
650: LUpivot = new int[pivot.length];
651: System.arraycopy(pivot, 0, LUpivot, 0, pivot.length);
652: return LU;
653: }
654:
655: /**
656: * Returns the LU decomposition of this matrix.
657: * Warning: no pivoting.
658: * @return an array with [0] containing the L-matrix
659: * and [1] containing the U-matrix.
660: * @jsci.planetmath LUDecomposition
661: */
662: public AbstractDoubleSquareMatrix[] luDecompose() {
663: final double arrayL[][] = new double[numRows][numCols];
664: final double arrayU[][] = new double[numRows][numCols];
665: // LU decomposition to arrayU
666: for (int j = 0; j < numCols; j++) {
667: for (int i = 0; i < j; i++) {
668: double tmp = getElement(i, j);
669: for (int k = 0; k < i; k++)
670: tmp -= arrayU[i][k] * arrayU[k][j];
671: arrayU[i][j] = tmp;
672: }
673: for (int i = j; i < numRows; i++) {
674: double tmp = getElement(i, j);
675: for (int k = 0; k < j; k++)
676: tmp -= arrayU[i][k] * arrayU[k][j];
677: arrayU[i][j] = tmp;
678: }
679: // divide
680: for (int i = j + 1; i < numRows; i++)
681: arrayU[i][j] /= arrayU[j][j];
682: }
683: // move lower triangular part to arrayL
684: for (int j = 0; j < numCols; j++) {
685: arrayL[j][j] = 1.0;
686: for (int i = j + 1; i < numRows; i++) {
687: arrayL[i][j] = arrayU[i][j];
688: arrayU[i][j] = 0.0;
689: }
690: }
691: AbstractDoubleSquareMatrix[] lu = new AbstractDoubleSquareMatrix[2];
692: lu[0] = new DoubleSquareMatrix(arrayL);
693: lu[1] = new DoubleSquareMatrix(arrayU);
694: return lu;
695: }
696:
697: // CHOLESKY DECOMPOSITION
698:
699: /**
700: * Returns the Cholesky decomposition of this matrix.
701: * Matrix must be symmetric and positive definite.
702: * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
703: */
704: public AbstractDoubleSquareMatrix[] choleskyDecompose() {
705: final double arrayL[][] = new double[numRows][numCols];
706: final double arrayU[][] = new double[numRows][numCols];
707: arrayL[0][0] = arrayU[0][0] = Math.sqrt(getElement(0, 0));
708: for (int i = 1; i < numRows; i++)
709: arrayL[i][0] = arrayU[0][i] = getElement(i, 0)
710: / arrayL[0][0];
711: for (int j = 1; j < numCols; j++) {
712: double tmp = getElement(j, j);
713: for (int i = 0; i < j; i++)
714: tmp -= arrayL[j][i] * arrayL[j][i];
715: arrayL[j][j] = arrayU[j][j] = Math.sqrt(tmp);
716: for (int i = j + 1; i < numRows; i++) {
717: tmp = getElement(i, j);
718: for (int k = 0; k < i; k++)
719: tmp -= arrayL[j][k] * arrayU[k][i];
720: arrayL[i][j] = arrayU[j][i] = tmp / arrayU[j][j];
721: }
722: }
723: final AbstractDoubleSquareMatrix lu[] = new AbstractDoubleSquareMatrix[2];
724: lu[0] = new DoubleSquareMatrix(arrayL);
725: lu[1] = new DoubleSquareMatrix(arrayU);
726: return lu;
727: }
728:
729: // MAP ELEMENTS
730:
731: /**
732: * Applies a function on all the matrix elements.
733: * @param f a user-defined function
734: * @return a double sparse matrix
735: */
736: public AbstractDoubleMatrix mapElements(final Mapping f) {
737: final DoubleSparseSquareMatrix ans = new DoubleSparseSquareMatrix(
738: numRows);
739: ans.elements = new double[elements.length];
740: ans.colPos = new int[colPos.length];
741: System.arraycopy(colPos, 0, ans.colPos, 0, colPos.length);
742: System.arraycopy(rows, 0, ans.rows, 0, rows.length);
743: for (int i = 0; i < colPos.length; i++)
744: ans.elements[i] = f.map(elements[i]);
745: return ans;
746: }
747: }
|