KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > JSci > maths > matrices > DoubleTridiagonalMatrix


1 /* AUTO-GENERATED */
2 package JSci.maths.matrices;
3
4 import JSci.maths.ArrayMath;
5 import JSci.maths.ExtraMath;
6 import JSci.maths.Mapping;
7 import JSci.maths.LinearMath;
8 import JSci.maths.MaximumIterationsExceededException;
9 import JSci.maths.DimensionException;
10 import JSci.maths.vectors.AbstractDoubleVector;
11 import JSci.maths.vectors.DoubleVector;
12 import JSci.maths.groups.AbelianGroup;
13 import JSci.maths.algebras.*;
14 import JSci.maths.fields.*;
15
16 /**
17 * The DoubleTridiagonalMatrix class provides an object for encapsulating double tridiagonal matrices.
18 * @version 2.2
19 * @author Mark Hale
20 */
21 public class DoubleTridiagonalMatrix extends AbstractDoubleSquareMatrix implements TridiagonalMatrix {
22         /**
23         * Tridiagonal data.
24         */
25         protected final double ldiag[];
26         protected final double diag[];
27         protected final double udiag[];
28         /**
29         * Constructs an empty matrix.
30         * @param size the number of rows/columns
31         */
32         public DoubleTridiagonalMatrix(final int size) {
33                 super(size);
34                 ldiag=new double[size];
35                 diag=new double[size];
36                 udiag=new double[size];
37         }
38         /**
39         * Constructs a matrix from an array.
40         * Any non-tridiagonal elements in the array are ignored.
41         * @param array an assigned value
42         * @exception MatrixDimensionException If the array is not square.
43         */
44         public DoubleTridiagonalMatrix(final double array[][]) {
45                 this(array.length);
46                 if(!ArrayMath.isSquare(array))
47                         throw new MatrixDimensionException("Array is not square.");
48                 diag[0]=array[0][0];
49                 udiag[0]=array[0][1];
50                 int i=1;
51                 for(;i<array.length-1;i++) {
52                         ldiag[i]=array[i][i-1];
53                         diag[i]=array[i][i];
54                         udiag[i]=array[i][i+1];
55                 }
56                 ldiag[i]=array[i][i-1];
57                 diag[i]=array[i][i];
58         }
59         /**
60         * Compares two ${nativeTyp} matrices for equality.
61         * @param m a double matrix
62         */
63         public boolean equals(AbstractDoubleMatrix m, double tol) {
64                 if(m instanceof TridiagonalMatrix) {
65                         if(numRows != m.rows() || numCols != m.columns())
66                                 return false;
67             double sumSqr = 0;
68             double ldelta;
69             double delta = diag[0] - m.getElement(0,0);
70             double udelta = udiag[0] - m.getElement(0,1);
71             sumSqr += delta*delta + udelta*udelta;
72                         int i=1;
73                         for(;i<numRows-1;i++) {
74                 ldelta = ldiag[i] - m.getElement(i,i-1);
75                 delta = diag[i] - m.getElement(i,i);
76                 udelta = udiag[i] - m.getElement(i,i+1);
77                 sumSqr += ldelta*ldelta + delta*delta + udelta*udelta;
78                         }
79             ldelta = ldiag[i] - m.getElement(i,i-1);
80             delta = diag[i] - m.getElement(i,i);
81             sumSqr += ldelta*ldelta + delta*delta;
82                         return (sumSqr <= tol*tol);
83                 } else {
84                         return false;
85                 }
86         }
87         /**
88         * Returns a string representing this matrix.
89         */
90         public String JavaDoc toString() {
91                 final StringBuffer JavaDoc buf=new StringBuffer JavaDoc(5*numRows*numCols);
92                 for(int i=0;i<numRows;i++) {
93                         for(int j=0;j<numCols;j++) {
94                                 buf.append(getElement(i,j));
95                                 buf.append(' ');
96                         }
97                         buf.append('\n');
98                 }
99                 return buf.toString();
100         }
101         /**
102         * Converts this matrix to an integer matrix.
103         * @return an integer matrix
104         */
105         public AbstractIntegerMatrix toIntegerMatrix() {
106                 final IntegerTridiagonalMatrix m=new IntegerTridiagonalMatrix(numRows);
107                 m.diag[0]=Math.round((float)diag[0]);
108                 m.udiag[0]=Math.round((float)udiag[0]);
109                 int i=1;
110                 for(;i<numRows-1;i++) {
111                         m.ldiag[i]=Math.round((float)ldiag[i]);
112                         m.diag[i]=Math.round((float)diag[i]);
113                         m.udiag[i]=Math.round((float)udiag[i]);
114                 }
115                 m.ldiag[i]=Math.round((float)ldiag[i]);
116                 m.diag[i]=Math.round((float)diag[i]);
117                 return m;
118         }
119         /**
120         * Converts this matrix to a complex matrix.
121         * @return a complex matrix
122         */
123         public AbstractComplexMatrix toComplexMatrix() {
124                 final ComplexTridiagonalMatrix m=new ComplexTridiagonalMatrix(numRows);
125                 m.diagRe[0]=diag[0];
126                 m.udiagRe[0]=udiag[0];
127                 int i=1;
128                 for(;i<numRows-1;i++) {
129                         m.ldiagRe[i]=ldiag[i];
130                         m.diagRe[i]=diag[i];
131                         m.udiagRe[i]=udiag[i];
132                 }
133                 m.ldiagRe[i]=ldiag[i];
134                 m.diagRe[i]=diag[i];
135                 return m;
136         }
137         /**
138         * Returns an element of the matrix.
139         * @param i row index of the element
140         * @param j column index of the element
141         * @exception MatrixDimensionException If attempting to access an invalid element.
142         */
143         public double getElement(int i, int j) {
144                 if(i>=0 && i<numRows && j>=0 && j<numCols) {
145                         if(j == i-1)
146                                 return ldiag[i];
147                         else if(j == i)
148                                 return diag[i];
149                         else if(j == i+1)
150                                 return udiag[i];
151                         else
152                                 return 0;
153                 } else
154                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
155         }
156         /**
157         * Sets the value of an element of the matrix.
158         * Should only be used to initialise this matrix.
159         * @param i row index of the element
160         * @param j column index of the element
161         * @param x a number
162         * @exception MatrixDimensionException If attempting to access an invalid element.
163         */
164         public void setElement(int i, int j, final double x) {
165                 if(i>=0 && i<numRows && j>=0 && j<numCols) {
166                         if(j == i-1)
167                                 ldiag[i] = x;
168                         else if(j == i)
169                                 diag[i] = x;
170                         else if(j == i+1)
171                                 udiag[i] = x;
172                         else
173                                 throw new MatrixDimensionException(getInvalidElementMsg(i,j));
174                 } else
175                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
176         }
177         /**
178         * Returns true if this matrix is symmetric.
179         */
180         public boolean isSymmetric() {
181                 if(ldiag[1]!=udiag[0])
182                         return false;
183                 for(int i=1;i<numRows-1;i++) {
184                         if(ldiag[i+1]!=udiag[i])
185                                 return false;
186                 }
187                 return true;
188         }
189         /**
190         * Returns the trace.
191         */
192         public double trace() {
193                 double tr=diag[0];
194                 for(int i=1;i<numRows;i++)
195                         tr+=diag[i];
196                 return tr;
197         }
198         /**
199         * Returns the l<sup><img border=0 alt="infinity" SRC="doc-files/infinity.gif"></sup>-norm.
200         * @author Taber Smith
201         */
202         public double infNorm() {
203                 double result=Math.abs(diag[0])+Math.abs(udiag[0]);
204                 double tmpResult;
205                 int i=1;
206                 for(;i<numRows-1;i++) {
207                         tmpResult=Math.abs(ldiag[i])+Math.abs(diag[i])+Math.abs(udiag[i]);
208                         if(tmpResult>result)
209                                 result=tmpResult;
210                 }
211                 tmpResult=Math.abs(ldiag[i])+Math.abs(diag[i]);
212                 if(tmpResult>result)
213                         result=tmpResult;
214                 return result;
215         }
216         /**
217         * Returns the Frobenius (l<sup>2</sup>) norm.
218         * @author Taber Smith
219         */
220         public double frobeniusNorm() {
221                 double result=diag[0]*diag[0]+udiag[0]*udiag[0];
222                 int i=1;
223                 for(;i<numRows-1;i++)
224                         result+=ldiag[i]*ldiag[i]+diag[i]*diag[i]+udiag[i]*udiag[i];
225                 result+=ldiag[i]*ldiag[i]+diag[i]*diag[i];
226                 return Math.sqrt(result);
227         }
228         /**
229         * Returns the operator norm.
230         * @exception MaximumIterationsExceededException If it takes more than 50 iterations to determine an eigenvalue.
231         */
232         public double operatorNorm() throws MaximumIterationsExceededException {
233                 return Math.sqrt(ArrayMath.max(LinearMath.eigenvalueSolveSymmetric((DoubleTridiagonalMatrix)(this.transpose().multiply(this)))));
234         }
235
236 //============
237 // OPERATIONS
238 //============
239
240 // ADDITION
241
242         /**
243         * Returns the addition of this matrix and another.
244         * @param m a double matrix
245         * @exception MatrixDimensionException If the matrices are different sizes.
246         */
247         public AbstractDoubleSquareMatrix add(final AbstractDoubleSquareMatrix m) {
248                 if(m instanceof DoubleTridiagonalMatrix)
249                         return add((DoubleTridiagonalMatrix)m);
250                 if(m instanceof TridiagonalMatrix)
251                         return addTridiagonal(m);
252                 if(m instanceof DoubleSquareMatrix)
253                         return add((DoubleSquareMatrix)m);
254
255                 if(numRows==m.rows() && numCols==m.columns()) {
256                         final double array[][]=new double[numRows][numCols];
257                         for(int i=0;i<numRows;i++) {
258                                 array[i][0]=getElement(i,0)+m.getElement(i,0);
259                                 for(int j=1;j<numCols;j++)
260                                         array[i][j]=getElement(i,j)+m.getElement(i,j);
261                         }
262                         return new DoubleSquareMatrix(array);
263                 } else {
264                         throw new MatrixDimensionException("Matrices are different sizes.");
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                         array[0][0]+=diag[0];
273                         array[0][1]+=udiag[0];
274                         int n=numRows-1;
275                         for(int i=1;i<n;i++) {
276                                 array[i][i-1]+=ldiag[i];
277                                 array[i][i]+=diag[i];
278                                 array[i][i+1]+=udiag[i];
279                         }
280                         array[n][n-1]+=ldiag[n];
281                         array[n][n]+=diag[n];
282                         return new DoubleSquareMatrix(array);
283                 } else
284                         throw new MatrixDimensionException("Matrices are different sizes.");
285         }
286         /**
287         * Returns the addition of this matrix and another.
288         * @param m a double tridiagonal matrix
289         * @exception MatrixDimensionException If the matrices are different sizes.
290         */
291         public DoubleTridiagonalMatrix add(final DoubleTridiagonalMatrix m) {
292                 int mRow=numRows;
293                 if(mRow==m.numRows) {
294                         final DoubleTridiagonalMatrix ans=new DoubleTridiagonalMatrix(mRow);
295                         ans.diag[0]=diag[0]+m.diag[0];
296                         ans.udiag[0]=udiag[0]+m.udiag[0];
297                         mRow--;
298                         for(int i=1;i<mRow;i++) {
299                                 ans.ldiag[i]=ldiag[i]+m.ldiag[i];
300                                 ans.diag[i]=diag[i]+m.diag[i];
301                                 ans.udiag[i]=udiag[i]+m.udiag[i];
302                         }
303                         ans.ldiag[mRow]=ldiag[mRow]+m.ldiag[mRow];
304                         ans.diag[mRow]=diag[mRow]+m.diag[mRow];
305                         return ans;
306                 } else
307                         throw new MatrixDimensionException("Matrices are different sizes.");
308         }
309         private DoubleTridiagonalMatrix addTridiagonal(final AbstractDoubleSquareMatrix m) {
310                 int mRow=numRows;
311                 if(mRow==m.rows()) {
312                         final DoubleTridiagonalMatrix ans=new DoubleTridiagonalMatrix(mRow);
313                         ans.diag[0]=diag[0]+m.getElement(0,0);
314                         ans.udiag[0]=udiag[0]+m.getElement(0,1);
315                         mRow--;
316                         for(int i=1;i<mRow;i++) {
317                                 ans.ldiag[i]=ldiag[i]+m.getElement(i,i-1);
318                                 ans.diag[i]=diag[i]+m.getElement(i,i);
319                                 ans.udiag[i]=udiag[i]+m.getElement(i,i+1);
320                         }
321                         ans.ldiag[mRow]=ldiag[mRow]+m.getElement(mRow,mRow-1);
322                         ans.diag[mRow]=diag[mRow]+m.getElement(mRow,mRow);
323                         return ans;
324                 } else {
325                         throw new MatrixDimensionException("Matrices are different sizes.");
326                 }
327         }
328
329 // SUBTRACTION
330
331         /**
332         * Returns the subtraction of this matrix by another.
333         * @param m a double matrix
334         * @exception MatrixDimensionException If the matrices are different sizes.
335         */
336         public AbstractDoubleSquareMatrix subtract(final AbstractDoubleSquareMatrix m) {
337                 if(m instanceof DoubleTridiagonalMatrix)
338                         return subtract((DoubleTridiagonalMatrix)m);
339                 if(m instanceof TridiagonalMatrix)
340                         return subtractTridiagonal(m);
341                 if(m instanceof DoubleSquareMatrix)
342                         return subtract((DoubleSquareMatrix)m);
343
344                 if(numRows==m.rows() && numCols==m.columns()) {
345                         final double array[][]=new double[numRows][numCols];
346                         for(int i=0;i<numRows;i++) {
347                                 array[i][0]=getElement(i,0)-m.getElement(i,0);
348                                 for(int j=1;j<numCols;j++)
349                                         array[i][j]=getElement(i,j)-m.getElement(i,j);
350                         }
351                         return new DoubleSquareMatrix(array);
352                 } else {
353                         throw new MatrixDimensionException("Matrices are different sizes.");
354                 }
355         }
356         public DoubleSquareMatrix subtract(final DoubleSquareMatrix m) {
357                 if(numRows==m.numRows && numCols==m.numCols) {
358                         final double array[][]=new double[numRows][numCols];
359                         for(int i=0;i<numRows;i++) {
360                                 array[i][0]=getElement(i,0)-m.matrix[i][0];
361                                 for(int j=1;j<numCols;j++)
362                                         array[i][j]=getElement(i,j)-m.matrix[i][j];
363                         }
364                         return new DoubleSquareMatrix(array);
365                 } else
366                         throw new MatrixDimensionException("Matrices are different sizes.");
367         }
368         /**
369         * Returns the subtraction of this matrix and another.
370         * @param m a double tridiagonal matrix
371         * @exception MatrixDimensionException If the matrices are different sizes.
372         */
373         public DoubleTridiagonalMatrix subtract(final DoubleTridiagonalMatrix m) {
374                 int mRow=numRows;
375                 if(mRow==m.numRows) {
376                         final DoubleTridiagonalMatrix ans=new DoubleTridiagonalMatrix(mRow);
377                         ans.diag[0]=diag[0]-m.diag[0];
378                         ans.udiag[0]=udiag[0]-m.udiag[0];
379                         mRow--;
380                         for(int i=1;i<mRow;i++) {
381                                 ans.ldiag[i]=ldiag[i]-m.ldiag[i];
382                                 ans.diag[i]=diag[i]-m.diag[i];
383                                 ans.udiag[i]=udiag[i]-m.udiag[i];
384                         }
385                         ans.ldiag[mRow]=ldiag[mRow]-m.ldiag[mRow];
386                         ans.diag[mRow]=diag[mRow]-m.diag[mRow];
387                         return ans;
388                 } else
389                         throw new MatrixDimensionException("Matrices are different sizes.");
390         }
391         private DoubleTridiagonalMatrix subtractTridiagonal(final AbstractDoubleSquareMatrix m) {
392                 int mRow=numRows;
393                 if(mRow==m.rows()) {
394                         final DoubleTridiagonalMatrix ans=new DoubleTridiagonalMatrix(mRow);
395                         ans.diag[0]=diag[0]-m.getElement(0,0);
396                         ans.udiag[0]=udiag[0]-m.getElement(0,1);
397                         mRow--;
398                         for(int i=1;i<mRow;i++) {
399                                 ans.ldiag[i]=ldiag[i]-m.getElement(i,i-1);
400                                 ans.diag[i]=diag[i]-m.getElement(i,i);
401                                 ans.udiag[i]=udiag[i]-m.getElement(i,i+1);
402                         }
403                         ans.ldiag[mRow]=ldiag[mRow]-m.getElement(mRow,mRow-1);
404                         ans.diag[mRow]=diag[mRow]-m.getElement(mRow,mRow);
405                         return ans;
406                 } else {
407                         throw new MatrixDimensionException("Matrices are different sizes.");
408                 }
409         }
410
411 // SCALAR MULTIPLICATION
412
413         /**
414         * Returns the multiplication of this matrix by a scalar.
415         * @param x a double.
416         * @return a double square matrix.
417         */
418         public AbstractDoubleMatrix scalarMultiply(final double x) {
419                 final DoubleTridiagonalMatrix ans=new DoubleTridiagonalMatrix(numRows);
420                 final int lastRow = numRows-1;
421                 ans.diag[0]=x*diag[0];
422                 ans.udiag[0]=x*udiag[0];
423                 for(int i=1;i<lastRow;i++) {
424                         ans.ldiag[i]=x*ldiag[i];
425                         ans.diag[i]=x*diag[i];
426                         ans.udiag[i]=x*udiag[i];
427                 }
428                 ans.ldiag[lastRow]=x*ldiag[lastRow];
429                 ans.diag[lastRow]=x*diag[lastRow];
430                 return ans;
431         }
432
433 // SCALAR DIVISON
434
435         /**
436         * Returns the division of this matrix by a scalar.
437         * @param x a double.
438         * @return a double square matrix.
439         */
440         public AbstractDoubleMatrix scalarDivide(final double x) {
441                 final DoubleTridiagonalMatrix ans=new DoubleTridiagonalMatrix(numRows);
442                 final int lastRow = numRows-1;
443                 ans.diag[0]=diag[0]/x;
444                 ans.udiag[0]=udiag[0]/x;
445                 for(int i=1;i<lastRow;i++) {
446                         ans.ldiag[i]=ldiag[i]/x;
447                         ans.diag[i]=diag[i]/x;
448                         ans.udiag[i]=udiag[i]/x;
449                 }
450                 ans.ldiag[lastRow]=ldiag[lastRow]/x;
451                 ans.diag[lastRow]=diag[lastRow]/x;
452                 return ans;
453         }
454
455 // SCALAR PRODUCT
456
457         /**
458         * Returns the scalar product of this matrix and another.
459         * @param m a double matrix.
460         * @exception MatrixDimensionException If the matrices are different sizes.
461         */
462         public double scalarProduct(final AbstractDoubleSquareMatrix m) {
463                 if(m instanceof DoubleTridiagonalMatrix)
464                         return scalarProduct((DoubleTridiagonalMatrix)m);
465                 if(m instanceof DoubleSquareMatrix)
466                         return scalarProduct((DoubleSquareMatrix)m);
467
468                 if(numRows==m.rows() && numCols==m.columns()) {
469                         double ans = diag[0]*m.getElement(0,0)+udiag[0]*m.getElement(0,1);
470                         final int lastRow = numRows-1;
471                         for(int i=1;i<lastRow;i++)
472                                 ans += ldiag[i]*m.getElement(i,i-1)+diag[i]*m.getElement(i,i)+udiag[i]*m.getElement(i,i+1);
473                         ans += ldiag[lastRow]*m.getElement(lastRow,lastRow-1)+diag[lastRow]*m.getElement(lastRow,lastRow);
474                         return ans;
475                 } else {
476                        throw new MatrixDimensionException("Matrices are different sizes.");
477                 }
478         }
479         public double scalarProduct(final DoubleSquareMatrix m) {
480                 if(numRows==m.numRows && numCols==m.numCols) {
481                         double ans = diag[0]*m.matrix[0][0]+udiag[0]*m.matrix[0][1];
482                         final int lastRow = numRows-1;
483                         for(int i=1;i<lastRow;i++)
484                                 ans += ldiag[i]*m.matrix[i][i-1]+diag[i]*m.matrix[i][i]+udiag[i]*m.matrix[i][i+1];
485                         ans += ldiag[lastRow]*m.matrix[lastRow][lastRow-1]+diag[lastRow]*m.matrix[lastRow][lastRow];
486                         return ans;
487                 } else
488                         throw new MatrixDimensionException("Matrices are different sizes.");
489         }
490         public double scalarProduct(final DoubleTridiagonalMatrix m) {
491                 if(numRows==m.numRows) {
492                         double ans = diag[0]*m.diag[0]+udiag[0]*m.udiag[0];
493                         final int lastRow = numRows-1;
494                         for(int i=1;i<lastRow;i++)
495                                 ans += ldiag[i]*m.ldiag[i]+diag[i]*m.diag[i]+udiag[i]*m.udiag[i];
496                         ans += ldiag[lastRow]*m.ldiag[lastRow]+diag[lastRow]*m.diag[lastRow];
497                         return ans;
498                 } else
499                         throw new MatrixDimensionException("Matrices are different sizes.");
500         }
501
502 // MATRIX MULTIPLICATION
503
504         /**
505         * Returns the multiplication of a vector by this matrix.
506         * @param v a double vector.
507         * @exception DimensionException If the matrix and vector are incompatible.
508         */
509         public AbstractDoubleVector multiply(final AbstractDoubleVector v) {
510                 if(numCols==v.dimension()) {
511                         final double array[]=new double[numRows];
512                         final int lastRow = numRows-1;
513                         array[0]=diag[0]*v.getComponent(0)+udiag[0]*v.getComponent(1);
514                         for(int i=1;i<lastRow;i++)
515                                 array[i]=ldiag[i]*v.getComponent(i-1)+diag[i]*v.getComponent(i)+udiag[i]*v.getComponent(i+1);
516                         array[lastRow]=ldiag[lastRow]*v.getComponent(lastRow-1)+diag[lastRow]*v.getComponent(lastRow);
517                         return new DoubleVector(array);
518                 } else {
519                         throw new DimensionException("Matrix and vector are incompatible.");
520                 }
521         }
522         /**
523         * Returns the multiplication of this matrix and another.
524         * @param m a double matrix
525         * @return a AbstractDoubleMatrix or a AbstractDoubleSquareMatrix as appropriate
526         * @exception MatrixDimensionException If the matrices are incompatible.
527         */
528         public AbstractDoubleSquareMatrix multiply(final AbstractDoubleSquareMatrix m) {
529                 if(m instanceof DoubleTridiagonalMatrix)
530                         return multiply((DoubleTridiagonalMatrix)m);
531                 if(m instanceof TridiagonalMatrix)
532                         return multiplyTridiagonal(m);
533                 if(m instanceof DoubleSquareMatrix)
534                         return multiply((DoubleSquareMatrix)m);
535
536                 if(numCols==m.rows()) {
537                         final int mColumns = m.columns();
538                         final double array[][]=new double[numRows][mColumns];
539                         final int lastRow = numRows-1;
540                         for(int j=0; j<numRows; j++) {
541                                 array[0][j]=diag[0]*m.getElement(0,j)+udiag[0]*m.getElement(1,j);
542                                 for(int i=1;i<lastRow;i++)
543                                         array[i][j]=ldiag[i]*m.getElement(i-1,j)+diag[i]*m.getElement(i,j)+udiag[i]*m.getElement(i+1,j);
544                                 array[lastRow][j]=ldiag[lastRow]*m.getElement(lastRow-1,j)+diag[lastRow]*m.getElement(lastRow,j);
545                         }
546                         return new DoubleSquareMatrix(array);
547                 } else {
548                         throw new MatrixDimensionException("Incompatible matrices.");
549                 }
550         }
551         public DoubleSquareMatrix multiply(final DoubleSquareMatrix m) {
552                 if(numCols==m.numRows) {
553                         final double array[][]=new double[numRows][m.numCols];
554                         final int lastRow = numRows-1;
555                         for(int j=0;j<numRows;j++) {
556                                 array[0][j]=diag[0]*m.matrix[0][j]+udiag[0]*m.matrix[1][j];
557                                 for(int i=1;i<lastRow;i++)
558                                         array[i][j]=ldiag[i]*m.matrix[i-1][j]+diag[i]*m.matrix[i][j]+udiag[i]*m.matrix[i+1][j];
559                                 array[lastRow][j]=ldiag[lastRow]*m.matrix[lastRow-1][j]+diag[lastRow]*m.matrix[lastRow][j];
560                         }
561                         return new DoubleSquareMatrix(array);
562                 } else
563                         throw new MatrixDimensionException("Incompatible matrices.");
564         }
565         public DoubleSquareMatrix multiply(final DoubleTridiagonalMatrix m) {
566                 int mRow=numRows;
567                 if(numCols==m.numRows) {
568                         final double array[][]=new double[mRow][mRow];
569                         array[0][0]=diag[0]*m.diag[0]+udiag[0]*m.ldiag[1];
570                         array[0][1]=diag[0]*m.udiag[0]+udiag[0]*m.diag[1];
571                         array[0][2]=udiag[0]*m.udiag[1];
572                         if(mRow>3) {
573                                 array[1][0]=ldiag[1]*m.diag[0]+diag[1]*m.ldiag[1];
574                                 array[1][1]=ldiag[1]*m.udiag[0]+diag[1]*m.diag[1]+udiag[1]*m.ldiag[2];
575                                 array[1][2]=diag[1]*m.udiag[1]+udiag[1]*m.diag[2];
576                                 array[1][3]=udiag[1]*m.udiag[2];
577                         }
578                         if(mRow==3) {
579                                 array[1][0]=ldiag[1]*m.diag[0]+diag[1]*m.ldiag[1];
580                                 array[1][1]=ldiag[1]*m.udiag[0]+diag[1]*m.diag[1]+udiag[1]*m.ldiag[2];
581                                 array[1][2]=diag[1]*m.udiag[1]+udiag[1]*m.diag[2];
582                         } else if(mRow>4) {
583                                 for(int i=2;i<mRow-2;i++) {
584                                         array[i][i-2]=ldiag[i]*m.ldiag[i-1];
585                                         array[i][i-1]=ldiag[i]*m.diag[i-1]+diag[i]*m.ldiag[i];
586                                         array[i][i]=ldiag[i]*m.udiag[i-1]+diag[i]*m.diag[i]+udiag[i]*m.ldiag[i+1];
587                                         array[i][i+1]=diag[i]*m.udiag[i]+udiag[i]*m.diag[i+1];
588                                         array[i][i+2]=udiag[i]*m.udiag[i+1];
589                                 }
590                         }
591                         if(mRow>3) {
592                                 array[mRow-2][mRow-4]=ldiag[mRow-2]*m.ldiag[mRow-3];
593                                 array[mRow-2][mRow-3]=ldiag[mRow-2]*m.diag[mRow-3]+diag[mRow-2]*m.ldiag[mRow-2];
594                                 array[mRow-2][mRow-2]=ldiag[mRow-2]*m.udiag[mRow-3]+diag[mRow-2]*m.diag[mRow-2]+udiag[mRow-2]*m.ldiag[mRow-1];
595                                 array[mRow-2][mRow-1]=diag[mRow-2]*m.udiag[mRow-2]+udiag[mRow-2]*m.diag[mRow-1];
596                         }
597                         mRow--;
598                         array[mRow][mRow-2]=ldiag[mRow]*m.ldiag[mRow-1];
599                         array[mRow][mRow-1]=ldiag[mRow]*m.diag[mRow-1]+diag[mRow]*m.ldiag[mRow];
600                         array[mRow][mRow]=ldiag[mRow]*m.udiag[mRow-1]+diag[mRow]*m.diag[mRow];
601                         return new DoubleSquareMatrix(array);
602                 } else
603                         throw new MatrixDimensionException("Incompatible matrices.");
604         }
605         private DoubleSquareMatrix multiplyTridiagonal(final AbstractDoubleSquareMatrix m) {
606                 int mRow=numRows;
607                 if(numCols==m.rows()) {
608                         final double array[][]=new double[mRow][mRow];
609                         array[0][0]=diag[0]*m.getElement(0,0)+udiag[0]*m.getElement(1,0);
610                         array[0][1]=diag[0]*m.getElement(0,1)+udiag[0]*m.getElement(1,1);
611                         array[0][2]=udiag[0]*m.getElement(1,2);
612                         if(mRow>3) {
613                                 array[1][0]=ldiag[1]*m.getElement(0,0)+diag[1]*m.getElement(1,0);
614                                 array[1][1]=ldiag[1]*m.getElement(0,1)+diag[1]*m.getElement(1,1)+udiag[1]*m.getElement(2,1);
615                                 array[1][2]=diag[1]*m.getElement(1,2)+udiag[1]*m.getElement(2,2);
616                                 array[1][3]=udiag[1]*m.getElement(2,3);
617                         }
618                         if(mRow==3) {
619                                 array[1][0]=ldiag[1]*m.getElement(0,0)+diag[1]*m.getElement(1,0);
620                                 array[1][1]=ldiag[1]*m.getElement(0,1)+diag[1]*m.getElement(1,1)+udiag[1]*m.getElement(2,1);
621                                 array[1][2]=diag[1]*m.getElement(1,2)+udiag[1]*m.getElement(2,2);
622                         } else if(mRow>4) {
623                                 for(int i=2;i<mRow-2;i++) {
624                                         array[i][i-2]=ldiag[i]*m.getElement(i-1,i-2);
625                                         array[i][i-1]=ldiag[i]*m.getElement(i-1,i-1)+diag[i]*m.getElement(i,i-1);
626                                         array[i][i]=ldiag[i]*m.getElement(i-1,i)+diag[i]*m.getElement(i,i)+udiag[i]*m.getElement(i+1,i);
627                                         array[i][i+1]=diag[i]*m.getElement(i,i+1)+udiag[i]*m.getElement(i+1,i+1);
628                                         array[i][i+2]=udiag[i]*m.getElement(i+1,i+2);
629                                 }
630                         }
631                         if(mRow>3) {
632                                 array[mRow-2][mRow-4]=ldiag[mRow-2]*m.getElement(mRow-3,mRow-4);
633                                 array[mRow-2][mRow-3]=ldiag[mRow-2]*m.getElement(mRow-3,mRow-3)+diag[mRow-2]*m.getElement(mRow-2,mRow-3);
634                                 array[mRow-2][mRow-2]=ldiag[mRow-2]*m.getElement(mRow-3,mRow-2)+diag[mRow-2]*m.getElement(mRow-2,mRow-2)+udiag[mRow-2]*m.getElement(mRow-1,mRow-2);
635                                 array[mRow-2][mRow-1]=diag[mRow-2]*m.getElement(mRow-2,mRow-1)+udiag[mRow-2]*m.getElement(mRow-1,mRow-1);
636                         }
637                         mRow--;
638                         array[mRow][mRow-2]=ldiag[mRow]*m.getElement(mRow-1,mRow-2);
639                         array[mRow][mRow-1]=ldiag[mRow]*m.getElement(mRow-1,mRow-1)+diag[mRow]*m.getElement(mRow,mRow-1);
640                         array[mRow][mRow]=ldiag[mRow]*m.getElement(mRow-1,mRow)+diag[mRow]*m.getElement(mRow,mRow);
641                         return new DoubleSquareMatrix(array);
642                 } else {
643                         throw new MatrixDimensionException("Incompatible matrices.");
644                 }
645         }
646
647 // TRANSPOSE
648
649         /**
650         * Returns the transpose of this matrix.
651         * @return a double matrix
652         */
653         public Matrix transpose() {
654                 final DoubleTridiagonalMatrix ans=new DoubleTridiagonalMatrix(numRows);
655                 System.arraycopy(ldiag,1,ans.udiag,0,ldiag.length-1);
656                 System.arraycopy(diag,0,ans.diag,0,diag.length);
657                 System.arraycopy(udiag,0,ans.ldiag,1,udiag.length-1);
658                 return ans;
659         }
660
661 // CHOLESKY DECOMPOSITION
662
663         /**
664         * Returns the Cholesky decomposition of this matrix.
665         * Matrix must be symmetric and positive definite.
666         * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
667         */
668         public AbstractDoubleSquareMatrix[] choleskyDecompose() {
669                 final int N=numRows;
670                 final double arrayL[][]=new double[N][N];
671                 final double arrayU[][]=new double[N][N];
672                 double tmp=Math.sqrt(diag[0]);
673                 arrayL[0][0]=arrayU[0][0]=tmp;
674                 arrayL[1][0]=arrayU[0][1]=ldiag[1]/tmp;
675                 for(int j=1;j<N;j++) {
676                         tmp=diag[j];
677                         for(int i=0;i<j;i++)
678                                 tmp-=arrayL[j][i]*arrayL[j][i];
679                         arrayL[j][j]=arrayU[j][j]=Math.sqrt(tmp);
680                         if(j+1<N) {
681                                 tmp=ldiag[j+1];
682                                 for(int k=0;k<j+1;k++)
683                                         tmp-=arrayL[j][k]*arrayU[k][j+1];
684                                 arrayL[j+1][j]=arrayU[j][j+1]=tmp/arrayU[j][j];
685                                 for(int i=j+2;i<N;i++) {
686                                         tmp=0.0;
687                                         for(int k=0;k<i;k++)
688                                                 tmp-=arrayL[j][k]*arrayU[k][i];
689                                         arrayL[i][j]=arrayU[j][i]=tmp/arrayU[j][j];
690                                 }
691                         }
692                 }
693                 final AbstractDoubleSquareMatrix lu[]=new AbstractDoubleSquareMatrix[2];
694                 lu[0]=new DoubleSquareMatrix(arrayL);
695                 lu[1]=new DoubleSquareMatrix(arrayU);
696                 return lu;
697         }
698
699 // QR DECOMPOSITION
700
701         /**
702         * Returns the QR decomposition of this matrix.
703         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
704         * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
705         * @jsci.planetmath QRDecomposition
706         */
707         public AbstractDoubleSquareMatrix[] qrDecompose() {
708                 final int N=numRows;
709                 final double array[][]=new double[N][N];
710                 final double arrayQ[][]=new double[N][N];
711                 final double arrayR[][]=new double[N][N];
712                 // copy matrix
713 array[0][0]=diag[0];
714                 array[0][1]=udiag[0];
715                 for(int i=1;i<N-1;i++) {
716                         array[i][i-1]=ldiag[i];
717                         array[i][i]=diag[i];
718                         array[i][i+1]=udiag[i];
719                 }
720                 array[N-1][N-2]=ldiag[N-1];
721                 array[N-1][N-1]=diag[N-1];
722                 for(int k=0; k<N; k++) {
723                         // compute l2-norm of kth column
724 double norm = array[k][k];
725                         for(int i=k+1; i<N; i++)
726                                 norm = ExtraMath.hypot(norm, array[i][k]);
727                         if(norm != 0.0) {
728                                 // form kth Householder vector
729 if(array[k][k] < 0.0)
730                                         norm = -norm;
731                                 for(int i=k; i<N; i++)
732                                         array[i][k] /= norm;
733                                 array[k][k] += 1.0;
734                                 // apply transformation to remaining columns
735 for(int j=k+1; j<N; j++) {
736                                         double s=array[k][k]*array[k][j];
737                                         for(int i=k+1; i<N; i++)
738                                                 s += array[i][k]*array[i][j];
739                                         s /= array[k][k];
740                                         for(int i=k; i<N; i++)
741                                                 array[i][j] -= s*array[i][k];
742                                 }
743                         }
744                         arrayR[k][k] = -norm;
745                 }
746                 for(int k=N-1; k>=0; k--) {
747                         arrayQ[k][k] = 1.0;
748                         for(int j=k; j<N; j++) {
749                                 if(array[k][k] != 0.0) {
750                                         double s = array[k][k]*arrayQ[k][j];
751                                         for(int i=k+1; i<N; i++)
752                                                 s += array[i][k]*arrayQ[i][j];
753                                         s /= array[k][k];
754                                         for(int i=k; i<N; i++)
755                                                 arrayQ[i][j] -= s*array[i][k];
756                                 }
757                         }
758                 }
759                 for(int i=0; i<N; i++) {
760                         for(int j=i+1; j<N; j++)
761                                 arrayR[i][j] = array[i][j];
762                 }
763                 final AbstractDoubleSquareMatrix qr[]=new AbstractDoubleSquareMatrix[2];
764                 qr[0]=new DoubleSquareMatrix(arrayQ);
765                 qr[1]=new DoubleSquareMatrix(arrayR);
766                 return qr;
767         }
768
769 // SINGULAR VALUE DECOMPOSITION
770
771         /**
772         * Returns the singular value decomposition of this matrix.
773         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
774         * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
775         */
776         public AbstractDoubleSquareMatrix[] singularValueDecompose() {
777                 final int N=numRows;
778                 final int Nm1=N-1;
779                 final double array[][]=new double[N][N];
780                 final double arrayU[][]=new double[N][N];
781                 final double arrayS[]=new double[N];
782                 final double arrayV[][]=new double[N][N];
783                 final double e[]=new double[N];
784                 final double work[]=new double[N];
785                 // copy matrix
786 array[0][0]=diag[0];
787                 array[0][1]=udiag[0];
788                 for(int i=1;i<Nm1;i++) {
789                         array[i][i-1]=ldiag[i];
790                         array[i][i]=diag[i];
791                         array[i][i+1]=udiag[i];
792                 }
793                 array[Nm1][Nm1-1]=ldiag[Nm1];
794                 array[Nm1][Nm1]=diag[Nm1];
795                 // reduce matrix to bidiagonal form
796 for(int k=0;k<Nm1;k++) {
797                         // compute the transformation for the kth column
798 // compute l2-norm of kth column
799 arrayS[k]=array[k][k];
800                         for(int i=k+1;i<N;i++)
801                                 arrayS[k]=ExtraMath.hypot(arrayS[k],array[i][k]);
802                         if(arrayS[k]!=0.0) {
803                                 if(array[k][k]<0.0)
804                                         arrayS[k]=-arrayS[k];
805                                 for(int i=k;i<N;i++)
806                                         array[i][k]/=arrayS[k];
807                                 array[k][k]+=1.0;
808                         }
809                         arrayS[k]=-arrayS[k];
810                         for(int j=k+1;j<N;j++) {
811                                 if(arrayS[k]!=0.0) {
812                                         // apply the transformation
813 double t=array[k][k]*array[k][j];
814                                         for(int i=k+1; i<N; i++)
815                                                 t+=array[i][k]*array[i][j];
816                                         t /= array[k][k];
817                                         for(int i=k; i<N; i++)
818                                                 array[i][j] -= t*array[i][k];
819                                 }
820                                 e[j]=array[k][j];
821                         }
822                         for(int i=k;i<N;i++)
823                                 arrayU[i][k]=array[i][k];
824                         if(k<N-2) {
825                                 // compute the kth row transformation
826 // compute l2-norm of kth column
827 e[k]=e[k+1];
828                                 for(int i=k+2;i<N;i++)
829                                         e[k]=ExtraMath.hypot(e[k],e[i]);
830                                 if(e[k]!=0.0) {
831                                         if(e[k+1]<0.0)
832                                                 e[k]=-e[k];
833                                         for(int i=k+1;i<N;i++)
834                                                 e[i]/=e[k];
835                                         e[k+1]+=1.0;
836                                 }
837                                 e[k]=-e[k];
838                                 if(e[k]!=0.0) {
839                                         // apply the transformation
840 for(int i=k+1;i<N;i++) {
841                                                 work[i]=0.0;
842                                                 for(int j=k+1;j<N;j++)
843                                                         work[i]+=e[j]*array[i][j];
844                                         }
845                                         for(int j=k+1;j<N;j++) {
846                                                 double t = e[j]/e[k+1];
847                                                 for(int i=k+1;i<N;i++)
848                                                         array[i][j] -= t*work[i];
849                                         }
850                                 }
851                                 for(int i=k+1;i<N;i++)
852                                         arrayV[i][k]=e[i];
853                         }
854                 }
855                 // setup the final bidiagonal matrix of order p
856 int p=N;
857                 arrayS[Nm1]=array[Nm1][Nm1];
858                 e[N-2]=array[N-2][Nm1];
859                 e[Nm1]=0.0;
860                 arrayU[Nm1][Nm1]=1.0;
861                 for(int k=N-2;k>=0;k--) {
862                         if(arrayS[k]!=0.0) {
863                                 for(int j=k+1;j<N;j++) {
864                                         double t=arrayU[k][k]*arrayU[k][j];
865                                         for(int i=k+1;i<N;i++)
866                                                 t+=arrayU[i][k]*arrayU[i][j];
867                                         t /= arrayU[k][k];
868                                         for(int i=k;i<N;i++)
869                                                 arrayU[i][j] -= t*arrayU[i][k];
870                                 }
871                                 for(int i=k;i<N;i++)
872                                         arrayU[i][k]=-arrayU[i][k];
873                                 arrayU[k][k]+=1.0;
874                                 for(int i=0;i<k-1;i++)
875                                         arrayU[i][k]=0.0;
876                         } else {
877                                 for(int i=0;i<N;i++)
878                                         arrayU[i][k]=0.0;
879                                 arrayU[k][k]=1.0;
880                         }
881                 }
882                 for(int k=Nm1;k>=0;k--) {
883                         if(k<N-2 && e[k]!=0.0) {
884                                 for(int j=k+1;j<N;j++) {
885                                         double t=arrayV[k+1][k]*arrayV[k+1][j];
886                                         for(int i=k+2;i<N;i++)
887                                                 t+=arrayV[i][k]*arrayV[i][j];
888                                         t /= arrayV[k+1][k];
889                                         for(int i=k+1;i<N;i++)
890                                                 arrayV[i][j] -= t*arrayV[i][k];
891                                 }
892                         }
893                         for(int i=0;i<N;i++)
894                                 arrayV[i][k]=0.0;
895                         arrayV[k][k]=1.0;
896                 }
897                 final double eps=Math.pow(2.0,-52.0);
898                 int iter=0;
899                 while(p>0) {
900                         int k, action;
901                         // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
902 // action = 2 if arrayS[k] is negligible and k<p
903 // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
904 // action = 4 if e[p-1] is negligible (convergence)
905 for(k=p-2;k>=-1;k--) {
906                                 if(k==-1)
907                                         break;
908                                 if(Math.abs(e[k])<=eps*(Math.abs(arrayS[k])+Math.abs(arrayS[k+1]))) {
909                                         e[k]=0.0;
910                                         break;
911                                 }
912                         }
913                         if(k==p-2) {
914                                 action=4;
915                         } else {
916                                 int ks;
917                                 for(ks=p-1;ks>=k;ks--) {
918                                         if(ks==k)
919                                                 break;
920                                         double t=(ks!=p ? Math.abs(e[ks]) : 0.0)+(ks!=k+1 ? Math.abs(e[ks-1]) : 0.0);
921                                         if(Math.abs(arrayS[ks])<=eps*t) {
922                                                 arrayS[ks]=0.0;
923                                                 break;
924                                         }
925                                 }
926                                 if(ks==k) {
927                                         action=3;
928                                 } else if(ks==p-1) {
929                                         action=1;
930                                 } else {
931                                         action=2;
932                                         k=ks;
933                                 }
934                         }
935                         k++;
936                         switch(action) {
937                                 // deflate negligible arrayS[p]
938 case 1: {
939                                         double f=e[p-2];
940                                         e[p-2]=0.0;
941                                         for(int j=p-2;j>=k;j--) {
942                                                 double t=ExtraMath.hypot(arrayS[j],f);
943                                                 final double cs=arrayS[j]/t;
944                                                 final double sn=f/t;
945                                                 arrayS[j]=t;
946                                                 if(j!=k) {
947                                                         f=-sn*e[j-1];
948                                                         e[j-1]*=cs;
949                                                 }
950                                                 for(int i=0;i<N;i++) {
951                                                         t=cs*arrayV[i][j]+sn*arrayV[i][p-1];
952                                                         arrayV[i][p-1]=-sn*arrayV[i][j]+cs*arrayV[i][p-1];
953                                                         arrayV[i][j]=t;
954                                                 }
955                                         }
956                                         } break;
957                                 // split at negligible arrayS[k]
958 case 2: {
959                                         double f=e[k-1];
960                                         e[k-1]=0.0;
961                                         for(int j=k;j<p;j++) {
962                                                 double t=ExtraMath.hypot(arrayS[j],f);
963                                                 final double cs=arrayS[j]/t;
964                                                 final double sn=f/t;
965                                                 arrayS[j]=t;
966                                                 f=-sn*e[j];
967                                                 e[j]*=cs;
968                                                 for(int i=0;i<N;i++) {
969                                                         t=cs*arrayU[i][j]+sn*arrayU[i][k-1];
970                                                         arrayU[i][k-1]=-sn*arrayU[i][j]+cs*arrayU[i][k-1];
971                                                         arrayU[i][j]=t;
972                                                 }
973                                         }
974                                         } break;
975                                 // perform one QR step
976 case 3: {
977                                         // calculate the shift
978 final double scale=Math.max(Math.max(Math.max(Math.max(
979                                                 Math.abs(arrayS[p-1]),Math.abs(arrayS[p-2])),Math.abs(e[p-2])),
980                                                 Math.abs(arrayS[k])),Math.abs(e[k]));
981                                         double sp=arrayS[p-1]/scale;
982                                         double spm1=arrayS[p-2]/scale;
983                                         double epm1=e[p-2]/scale;
984                                         double sk=arrayS[k]/scale;
985                                         double ek=e[k]/scale;
986                                         double b=((spm1+sp)*(spm1-sp)+epm1*epm1)/2.0;
987                                         double c=(sp*epm1)*(sp*epm1);
988                                         double shift=0.0;
989                                         if(b!=0.0 || c!=0.0) {
990                                                 shift=Math.sqrt(b*b+c);
991                                                 if(b<0.0)
992                                                         shift=-shift;
993                                                 shift=c/(b+shift);
994                                         }
995                                         double f=(sk+sp)*(sk-sp)+shift;
996                                         double g=sk*ek;
997                                         // chase zeros
998 for(int j=k;j<p-1;j++) {
999                                                 double t=ExtraMath.hypot(f,g);
1000                                                double cs=f/t;
1001                                                double sn=g/t;
1002                                                if(j!=k)
1003                                                        e[j-1]=t;
1004                                                f=cs*arrayS[j]+sn*e[j];
1005                                                e[j]=cs*e[j]-sn*arrayS[j];
1006                                                g=sn*arrayS[j+1];
1007                                                arrayS[j+1]*=cs;
1008                                                for(int i=0;i<N;i++) {
1009                                                        t=cs*arrayV[i][j]+sn*arrayV[i][j+1];
1010                                                        arrayV[i][j+1]=-sn*arrayV[i][j]+cs*arrayV[i][j+1];
1011                                                        arrayV[i][j]=t;
1012                                                }
1013                                                t=ExtraMath.hypot(f,g);
1014                                                cs=f/t;
1015                                                sn=g/t;
1016                                                arrayS[j]=t;
1017                                                f=cs*e[j]+sn*arrayS[j+1];
1018                                                arrayS[j+1]=-sn*e[j]+cs*arrayS[j+1];
1019                                                g=sn*e[j+1];
1020                                                e[j+1]*=cs;
1021                                                if(j<Nm1) {
1022                                                        for(int i=0;i<N;i++) {
1023                                                                t=cs*arrayU[i][j]+sn*arrayU[i][j+1];
1024                                                                arrayU[i][j+1]=-sn*arrayU[i][j]+cs*arrayU[i][j+1];
1025                                                                arrayU[i][j]=t;
1026                                                        }
1027                                                }
1028                                        }
1029                                        e[p-2]=f;
1030                                        iter++;
1031                                        } break;
1032                                // convergence
1033 case 4: {
1034                                        // make the singular values positive
1035 if(arrayS[k]<=0.0) {
1036                                                arrayS[k]=-arrayS[k];
1037                                                for(int i=0;i<p;i++)
1038                                                        arrayV[i][k]=-arrayV[i][k];
1039                                        }
1040                                        // order the singular values
1041 while(k<p-1) {
1042                                                if(arrayS[k]>=arrayS[k+1])
1043                                                        break;
1044                                                double tmp=arrayS[k];
1045                                                arrayS[k]=arrayS[k+1];
1046                                                arrayS[k+1]=tmp;
1047                                                if(k<Nm1) {
1048                                                        for(int i=0;i<N;i++) {
1049                                                                tmp=arrayU[i][k+1];
1050                                                                arrayU[i][k+1]=arrayU[i][k];
1051                                                                arrayU[i][k]=tmp;
1052                                                                tmp=arrayV[i][k+1];
1053                                                                arrayV[i][k+1]=arrayV[i][k];
1054                                                                arrayV[i][k]=tmp;
1055                                                        }
1056                                                }
1057                                                k++;
1058                                        }
1059                                        iter=0;
1060                                        p--;
1061                                        } break;
1062                        }
1063                }
1064                final AbstractDoubleSquareMatrix svd[]=new AbstractDoubleSquareMatrix[3];
1065                svd[0]=new DoubleSquareMatrix(arrayU);
1066                svd[1]=new DoubleDiagonalMatrix(arrayS);
1067                svd[2]=new DoubleSquareMatrix(arrayV);
1068                return svd;
1069        }
1070
1071// MAP ELEMENTS
1072
1073        /**
1074        * Applies a function on all the matrix elements.
1075        * @param f a user-defined function
1076        * @return a double matrix
1077        */
1078        public AbstractDoubleMatrix mapElements(final Mapping f) {
1079        double zeroValue = f.map(0.0);
1080        if(Math.abs(zeroValue) <= JSci.GlobalSettings.ZERO_TOL)
1081            return tridiagonalMap(f);
1082        else
1083            return generalMap(f, zeroValue);
1084    }
1085    private AbstractDoubleMatrix tridiagonalMap(Mapping f) {
1086                int mRow=numRows;
1087                final DoubleTridiagonalMatrix ans=new DoubleTridiagonalMatrix(mRow);
1088                ans.diag[0]=f.map(diag[0]);
1089                ans.udiag[0]=f.map(udiag[0]);
1090                mRow--;
1091                for(int i=1;i<mRow;i++) {
1092                        ans.ldiag[i]=f.map(ldiag[i]);
1093                        ans.diag[i]=f.map(diag[i]);
1094                        ans.udiag[i]=f.map(udiag[i]);
1095                }
1096                ans.ldiag[mRow]=f.map(ldiag[mRow]);
1097                ans.diag[mRow]=f.map(diag[mRow]);
1098                return ans;
1099        }
1100    private AbstractDoubleMatrix generalMap(Mapping f, double zeroValue) {
1101                final double array[][]=new double[numRows][numRows];
1102        for(int i=0; i<numRows; i++) {
1103            for(int j=0; j<numRows; j++) {
1104                array[i][j] = zeroValue;
1105            }
1106        }
1107                int mRow=numRows;
1108                array[0][0]=f.map(diag[0]);
1109                array[0][1]=f.map(udiag[0]);
1110                mRow--;
1111                for(int i=1;i<mRow;i++) {
1112                        array[i][i-1]=f.map(ldiag[i]);
1113                        array[i][i]=f.map(diag[i]);
1114                        array[i][i+1]=f.map(udiag[i]);
1115                }
1116                array[mRow][mRow-1]=f.map(ldiag[mRow]);
1117                array[mRow][mRow]=f.map(diag[mRow]);
1118                return new DoubleSquareMatrix(array);
1119        }
1120}
1121
Popular Tags