DoubleSparseSquareMatrix


1   package JSci.maths.matrices;
2   
3   import JSci.GlobalSettings;
4   import JSci.maths.Mapping;
5   import JSci.maths.DimensionException;
6   import JSci.maths.vectors.AbstractDoubleVector;
7   import JSci.maths.vectors.DoubleVector;
8   
9   /**
10  * The DoubleSparseSquareMatrix class provides an object for encapsulating sparse square matrices.
11  * Uses compressed row storage.
12  * @version 1.2
13  * @author Mark Hale
14  */
15  public final class DoubleSparseSquareMatrix extends AbstractDoubleSquareMatrix {
16          /**
17          * Matrix elements.
18          */
19          private double elements[];
20          /**
21          * Sparse indexing data.
22          * Contains the column positions of each element,
23          * e.g. <code>colPos[n]</code> is the column position
24          * of the <code>n</code>th element.
25          */
26          private int colPos[];
27          /**
28          * Sparse indexing data.
29          * Contains the indices of the start of each row,
30          * e.g. <code>rows[i]</code> is the index
31          * where the <code>i</code>th row starts.
32          */
33          private int rows[];
34          /**
35          * Constructs an empty matrix.
36          * @param size the number of rows/columns
37          */
38          public DoubleSparseSquareMatrix(final int size) {
39                  super(size);
40                  elements=new double[0];
41                  colPos=new int[0];
42                  rows=new int[numRows+1];
43          }
44          /**
45          * Constructs a matrix from an array.
46          * @param array an assigned value
47          * @exception MatrixDimensionException If the array is not square.
48          */
49          public DoubleSparseSquareMatrix(final double array[][]) {
50                  super(array.length);
51                  rows=new int[numRows+1];
52                  int n=0;
53                  for(int i=0;i<numRows;i++) {
54                          if(array[i].length != array.length)
55                                  throw new MatrixDimensionException("Array is not square.");
56                          for(int j=0;j<numCols;j++) {
57                                  if(Math.abs(array[i][j])>GlobalSettings.ZERO_TOL)
58                                          n++;
59                          }
60                  }
61                  elements=new double[n];
62                  colPos=new int[n];
63                  n=0;
64                  for(int i=0;i<numRows;i++) {
65                          rows[i]=n;
66                          for(int j=0;j<numCols;j++) {
67                                  if(Math.abs(array[i][j])>GlobalSettings.ZERO_TOL) {
68                                          elements[n]=array[i][j];
69                                          colPos[n]=j;
70                                          n++;
71                                  }
72                          }
73                  }
74                  rows[numRows]=n;
75          }
76          /**
77          * Compares two double sparse square matrices for equality.
78          * @param m a double matrix
79          */
80          public boolean equals(AbstractDoubleSquareMatrix m, double tol) {
81                  if(numRows==m.numRows && numCols==m.numCols) {
82                          if(m instanceof DoubleSparseSquareMatrix) {
83                                  return this.equals((DoubleSparseSquareMatrix)m);
84                          } else {
85                      double sumSqr = 0;
86                                  for(int i=0;i<numRows;i++) {
87                                          for(int j=0;j<numCols;j++) {
88                              double delta = getElement(i,j)-m.getElement(i,j);
89                              sumSqr += delta*delta;
90                                          }
91                                  }
92                                  return (sumSqr <= tol*tol);
93                          }
94                  } else
95                          return false;
96          }
97          public final boolean equals(DoubleSparseSquareMatrix m) {
98          return equals(m, GlobalSettings.ZERO_TOL);
99      }
100         public boolean equals(DoubleSparseSquareMatrix m, double tol) {
101                 if(numRows==m.numRows && numCols==m.numCols) {
102                         if(colPos.length!=m.colPos.length)
103                                 return false;
104                         for(int i=1;i<rows.length;i++) {
105                                 if(rows[i]!=m.rows[i])
106                                         return false;
107                         }
108                         double sumSqr = 0.0;
109                         for(int i=1;i<colPos.length;i++) {
110                                 if(colPos[i]!=m.colPos[i])
111                                         return false;
112                                 double delta = elements[i] - m.elements[i];
113                                 sumSqr += delta*delta;
114                         }
115                         return (sumSqr <= tol*tol);
116                 } else
117                         return false;
118         }
119         /**
120         * Returns a string representing this matrix.
121         */
122         public String   toString() {
123                 final StringBuffer   buf=new StringBuffer  (numRows*numCols);
124                 for(int i=0;i<numRows;i++) {
125                         for(int j=0;j<numCols;j++) {
126                                 buf.append(getElement(i,j));
127                                 buf.append(' ');
128                         }
129                         buf.append('\n');
130                 }
131                 return buf.toString();
132         }
133         /**
134         * Converts this matrix to an integer matrix.
135         * @return an integer square matrix
136         */
137         public AbstractIntegerMatrix toIntegerMatrix() {
138                 final int ans[][]=new int[numRows][numCols];
139                 for(int i=0;i<numRows;i++) {
140                         for(int j=0;j<numCols;j++)
141                                 ans[i][j]=Math.round((float)getElement(i,j));
142                 }
143                 return new IntegerSquareMatrix(ans);
144         }
145         /**
146         * Converts this matrix to a complex matrix.
147         * @return a complex square matrix
148         */
149         public AbstractComplexMatrix toComplexMatrix() {
150                 final double arrayRe[][]=new double[numRows][numCols];
151                 for(int i=0;i<numRows;i++) {
152                         for(int j=0;j<numCols;j++)
153                                 arrayRe[i][j]=getElement(i,j);
154                 }
155                 return new ComplexSquareMatrix(arrayRe,new double[numRows][numCols]);
156         }
157         /**
158         * Returns an element of the matrix.
159         * @param i row index of the element
160         * @param j column index of the element
161         * @exception MatrixDimensionException If attempting to access an invalid element.
162         */
163         public double getElement(final int i, final int j) {
164                 if(i>=0 && i<numRows && j>=0 && j<numCols) {
165                         for(int k=rows[i];k<rows[i+1];k++) {
166                                 if(colPos[k]==j)
167                                         return elements[k];
168                         }
169                         return 0.0;
170                 } else
171                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
172         }
173         /**
174         * Sets the value of an element of the matrix.
175         * @param i row index of the element
176         * @param j column index of the element
177         * @param x a number
178         * @exception MatrixDimensionException If attempting to access an invalid element.
179         */
180         public void setElement(final int i, final int j, final double x) {
181                 if(i>=0 && i<numRows && j>=0 && j<numCols) {
182                         if(Math.abs(x)<=GlobalSettings.ZERO_TOL)
183                                 return;
184                         for(int k=rows[i];k<rows[i+1];k++) {
185                                 if(colPos[k]==j) {
186                                         elements[k]=x;
187                                         return;
188                                 }
189                         }
190                         final double oldMatrix[]=elements;
191                         final int oldColPos[]=colPos;
192                         elements=new double[oldMatrix.length+1];
193                         colPos=new int[oldColPos.length+1];
194                         System.arraycopy(oldMatrix,0,elements,0,rows[i]);
195                         System.arraycopy(oldColPos,0,colPos,0,rows[i]);
196                         int k;
197                         for(k=rows[i];k<rows[i+1] && oldColPos[k]<j;k++) {
198                                 elements[k]=oldMatrix[k];
199                                 colPos[k]=oldColPos[k];
200                         }
201                         elements[k]=x;
202                         colPos[k]=j;
203                         System.arraycopy(oldMatrix,k,elements,k+1,oldMatrix.length-k);
204                         System.arraycopy(oldColPos,k,colPos,k+1,oldColPos.length-k);
205                         for(k=i+1;k<rows.length;k++)
206                                 rows[k]++;
207                 } else
208                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
209         }
210         /**
211         * Returns the determinant.
212         */
213         public double det() {
214                 final AbstractDoubleSquareMatrix lu[]=this.luDecompose(null);
215                 double det=lu[1].getElement(0,0);
216                 for(int i=1;i<numRows;i++)
217                         det*=lu[1].getElement(i,i);
218                 return det*LUpivot[numRows];
219         }
220         /**
221         * Returns the trace.
222         */
223         public double trace() {
224                 double result=getElement(0,0);
225                 for(int i=1;i<numRows;i++)
226                         result+=getElement(i,i);
227                 return result;
228         }
229         /**
230         * Returns the l<sup><img border=0 alt="infinity" SRC="doc-files/infinity.gif"></sup>-norm.
231         */
232         public double infNorm() {
233                 double result=0.0,tmpResult;
234                 for(int j,i=0;i<numRows;i++) {
235                         tmpResult=0.0;
236                         for(j=rows[i];j<rows[i+1];j++)
237                                 tmpResult+=Math.abs(elements[j]);
238                         if(tmpResult>result)
239                                 result=tmpResult;
240                 }
241                 return result;
242         }
243         /**
244         * Returns the Frobenius (l<sup>2</sup>) norm.
245         */
246         public double frobeniusNorm() {
247                 double result=0.0;
248                 for(int i=0;i<colPos.length;i++)
249                         result+=elements[i]*elements[i];
250                 return Math.sqrt(result);
251         }
252 
253 //============
254 // OPERATIONS
255 //============
256 
257 // ADDITION
258 
259         /**
260         * Returns the addition of this matrix and another.
261         * @param m a double matrix
262         * @exception MatrixDimensionException If the matrices are different sizes.
263         */
264         public AbstractDoubleSquareMatrix add(final AbstractDoubleSquareMatrix m) {
265                 if(m instanceof DoubleSparseSquareMatrix)
266                         return add((DoubleSparseSquareMatrix)m);
267                 if(m instanceof DoubleSquareMatrix)
268                         return add((DoubleSquareMatrix)m);
269 
270                 if(numRows==m.rows() && numCols==m.columns()) {
271                         final double array[][]=new double[numRows][numCols];
272                         for(int i=0;i<numRows;i++) {
273                                 for(int j=rows[i];j<rows[i+1];j++)
274                                         array[i][colPos[j]]=elements[j];
275                                 array[i][0]+=m.getElement(i,0);
276                                 for(int j=1;j<numCols;j++)
277                                         array[i][j]+=m.getElement(i,j);
278                         }
279                         return new DoubleSquareMatrix(array);
280                 } else {
281                         throw new MatrixDimensionException("Matrices are different sizes.");
282                 }
283         }
284         public DoubleSquareMatrix add(final DoubleSquareMatrix m) {
285                 if(numRows==m.numRows && numCols==m.numCols) {
286                         final double array[][]=new double[numRows][numCols];
287                         for(int i=0;i<numRows;i++) {
288                                 for(int j=rows[i];j<rows[i+1];j++)
289                                         array[i][colPos[j]]=elements[j];
290                                 array[i][0]+=m.matrix[i][0];
291                                 for(int j=1;j<numCols;j++)
292                                         array[i][j]+=m.matrix[i][j];
293                         }
294                         return new DoubleSquareMatrix(array);
295                 } else
296                         throw new MatrixDimensionException("Matrices are different sizes.");
297         }
298         /**
299         * Returns the addition of this matrix and another.
300         * @param m a double sparse matrix
301         * @exception MatrixDimensionException If the matrices are different sizes.
302         */
303         public DoubleSparseSquareMatrix add(final DoubleSparseSquareMatrix m) {
304                 if(numRows==m.numRows && numCols==m.numCols) {
305                         DoubleSparseSquareMatrix ans=new DoubleSparseSquareMatrix(numRows);
306                         for(int i=0;i<numRows;i++) {
307                                 for(int j=0;j<numCols;j++)
308                                         ans.setElement(i,j,getElement(i,j)+m.getElement(i,j));
309                         }
310                         return ans;
311                 } else
312                         throw new MatrixDimensionException("Matrices are different sizes.");
313         }
314 
315 // SUBTRACTION
316 
317         /**
318         * Returns the subtraction of this matrix and another.
319         * @param m a double matrix
320         * @exception MatrixDimensionException If the matrices are different sizes.
321         */
322         public AbstractDoubleSquareMatrix subtract(final AbstractDoubleSquareMatrix m) {
323                 if(m instanceof DoubleSparseSquareMatrix)
324                         return subtract((DoubleSparseSquareMatrix)m);
325                 if(m instanceof DoubleSquareMatrix)
326                         return subtract((DoubleSquareMatrix)m);
327 
328                 if(numRows==m.rows() && numCols==m.columns()) {
329                         final double array[][]=new double[numRows][numCols];
330                         for(int i=0;i<numRows;i++) {
331                                 for(int j=rows[i];j<rows[i+1];j++)
332                                         array[i][colPos[j]]=elements[j];
333                                 array[i][0]-=m.getElement(i,0);
334                                 for(int j=1;j<numCols;j++)
335                                         array[i][j]-=m.getElement(i,j);
336                         }
337                         return new DoubleSquareMatrix(array);
338                 } else {
339                         throw new MatrixDimensionException("Matrices are different sizes.");
340                 }
341         }
342         public DoubleSquareMatrix subtract(final DoubleSquareMatrix m) {
343                 if(numRows==m.numRows && numCols==m.numCols) {
344                         final double array[][]=new double[numRows][numCols];
345                         for(int i=0;i<numRows;i++) {
346                                 for(int j=rows[i];j<rows[i+1];j++)
347                                         array[i][colPos[j]]=elements[j];
348                                 array[i][0]-=m.matrix[i][0];
349                                 for(int j=1;j<numCols;j++)
350                                         array[i][j]-=m.matrix[i][j];
351                         }
352                         return new DoubleSquareMatrix(array);
353                 } else
354                         throw new MatrixDimensionException("Matrices are different sizes.");
355         }
356         /**
357         * Returns the addition of this matrix and another.
358         * @param m a double sparse matrix
359         * @exception MatrixDimensionException If the matrices are different sizes.
360         */
361         public DoubleSparseSquareMatrix subtract(final DoubleSparseSquareMatrix m) {
362                 if(numRows==m.numRows && numCols==m.numCols) {
363                         DoubleSparseSquareMatrix ans=new DoubleSparseSquareMatrix(numRows);
364                         for(int i=0;i<numRows;i++) {
365                                 for(int j=0;j<numCols;j++)
366                                         ans.setElement(i,j,getElement(i,j)-m.getElement(i,j));
367                         }
368                         return ans;
369                 } else
370                         throw new MatrixDimensionException("Matrices are different sizes.");
371         }
372 
373 // SCALAR MULTIPLICATION
374 
375         /**
376         * Returns the multiplication of this matrix by a scalar.
377         * @param x a double
378         * @return a double sparse matrix
379         */
380         public AbstractDoubleMatrix scalarMultiply(final double x) {
381                 final DoubleSparseSquareMatrix ans=new DoubleSparseSquareMatrix(numRows);
382                 ans.elements=new double[elements.length];
383                 ans.colPos=new int[colPos.length];
384                 System.arraycopy(colPos,0,ans.colPos,0,colPos.length);
385                 System.arraycopy(rows,0,ans.rows,0,rows.length);
386                 for(int i=0;i<colPos.length;i++)
387                         ans.elements[i]=x*elements[i];
388                 return ans;
389         }
390 
391         public AbstractDoubleMatrix scalarDivide(final double x) {
392                 final DoubleSparseSquareMatrix ans=new DoubleSparseSquareMatrix(numRows);
393                 ans.elements=new double[elements.length];
394                 ans.colPos=new int[colPos.length];
395                 System.arraycopy(colPos,0,ans.colPos,0,colPos.length);
396                 System.arraycopy(rows,0,ans.rows,0,rows.length);
397                 for(int i=0;i<colPos.length;i++)
398                         ans.elements[i]=elements[i]/x;
399                 return ans;
400         }
401 
402 // SCALAR PRODUCT
403 
404         /**
405         * Returns the scalar product of this matrix and another.
406         * @param m a double matrix.
407         * @exception MatrixDimensionException If the matrices are different sizes.
408         */
409         public double scalarProduct(final AbstractDoubleSquareMatrix m) {
410                 if(m instanceof DoubleSquareMatrix)
411                         return scalarProduct((DoubleSquareMatrix)m);
412 
413                 if(numRows==m.numRows && numCols==m.numCols) {
414                         double ans=0.0;
415                         for(int i=0;i<numRows;i++) {
416                                 for(int j=rows[i];j<rows[i+1];j++)
417                                         ans+=elements[j]*m.getElement(i,colPos[j]);
418                         }
419                         return ans;
420                 } else
421                         throw new MatrixDimensionException("Matrices are different sizes.");
422         }
423         public double scalarProduct(final DoubleSquareMatrix m) {
424                 if(numRows==m.numRows && numCols==m.numCols) {
425                         double ans=0.0;
426                         for(int i=0;i<numRows;i++) {
427                                 for(int j=rows[i];j<rows[i+1];j++)
428                                         ans+=elements[j]*m.matrix[i][colPos[j]];
429                         }
430                         return ans;
431                 } else
432                         throw new MatrixDimensionException("Matrices are different sizes.");
433         }
434 
435 // MATRIX MULTIPLICATION
436 
437         /**
438         * Returns the multiplication of a vector by this matrix.
439         * @param v a double vector
440         * @exception DimensionException If the matrix and vector are incompatible.
441         */
442         public AbstractDoubleVector multiply(final AbstractDoubleVector v) {
443                 if(numCols==v.dimension()) {
444                         final double array[]=new double[numRows];
445                         for(int i=0;i<numRows;i++) {
446                                 for(int j=rows[i];j<rows[i+1];j++)
447                                         array[i]+=elements[j]*v.getComponent(colPos[j]);
448                         }
449                         return new DoubleVector(array);
450                 } else
451                         throw new DimensionException("Matrix and vector are incompatible.");
452         }
453         /**
454         * Returns the multiplication of this matrix and another.
455         * @param m a double matrix
456         * @exception MatrixDimensionException If the matrices are incompatible.
457         */
458         public AbstractDoubleSquareMatrix multiply(final AbstractDoubleSquareMatrix m) {
459                 if(m instanceof DoubleSparseSquareMatrix)
460                         return multiply((DoubleSparseSquareMatrix)m);
461                 if(m instanceof DoubleSquareMatrix)
462                         return multiply((DoubleSquareMatrix)m);
463 
464                 if(numCols==m.numRows) {
465                         final double array[][]=new double[numRows][m.numCols];
466                         for(int j=0;j<numRows;j++) {
467                                 for(int k=0;k<m.numCols;k++) {
468                                         for(int n=rows[j];n<rows[j+1];n++)
469                                                 array[j][k]+=elements[n]*m.getElement(colPos[n],k);
470                                 }
471                         }
472                         return new DoubleSquareMatrix(array);
473                 } else
474                         throw new MatrixDimensionException("Incompatible matrices.");
475         }
476         public DoubleSquareMatrix multiply(final DoubleSquareMatrix m) {
477                 if(numCols==m.numRows) {
478                         final double array[][]=new double[numRows][m.numCols];
479                         for(int j=0;j<numRows;j++) {
480                                 for(int k=0;k<m.numCols;k++) {
481                                         for(int n=rows[j];n<rows[j+1];n++)
482                                                 array[j][k]+=elements[n]*m.matrix[colPos[n]][k];
483                                 }
484                         }
485                         return new DoubleSquareMatrix(array);
486                 } else
487                         throw new MatrixDimensionException("Incompatible matrices.");
488         }
489         /**
490         * Returns the multiplication of this matrix and another.
491         * @param m a double sparse matrix
492         * @exception MatrixDimensionException If the matrices are incompatible.
493         */
494         public DoubleSparseSquareMatrix multiply(final DoubleSparseSquareMatrix m) {
495                 if(numCols==m.numRows) {
496                         DoubleSparseSquareMatrix ans=new DoubleSparseSquareMatrix(numRows);
497                         for(int j=0;j<numRows;j++) {
498                                 for(int k=0;k<numCols;k++) {
499                                         double tmp=0.0;
500                                         for(int n=rows[j];n<rows[j+1];n++)
501                                                 tmp+=elements[n]*m.getElement(colPos[n],k);
502                                         ans.setElement(j,k,tmp);
503                                 }
504                         }
505                         return ans;
506                 } else
507                         throw new MatrixDimensionException("Incompatible matrices.");
508         }
509 
510 // TRANSPOSE
511 
512         /**
513         * Returns the transpose of this matrix.
514         * @return a double sparse matrix
515         */
516         public Matrix transpose() {
517                 final DoubleSparseSquareMatrix ans=new DoubleSparseSquareMatrix(numRows);
518                 for(int i=0;i<numRows;i++) {
519                         ans.setElement(0,i,getElement(i,0));
520                         for(int j=1;j<numCols;j++)
521                                 ans.setElement(j,i,getElement(i,j));
522                 }
523                 return ans;
524         }
525 
526 // LU DECOMPOSITION
527 
528         /**
529         * Returns the LU decomposition of this matrix.
530         * @param pivot an empty array of length <code>rows()+1</code>
531         * to hold the pivot information (null if not interested)
532         * @return an array with [0] containing the L-matrix
533         * and [1] containing the U-matrix.
534         */
535         public AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
536                 if(LU!=null) {
537                         if(pivot!=null)
538                                 System.arraycopy(LUpivot,0,pivot,0,pivot.length);
539                         return LU;
540                 }
541                 final double arrayL[][]=new double[numRows][numCols];
542                 final double arrayU[][]=new double[numRows][numCols];
543                 if(pivot==null)
544                         pivot=new int[numRows+1];
545                 for(int i=0;i<numRows;i++)
546                         pivot[i]=i;
547                 pivot[numRows]=1;
548         // LU decomposition to arrayU
549                 for(int j=0;j<numCols;j++) {
550                         for(int i=0;i<j;i++) {
551                                 double tmp=getElement(pivot[i],j);
552                                 for(int k=0;k<i;k++)
553                                         tmp-=arrayU[i][k]*arrayU[k][j];
554                                 arrayU[i][j]=tmp;
555                         }
556                         double max=0.0;
557                         int pivotrow=j;
558                         for(int i=j;i<numRows;i++) {
559                                 double tmp=getElement(pivot[i],j);
560                                 for(int k=0;k<j;k++)
561                                         tmp-=arrayU[i][k]*arrayU[k][j];
562                                 arrayU[i][j]=tmp;
563                         // while we're here search for a pivot for arrayU[j][j]
564                                 tmp=Math.abs(tmp);
565                                 if(tmp>max) {
566                                         max=tmp;
567                                         pivotrow=i;
568                                 }
569                         }
570                 // swap row j with pivotrow
571                         if(pivotrow!=j) {
572                                 double[] tmprow = arrayU[j];
573                                 arrayU[j] = arrayU[pivotrow];
574                                 arrayU[pivotrow] = tmprow;
575                                 int k=pivot[j];
576                                 pivot[j]=pivot[pivotrow];
577                                 pivot[pivotrow]=k;
578                                 // update parity
579                                 pivot[numRows]=-pivot[numRows];
580                         }
581                 // divide by pivot
582                         for(int i=j+1;i<numRows;i++)
583                                 arrayU[i][j]/=arrayU[j][j];
584                 }
585                 // move lower triangular part to arrayL
586                 for(int j=0;j<numCols;j++) {
587                         arrayL[j][j]=1.0;
588                         for(int i=j+1;i<numRows;i++) {
589                                 arrayL[i][j]=arrayU[i][j];
590                                 arrayU[i][j]=0.0;
591                         }
592                 }
593                 LU=new AbstractDoubleSquareMatrix[2];
594                 LU[0]=new DoubleSquareMatrix(arrayL);
595                 LU[1]=new DoubleSquareMatrix(arrayU);
596                 LUpivot=new int[pivot.length];
597                 System.arraycopy(pivot,0,LUpivot,0,pivot.length);
598                 return LU;
599         }
600         /**
601         * Returns the LU decomposition of this matrix.
602         * Warning: no pivoting.
603         * @return an array with [0] containing the L-matrix
604         * and [1] containing the U-matrix.
605         * @jsci.planetmath LUDecomposition
606         */
607         public AbstractDoubleSquareMatrix[] luDecompose() {
608                 final double arrayL[][]=new double[numRows][numCols];
609                 final double arrayU[][]=new double[numRows][numCols];
610         // LU decomposition to arrayU
611                 for(int j=0;j<numCols;j++) {
612                         for(int i=0;i<j;i++) {
613                                 double tmp=getElement(i,j);
614                                 for(int k=0;k<i;k++)
615                                         tmp-=arrayU[i][k]*arrayU[k][j];
616                                 arrayU[i][j]=tmp;
617                         }
618                         for(int i=j;i<numRows;i++) {
619                                 double tmp=getElement(i,j);
620                                 for(int k=0;k<j;k++)
621                                         tmp-=arrayU[i][k]*arrayU[k][j];
622                                 arrayU[i][j]=tmp;
623                         }
624                 // divide
625                         for(int i=j+1;i<numRows;i++)
626                                 arrayU[i][j]/=arrayU[j][j];
627                 }
628                 // move lower triangular part to arrayL
629                 for(int j=0;j<numCols;j++) {
630                         arrayL[j][j]=1.0;
631                         for(int i=j+1;i<numRows;i++) {
632                                 arrayL[i][j]=arrayU[i][j];
633                                 arrayU[i][j]=0.0;
634                         }
635                 }
636                 AbstractDoubleSquareMatrix[] lu=new AbstractDoubleSquareMatrix[2];
637                 lu[0]=new DoubleSquareMatrix(arrayL);
638                 lu[1]=new DoubleSquareMatrix(arrayU);
639                 return lu;
640         }
641 
642 // CHOLESKY DECOMPOSITION
643 
644         /**
645         * Returns the Cholesky decomposition of this matrix.
646         * Matrix must be symmetric and positive definite.
647         * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
648         */
649         public AbstractDoubleSquareMatrix[] choleskyDecompose() {
650                 final double arrayL[][]=new double[numRows][numCols];
651                 final double arrayU[][]=new double[numRows][numCols];
652                 arrayL[0][0]=arrayU[0][0]=Math.sqrt(getElement(0,0));
653                 for(int i=1;i<numRows;i++)
654                         arrayL[i][0]=arrayU[0][i]=getElement(i,0)/arrayL[0][0];
655                 for(int j=1;j<numCols;j++) {
656                         double tmp=getElement(j,j);
657                         for(int i=0;i<j;i++)
658                                 tmp-=arrayL[j][i]*arrayL[j][i];
659                         arrayL[j][j]=arrayU[j][j]=Math.sqrt(tmp);
660                         for(int i=j+1;i<numRows;i++) {
661                                 tmp=getElement(i,j);
662                                 for(int k=0;k<i;k++)
663                                         tmp-=arrayL[j][k]*arrayU[k][i];
664                                 arrayL[i][j]=arrayU[j][i]=tmp/arrayU[j][j];
665                         }
666                 }
667                 final AbstractDoubleSquareMatrix lu[]=new AbstractDoubleSquareMatrix[2];
668                 lu[0]=new DoubleSquareMatrix(arrayL);
669                 lu[1]=new DoubleSquareMatrix(arrayU);
670                 return lu;
671         }
672 
673 // MAP ELEMENTS
674 
675         /**
676         * Applies a function on all the matrix elements.
677         * @param f a user-defined function
678         * @return a double sparse matrix
679         */
680         public AbstractDoubleMatrix mapElements(final Mapping f) {
681                 final DoubleSparseSquareMatrix ans=new DoubleSparseSquareMatrix(numRows);
682                 ans.elements=new double[elements.length];
683                 ans.colPos=new int[colPos.length];
684                 System.arraycopy(colPos,0,ans.colPos,0,colPos.length);
685                 System.arraycopy(rows,0,ans.rows,0,rows.length);
686                 for(int i=0;i<colPos.length;i++)
687                         ans.elements[i]=f.map(elements[i]);
688                 return ans;
689         }
690 }
691 
692
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags