KickJava   Java API By Example, From Geeks To Geeks.

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


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.AbstractIntegerVector;
11 import JSci.maths.vectors.IntegerVector;
12 import JSci.maths.groups.AbelianGroup;
13 import JSci.maths.algebras.*;
14 import JSci.maths.fields.*;
15
16 /**
17 * The IntegerTridiagonalMatrix class provides an object for encapsulating integer tridiagonal matrices.
18 * @version 2.2
19 * @author Mark Hale
20 */
21 public class IntegerTridiagonalMatrix extends AbstractIntegerSquareMatrix implements TridiagonalMatrix {
22         /**
23         * Tridiagonal data.
24         */
25         protected final int ldiag[];
26         protected final int diag[];
27         protected final int udiag[];
28         /**
29         * Constructs an empty matrix.
30         * @param size the number of rows/columns
31         */
32         public IntegerTridiagonalMatrix(final int size) {
33                 super(size);
34                 ldiag=new int[size];
35                 diag=new int[size];
36                 udiag=new int[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 IntegerTridiagonalMatrix(final int 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 int matrix
62         */
63         public boolean equals(AbstractIntegerMatrix m, double tol) {
64                 if(m instanceof TridiagonalMatrix) {
65                         if(numRows != m.rows() || numCols != m.columns())
66                                 return false;
67             int sumSqr = 0;
68             int ldelta;
69             int delta = diag[0] - m.getElement(0,0);
70             int 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 a double matrix.
103         * @return a double matrix
104         */
105         public AbstractDoubleMatrix toDoubleMatrix() {
106                 final DoubleTridiagonalMatrix m=new DoubleTridiagonalMatrix(numRows);
107                 m.diag[0]=diag[0];
108                 m.udiag[0]=udiag[0];
109                 int i=1;
110                 for(;i<numRows-1;i++) {
111                         m.ldiag[i]=ldiag[i];
112                         m.diag[i]=diag[i];
113                         m.udiag[i]=udiag[i];
114                 }
115                 m.ldiag[i]=ldiag[i];
116                 m.diag[i]=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 int 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 int 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 int trace() {
193                 int 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 int infNorm() {
203                 int result=Math.abs(diag[0])+Math.abs(udiag[0]);
204                 int 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 int matrix
245         * @exception MatrixDimensionException If the matrices are different sizes.
246         */
247         public AbstractIntegerSquareMatrix add(final AbstractIntegerSquareMatrix m) {
248                 if(m instanceof IntegerTridiagonalMatrix)
249                         return add((IntegerTridiagonalMatrix)m);
250                 if(m instanceof TridiagonalMatrix)
251                         return addTridiagonal(m);
252                 if(m instanceof IntegerSquareMatrix)
253                         return add((IntegerSquareMatrix)m);
254
255                 if(numRows==m.rows() && numCols==m.columns()) {
256                         final int array[][]=new int[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 IntegerSquareMatrix(array);
263                 } else {
264                         throw new MatrixDimensionException("Matrices are different sizes.");
265                 }
266         }
267         public IntegerSquareMatrix add(final IntegerSquareMatrix m) {
268                 if(numRows==m.numRows && numCols==m.numCols) {
269                         final int array[][]=new int[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 IntegerSquareMatrix(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 int tridiagonal matrix
289         * @exception MatrixDimensionException If the matrices are different sizes.
290         */
291         public IntegerTridiagonalMatrix add(final IntegerTridiagonalMatrix m) {
292                 int mRow=numRows;
293                 if(mRow==m.numRows) {
294                         final IntegerTridiagonalMatrix ans=new IntegerTridiagonalMatrix(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 IntegerTridiagonalMatrix addTridiagonal(final AbstractIntegerSquareMatrix m) {
310                 int mRow=numRows;
311                 if(mRow==m.rows()) {
312                         final IntegerTridiagonalMatrix ans=new IntegerTridiagonalMatrix(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 int matrix
334         * @exception MatrixDimensionException If the matrices are different sizes.
335         */
336         public AbstractIntegerSquareMatrix subtract(final AbstractIntegerSquareMatrix m) {
337                 if(m instanceof IntegerTridiagonalMatrix)
338                         return subtract((IntegerTridiagonalMatrix)m);
339                 if(m instanceof TridiagonalMatrix)
340                         return subtractTridiagonal(m);
341                 if(m instanceof IntegerSquareMatrix)
342                         return subtract((IntegerSquareMatrix)m);
343
344                 if(numRows==m.rows() && numCols==m.columns()) {
345                         final int array[][]=new int[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 IntegerSquareMatrix(array);
352                 } else {
353                         throw new MatrixDimensionException("Matrices are different sizes.");
354                 }
355         }
356         public IntegerSquareMatrix subtract(final IntegerSquareMatrix m) {
357                 if(numRows==m.numRows && numCols==m.numCols) {
358                         final int array[][]=new int[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 IntegerSquareMatrix(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 int tridiagonal matrix
371         * @exception MatrixDimensionException If the matrices are different sizes.
372         */
373         public IntegerTridiagonalMatrix subtract(final IntegerTridiagonalMatrix m) {
374                 int mRow=numRows;
375                 if(mRow==m.numRows) {
376                         final IntegerTridiagonalMatrix ans=new IntegerTridiagonalMatrix(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 IntegerTridiagonalMatrix subtractTridiagonal(final AbstractIntegerSquareMatrix m) {
392                 int mRow=numRows;
393                 if(mRow==m.rows()) {
394                         final IntegerTridiagonalMatrix ans=new IntegerTridiagonalMatrix(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 int.
416         * @return a int square matrix.
417         */
418         public AbstractIntegerMatrix scalarMultiply(final int x) {
419                 final IntegerTridiagonalMatrix ans=new IntegerTridiagonalMatrix(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 // SCALAR PRODUCT
437
438         /**
439         * Returns the scalar product of this matrix and another.
440         * @param m a int matrix.
441         * @exception MatrixDimensionException If the matrices are different sizes.
442         */
443         public int scalarProduct(final AbstractIntegerSquareMatrix m) {
444                 if(m instanceof IntegerTridiagonalMatrix)
445                         return scalarProduct((IntegerTridiagonalMatrix)m);
446                 if(m instanceof IntegerSquareMatrix)
447                         return scalarProduct((IntegerSquareMatrix)m);
448
449                 if(numRows==m.rows() && numCols==m.columns()) {
450                         int ans = diag[0]*m.getElement(0,0)+udiag[0]*m.getElement(0,1);
451                         final int lastRow = numRows-1;
452                         for(int i=1;i<lastRow;i++)
453                                 ans += ldiag[i]*m.getElement(i,i-1)+diag[i]*m.getElement(i,i)+udiag[i]*m.getElement(i,i+1);
454                         ans += ldiag[lastRow]*m.getElement(lastRow,lastRow-1)+diag[lastRow]*m.getElement(lastRow,lastRow);
455                         return ans;
456                 } else {
457                        throw new MatrixDimensionException("Matrices are different sizes.");
458                 }
459         }
460         public int scalarProduct(final IntegerSquareMatrix m) {
461                 if(numRows==m.numRows && numCols==m.numCols) {
462                         int ans = diag[0]*m.matrix[0][0]+udiag[0]*m.matrix[0][1];
463                         final int lastRow = numRows-1;
464                         for(int i=1;i<lastRow;i++)
465                                 ans += ldiag[i]*m.matrix[i][i-1]+diag[i]*m.matrix[i][i]+udiag[i]*m.matrix[i][i+1];
466                         ans += ldiag[lastRow]*m.matrix[lastRow][lastRow-1]+diag[lastRow]*m.matrix[lastRow][lastRow];
467                         return ans;
468                 } else
469                         throw new MatrixDimensionException("Matrices are different sizes.");
470         }
471         public int scalarProduct(final IntegerTridiagonalMatrix m) {
472                 if(numRows==m.numRows) {
473                         int ans = diag[0]*m.diag[0]+udiag[0]*m.udiag[0];
474                         final int lastRow = numRows-1;
475                         for(int i=1;i<lastRow;i++)
476                                 ans += ldiag[i]*m.ldiag[i]+diag[i]*m.diag[i]+udiag[i]*m.udiag[i];
477                         ans += ldiag[lastRow]*m.ldiag[lastRow]+diag[lastRow]*m.diag[lastRow];
478                         return ans;
479                 } else
480                         throw new MatrixDimensionException("Matrices are different sizes.");
481         }
482
483 // MATRIX MULTIPLICATION
484
485         /**
486         * Returns the multiplication of a vector by this matrix.
487         * @param v a int vector.
488         * @exception DimensionException If the matrix and vector are incompatible.
489         */
490         public AbstractIntegerVector multiply(final AbstractIntegerVector v) {
491                 if(numCols==v.dimension()) {
492                         final int array[]=new int[numRows];
493                         final int lastRow = numRows-1;
494                         array[0]=diag[0]*v.getComponent(0)+udiag[0]*v.getComponent(1);
495                         for(int i=1;i<lastRow;i++)
496                                 array[i]=ldiag[i]*v.getComponent(i-1)+diag[i]*v.getComponent(i)+udiag[i]*v.getComponent(i+1);
497                         array[lastRow]=ldiag[lastRow]*v.getComponent(lastRow-1)+diag[lastRow]*v.getComponent(lastRow);
498                         return new IntegerVector(array);
499                 } else {
500                         throw new DimensionException("Matrix and vector are incompatible.");
501                 }
502         }
503         /**
504         * Returns the multiplication of this matrix and another.
505         * @param m a int matrix
506         * @return a AbstractIntegerMatrix or a AbstractIntegerSquareMatrix as appropriate
507         * @exception MatrixDimensionException If the matrices are incompatible.
508         */
509         public AbstractIntegerSquareMatrix multiply(final AbstractIntegerSquareMatrix m) {
510                 if(m instanceof IntegerTridiagonalMatrix)
511                         return multiply((IntegerTridiagonalMatrix)m);
512                 if(m instanceof TridiagonalMatrix)
513                         return multiplyTridiagonal(m);
514                 if(m instanceof IntegerSquareMatrix)
515                         return multiply((IntegerSquareMatrix)m);
516
517                 if(numCols==m.rows()) {
518                         final int mColumns = m.columns();
519                         final int array[][]=new int[numRows][mColumns];
520                         final int lastRow = numRows-1;
521                         for(int j=0; j<numRows; j++) {
522                                 array[0][j]=diag[0]*m.getElement(0,j)+udiag[0]*m.getElement(1,j);
523                                 for(int i=1;i<lastRow;i++)
524                                         array[i][j]=ldiag[i]*m.getElement(i-1,j)+diag[i]*m.getElement(i,j)+udiag[i]*m.getElement(i+1,j);
525                                 array[lastRow][j]=ldiag[lastRow]*m.getElement(lastRow-1,j)+diag[lastRow]*m.getElement(lastRow,j);
526                         }
527                         return new IntegerSquareMatrix(array);
528                 } else {
529                         throw new MatrixDimensionException("Incompatible matrices.");
530                 }
531         }
532         public IntegerSquareMatrix multiply(final IntegerSquareMatrix m) {
533                 if(numCols==m.numRows) {
534                         final int array[][]=new int[numRows][m.numCols];
535                         final int lastRow = numRows-1;
536                         for(int j=0;j<numRows;j++) {
537                                 array[0][j]=diag[0]*m.matrix[0][j]+udiag[0]*m.matrix[1][j];
538                                 for(int i=1;i<lastRow;i++)
539                                         array[i][j]=ldiag[i]*m.matrix[i-1][j]+diag[i]*m.matrix[i][j]+udiag[i]*m.matrix[i+1][j];
540                                 array[lastRow][j]=ldiag[lastRow]*m.matrix[lastRow-1][j]+diag[lastRow]*m.matrix[lastRow][j];
541                         }
542                         return new IntegerSquareMatrix(array);
543                 } else
544                         throw new MatrixDimensionException("Incompatible matrices.");
545         }
546         public IntegerSquareMatrix multiply(final IntegerTridiagonalMatrix m) {
547                 int mRow=numRows;
548                 if(numCols==m.numRows) {
549                         final int array[][]=new int[mRow][mRow];
550                         array[0][0]=diag[0]*m.diag[0]+udiag[0]*m.ldiag[1];
551                         array[0][1]=diag[0]*m.udiag[0]+udiag[0]*m.diag[1];
552                         array[0][2]=udiag[0]*m.udiag[1];
553                         if(mRow>3) {
554                                 array[1][0]=ldiag[1]*m.diag[0]+diag[1]*m.ldiag[1];
555                                 array[1][1]=ldiag[1]*m.udiag[0]+diag[1]*m.diag[1]+udiag[1]*m.ldiag[2];
556                                 array[1][2]=diag[1]*m.udiag[1]+udiag[1]*m.diag[2];
557                                 array[1][3]=udiag[1]*m.udiag[2];
558                         }
559                         if(mRow==3) {
560                                 array[1][0]=ldiag[1]*m.diag[0]+diag[1]*m.ldiag[1];
561                                 array[1][1]=ldiag[1]*m.udiag[0]+diag[1]*m.diag[1]+udiag[1]*m.ldiag[2];
562                                 array[1][2]=diag[1]*m.udiag[1]+udiag[1]*m.diag[2];
563                         } else if(mRow>4) {
564                                 for(int i=2;i<mRow-2;i++) {
565                                         array[i][i-2]=ldiag[i]*m.ldiag[i-1];
566                                         array[i][i-1]=ldiag[i]*m.diag[i-1]+diag[i]*m.ldiag[i];
567                                         array[i][i]=ldiag[i]*m.udiag[i-1]+diag[i]*m.diag[i]+udiag[i]*m.ldiag[i+1];
568                                         array[i][i+1]=diag[i]*m.udiag[i]+udiag[i]*m.diag[i+1];
569                                         array[i][i+2]=udiag[i]*m.udiag[i+1];
570                                 }
571                         }
572                         if(mRow>3) {
573                                 array[mRow-2][mRow-4]=ldiag[mRow-2]*m.ldiag[mRow-3];
574                                 array[mRow-2][mRow-3]=ldiag[mRow-2]*m.diag[mRow-3]+diag[mRow-2]*m.ldiag[mRow-2];
575                                 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];
576                                 array[mRow-2][mRow-1]=diag[mRow-2]*m.udiag[mRow-2]+udiag[mRow-2]*m.diag[mRow-1];
577                         }
578                         mRow--;
579                         array[mRow][mRow-2]=ldiag[mRow]*m.ldiag[mRow-1];
580                         array[mRow][mRow-1]=ldiag[mRow]*m.diag[mRow-1]+diag[mRow]*m.ldiag[mRow];
581                         array[mRow][mRow]=ldiag[mRow]*m.udiag[mRow-1]+diag[mRow]*m.diag[mRow];
582                         return new IntegerSquareMatrix(array);
583                 } else
584                         throw new MatrixDimensionException("Incompatible matrices.");
585         }
586         private IntegerSquareMatrix multiplyTridiagonal(final AbstractIntegerSquareMatrix m) {
587                 int mRow=numRows;
588                 if(numCols==m.rows()) {
589                         final int array[][]=new int[mRow][mRow];
590                         array[0][0]=diag[0]*m.getElement(0,0)+udiag[0]*m.getElement(1,0);
591                         array[0][1]=diag[0]*m.getElement(0,1)+udiag[0]*m.getElement(1,1);
592                         array[0][2]=udiag[0]*m.getElement(1,2);
593                         if(mRow>3) {
594                                 array[1][0]=ldiag[1]*m.getElement(0,0)+diag[1]*m.getElement(1,0);
595                                 array[1][1]=ldiag[1]*m.getElement(0,1)+diag[1]*m.getElement(1,1)+udiag[1]*m.getElement(2,1);
596                                 array[1][2]=diag[1]*m.getElement(1,2)+udiag[1]*m.getElement(2,2);
597                                 array[1][3]=udiag[1]*m.getElement(2,3);
598                         }
599                         if(mRow==3) {
600                                 array[1][0]=ldiag[1]*m.getElement(0,0)+diag[1]*m.getElement(1,0);
601                                 array[1][1]=ldiag[1]*m.getElement(0,1)+diag[1]*m.getElement(1,1)+udiag[1]*m.getElement(2,1);
602                                 array[1][2]=diag[1]*m.getElement(1,2)+udiag[1]*m.getElement(2,2);
603                         } else if(mRow>4) {
604                                 for(int i=2;i<mRow-2;i++) {
605                                         array[i][i-2]=ldiag[i]*m.getElement(i-1,i-2);
606                                         array[i][i-1]=ldiag[i]*m.getElement(i-1,i-1)+diag[i]*m.getElement(i,i-1);
607                                         array[i][i]=ldiag[i]*m.getElement(i-1,i)+diag[i]*m.getElement(i,i)+udiag[i]*m.getElement(i+1,i);
608                                         array[i][i+1]=diag[i]*m.getElement(i,i+1)+udiag[i]*m.getElement(i+1,i+1);
609                                         array[i][i+2]=udiag[i]*m.getElement(i+1,i+2);
610                                 }
611                         }
612                         if(mRow>3) {
613                                 array[mRow-2][mRow-4]=ldiag[mRow-2]*m.getElement(mRow-3,mRow-4);
614                                 array[mRow-2][mRow-3]=ldiag[mRow-2]*m.getElement(mRow-3,mRow-3)+diag[mRow-2]*m.getElement(mRow-2,mRow-3);
615                                 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);
616                                 array[mRow-2][mRow-1]=diag[mRow-2]*m.getElement(mRow-2,mRow-1)+udiag[mRow-2]*m.getElement(mRow-1,mRow-1);
617                         }
618                         mRow--;
619                         array[mRow][mRow-2]=ldiag[mRow]*m.getElement(mRow-1,mRow-2);
620                         array[mRow][mRow-1]=ldiag[mRow]*m.getElement(mRow-1,mRow-1)+diag[mRow]*m.getElement(mRow,mRow-1);
621                         array[mRow][mRow]=ldiag[mRow]*m.getElement(mRow-1,mRow)+diag[mRow]*m.getElement(mRow,mRow);
622                         return new IntegerSquareMatrix(array);
623                 } else {
624                         throw new MatrixDimensionException("Incompatible matrices.");
625                 }
626         }
627
628 // TRANSPOSE
629
630         /**
631         * Returns the transpose of this matrix.
632         * @return a int matrix
633         */
634         public Matrix transpose() {
635                 final IntegerTridiagonalMatrix ans=new IntegerTridiagonalMatrix(numRows);
636                 System.arraycopy(ldiag,1,ans.udiag,0,ldiag.length-1);
637                 System.arraycopy(diag,0,ans.diag,0,diag.length);
638                 System.arraycopy(udiag,0,ans.ldiag,1,udiag.length-1);
639                 return ans;
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 int N=numRows;
651                 final double arrayL[][]=new double[N][N];
652                 final double arrayU[][]=new double[N][N];
653                 double tmp=Math.sqrt(diag[0]);
654                 arrayL[0][0]=arrayU[0][0]=tmp;
655                 arrayL[1][0]=arrayU[0][1]=ldiag[1]/tmp;
656                 for(int j=1;j<N;j++) {
657                         tmp=diag[j];
658                         for(int i=0;i<j;i++)
659                                 tmp-=arrayL[j][i]*arrayL[j][i];
660                         arrayL[j][j]=arrayU[j][j]=Math.sqrt(tmp);
661                         if(j+1<N) {
662                                 tmp=ldiag[j+1];
663                                 for(int k=0;k<j+1;k++)
664                                         tmp-=arrayL[j][k]*arrayU[k][j+1];
665                                 arrayL[j+1][j]=arrayU[j][j+1]=tmp/arrayU[j][j];
666                                 for(int i=j+2;i<N;i++) {
667                                         tmp=0.0;
668                                         for(int k=0;k<i;k++)
669                                                 tmp-=arrayL[j][k]*arrayU[k][i];
670                                         arrayL[i][j]=arrayU[j][i]=tmp/arrayU[j][j];
671                                 }
672                         }
673                 }
674                 final AbstractDoubleSquareMatrix lu[]=new AbstractDoubleSquareMatrix[2];
675                 lu[0]=new DoubleSquareMatrix(arrayL);
676                 lu[1]=new DoubleSquareMatrix(arrayU);
677                 return lu;
678         }
679
680 // QR DECOMPOSITION
681
682         /**
683         * Returns the QR decomposition of this matrix.
684         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
685         * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
686         * @jsci.planetmath QRDecomposition
687         */
688         public AbstractDoubleSquareMatrix[] qrDecompose() {
689                 final int N=numRows;
690                 final double array[][]=new double[N][N];
691                 final double arrayQ[][]=new double[N][N];
692                 final double arrayR[][]=new double[N][N];
693                 // copy matrix
694 array[0][0]=diag[0];
695                 array[0][1]=udiag[0];
696                 for(int i=1;i<N-1;i++) {
697                         array[i][i-1]=ldiag[i];
698                         array[i][i]=diag[i];
699                         array[i][i+1]=udiag[i];
700                 }
701                 array[N-1][N-2]=ldiag[N-1];
702                 array[N-1][N-1]=diag[N-1];
703                 for(int k=0; k<N; k++) {
704                         // compute l2-norm of kth column
705 double norm = array[k][k];
706                         for(int i=k+1; i<N; i++)
707                                 norm = ExtraMath.hypot(norm, array[i][k]);
708                         if(norm != 0.0) {
709                                 // form kth Householder vector
710 if(array[k][k] < 0.0)
711                                         norm = -norm;
712                                 for(int i=k; i<N; i++)
713                                         array[i][k] /= norm;
714                                 array[k][k] += 1.0;
715                                 // apply transformation to remaining columns
716 for(int j=k+1; j<N; j++) {
717                                         double s=array[k][k]*array[k][j];
718                                         for(int i=k+1; i<N; i++)
719                                                 s += array[i][k]*array[i][j];
720                                         s /= array[k][k];
721                                         for(int i=k; i<N; i++)
722                                                 array[i][j] -= s*array[i][k];
723                                 }
724                         }
725                         arrayR[k][k] = -norm;
726                 }
727                 for(int k=N-1; k>=0; k--) {
728                         arrayQ[k][k] = 1.0;
729                         for(int j=k; j<N; j++) {
730                                 if(array[k][k] != 0.0) {
731                                         double s = array[k][k]*arrayQ[k][j];
732                                         for(int i=k+1; i<N; i++)
733                                                 s += array[i][k]*arrayQ[i][j];
734                                         s /= array[k][k];
735                                         for(int i=k; i<N; i++)
736                                                 arrayQ[i][j] -= s*array[i][k];
737                                 }
738                         }
739                 }
740                 for(int i=0; i<N; i++) {
741                         for(int j=i+1; j<N; j++)
742                                 arrayR[i][j] = array[i][j];
743                 }
744                 final AbstractDoubleSquareMatrix qr[]=new AbstractDoubleSquareMatrix[2];
745                 qr[0]=new DoubleSquareMatrix(arrayQ);
746                 qr[1]=new DoubleSquareMatrix(arrayR);
747                 return qr;
748         }
749
750 // SINGULAR VALUE DECOMPOSITION
751
752         /**
753         * Returns the singular value decomposition of this matrix.
754         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
755         * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
756         */
757         public AbstractDoubleSquareMatrix[] singularValueDecompose() {
758                 final int N=numRows;
759                 final int Nm1=N-1;
760                 final double array[][]=new double[N][N];
761                 final double arrayU[][]=new double[N][N];
762                 final double arrayS[]=new double[N];
763                 final double arrayV[][]=new double[N][N];
764                 final double e[]=new double[N];
765                 final double work[]=new double[N];
766                 // copy matrix
767 array[0][0]=diag[0];
768                 array[0][1]=udiag[0];
769                 for(int i=1;i<Nm1;i++) {
770                         array[i][i-1]=ldiag[i];
771                         array[i][i]=diag[i];
772                         array[i][i+1]=udiag[i];
773                 }
774                 array[Nm1][Nm1-1]=ldiag[Nm1];
775                 array[Nm1][Nm1]=diag[Nm1];
776                 // reduce matrix to bidiagonal form
777 for(int k=0;k<Nm1;k++) {
778                         // compute the transformation for the kth column
779 // compute l2-norm of kth column
780 arrayS[k]=array[k][k];
781                         for(int i=k+1;i<N;i++)
782                                 arrayS[k]=ExtraMath.hypot(arrayS[k],array[i][k]);
783                         if(arrayS[k]!=0.0) {
784                                 if(array[k][k]<0.0)
785                                         arrayS[k]=-arrayS[k];
786                                 for(int i=k;i<N;i++)
787                                         array[i][k]/=arrayS[k];
788                                 array[k][k]+=1.0;
789                         }
790                         arrayS[k]=-arrayS[k];
791                         for(int j=k+1;j<N;j++) {
792                                 if(arrayS[k]!=0.0) {
793                                         // apply the transformation
794 double t=array[k][k]*array[k][j];
795                                         for(int i=k+1; i<N; i++)
796                                                 t+=array[i][k]*array[i][j];
797                                         t /= array[k][k];
798                                         for(int i=k; i<N; i++)
799                                                 array[i][j] -= t*array[i][k];
800                                 }
801                                 e[j]=array[k][j];
802                         }
803                         for(int i=k;i<N;i++)
804                                 arrayU[i][k]=array[i][k];
805                         if(k<N-2) {
806                                 // compute the kth row transformation
807 // compute l2-norm of kth column
808 e[k]=e[k+1];
809                                 for(int i=k+2;i<N;i++)
810                                         e[k]=ExtraMath.hypot(e[k],e[i]);
811                                 if(e[k]!=0.0) {
812                                         if(e[k+1]<0.0)
813                                                 e[k]=-e[k];
814                                         for(int i=k+1;i<N;i++)
815                                                 e[i]/=e[k];
816                                         e[k+1]+=1.0;
817                                 }
818                                 e[k]=-e[k];
819                                 if(e[k]!=0.0) {
820                                         // apply the transformation
821 for(int i=k+1;i<N;i++) {
822                                                 work[i]=0.0;
823                                                 for(int j=k+1;j<N;j++)
824                                                         work[i]+=e[j]*array[i][j];
825                                         }
826                                         for(int j=k+1;j<N;j++) {
827                                                 double t = e[j]/e[k+1];
828                                                 for(int i=k+1;i<N;i++)
829                                                         array[i][j] -= t*work[i];
830                                         }
831                                 }
832                                 for(int i=k+1;i<N;i++)
833                                         arrayV[i][k]=e[i];
834                         }
835                 }
836                 // setup the final bidiagonal matrix of order p
837 int p=N;
838                 arrayS[Nm1]=array[Nm1][Nm1];
839                 e[N-2]=array[N-2][Nm1];
840                 e[Nm1]=0.0;
841                 arrayU[Nm1][Nm1]=1.0;
842                 for(int k=N-2;k>=0;k--) {
843                         if(arrayS[k]!=0.0) {
844                                 for(int j=k+1;j<N;j++) {
845                                         double t=arrayU[k][k]*arrayU[k][j];
846                                         for(int i=k+1;i<N;i++)
847                                                 t+=arrayU[i][k]*arrayU[i][j];
848                                         t /= arrayU[k][k];
849                                         for(int i=k;i<N;i++)
850                                                 arrayU[i][j] -= t*arrayU[i][k];
851                                 }
852                                 for(int i=k;i<N;i++)
853                                         arrayU[i][k]=-arrayU[i][k];
854                                 arrayU[k][k]+=1.0;
855                                 for(int i=0;i<k-1;i++)
856                                         arrayU[i][k]=0.0;
857                         } else {
858                                 for(int i=0;i<N;i++)
859                                         arrayU[i][k]=0.0;
860                                 arrayU[k][k]=1.0;
861                         }
862                 }
863                 for(int k=Nm1;k>=0;k--) {
864                         if(k<N-2 && e[k]!=0.0) {
865                                 for(int j=k+1;j<N;j++) {
866                                         double t=arrayV[k+1][k]*arrayV[k+1][j];
867                                         for(int i=k+2;i<N;i++)
868                                                 t+=arrayV[i][k]*arrayV[i][j];
869                                         t /= arrayV[k+1][k];
870                                         for(int i=k+1;i<N;i++)
871                                                 arrayV[i][j] -= t*arrayV[i][k];
872                                 }
873                         }
874                         for(int i=0;i<N;i++)
875                                 arrayV[i][k]=0.0;
876                         arrayV[k][k]=1.0;
877                 }
878                 final double eps=Math.pow(2.0,-52.0);
879                 int iter=0;
880                 while(p>0) {
881                         int k, action;
882                         // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
883 // action = 2 if arrayS[k] is negligible and k<p
884 // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
885 // action = 4 if e[p-1] is negligible (convergence)
886 for(k=p-2;k>=-1;k--) {
887                                 if(k==-1)
888                                         break;
889                                 if(Math.abs(e[k])<=eps*(Math.abs(arrayS[k])+Math.abs(arrayS[k+1]))) {
890                                         e[k]=0.0;
891                                         break;
892                                 }
893                         }
894                         if(k==p-2) {
895                                 action=4;
896                         } else {
897                                 int ks;
898                                 for(ks=p-1;ks>=k;ks--) {
899                                         if(ks==k)
900                                                 break;
901                                         double t=(ks!=p ? Math.abs(e[ks]) : 0.0)+(ks!=k+1 ? Math.abs(e[ks-1]) : 0.0);
902                                         if(Math.abs(arrayS[ks])<=eps*t) {
903                                                 arrayS[ks]=0.0;
904                                                 break;
905                                         }
906                                 }
907                                 if(ks==k) {
908                                         action=3;
909                                 } else if(ks==p-1) {
910                                         action=1;
911                                 } else {
912                                         action=2;
913                                         k=ks;
914                                 }
915                         }
916                         k++;
917                         switch(action) {
918                                 // deflate negligible arrayS[p]
919 case 1: {
920                                         double f=e[p-2];
921                                         e[p-2]=0.0;
922                                         for(int j=p-2;j>=k;j--) {
923                                                 double t=ExtraMath.hypot(arrayS[j],f);
924                                                 final double cs=arrayS[j]/t;
925                                                 final double sn=f/t;
926                                                 arrayS[j]=t;
927                                                 if(j!=k) {
928                                                         f=-sn*e[j-1];
929                                                         e[j-1]*=cs;
930                                                 }
931                                                 for(int i=0;i<N;i++) {
932                                                         t=cs*arrayV[i][j]+sn*arrayV[i][p-1];
933                                                         arrayV[i][p-1]=-sn*arrayV[i][j]+cs*arrayV[i][p-1];
934                                                         arrayV[i][j]=t;
935                                                 }
936                                         }
937                                         } break;
938                                 // split at negligible arrayS[k]
939 case 2: {
940                                         double f=e[k-1];
941                                         e[k-1]=0.0;
942                                         for(int j=k;j<p;j++) {
943                                                 double t=ExtraMath.hypot(arrayS[j],f);
944                                                 final double cs=arrayS[j]/t;
945                                                 final double sn=f/t;
946                                                 arrayS[j]=t;
947                                                 f=-sn*e[j];
948                                                 e[j]*=cs;
949                                                 for(int i=0;i<N;i++) {
950                                                         t=cs*arrayU[i][j]+sn*arrayU[i][k-1];
951                                                         arrayU[i][k-1]=-sn*arrayU[i][j]+cs*arrayU[i][k-1];
952                                                         arrayU[i][j]=t;
953                                                 }
954                                         }
955                                         } break;
956                                 // perform one QR step
957 case 3: {
958                                         // calculate the shift
959 final double scale=Math.max(Math.max(Math.max(Math.max(
960                                                 Math.abs(arrayS[p-1]),Math.abs(arrayS[p-2])),Math.abs(e[p-2])),
961                                                 Math.abs(arrayS[k])),Math.abs(e[k]));
962                                         double sp=arrayS[p-1]/scale;
963                                         double spm1=arrayS[p-2]/scale;
964                                         double epm1=e[p-2]/scale;
965                                         double sk=arrayS[k]/scale;
966                                         double ek=e[k]/scale;
967                                         double b=((spm1+sp)*(spm1-sp)+epm1*epm1)/2.0;
968                                         double c=(sp*epm1)*(sp*epm1);
969                                         double shift=0.0;
970                                         if(b!=0.0 || c!=0.0) {
971                                                 shift=Math.sqrt(b*b+c);
972                                                 if(b<0.0)
973                                                         shift=-shift;
974                                                 shift=c/(b+shift);
975                                         }
976                                         double f=(sk+sp)*(sk-sp)+shift;
977                                         double g=sk*ek;
978                                         // chase zeros
979 for(int j=k;j<p-1;j++) {
980                                                 double t=ExtraMath.hypot(f,g);
981                                                 double cs=f/t;
982                                                 double sn=g/t;
983                                                 if(j!=k)
984                                                         e[j-1]=t;
985                                                 f=cs*arrayS[j]+sn*e[j];
986                                                 e[j]=cs*e[j]-sn*arrayS[j];
987                                                 g=sn*arrayS[j+1];
988                                                 arrayS[j+1]*=cs;
989                                                 for(int i=0;i<N;i++) {
990                                                         t=cs*arrayV[i][j]+sn*arrayV[i][j+1];
991                                                         arrayV[i][j+1]=-sn*arrayV[i][j]+cs*arrayV[i][j+1];
992                                                         arrayV[i][j]=t;
993                                                 }
994                                                 t=ExtraMath.hypot(f,g);
995                                                 cs=f/t;
996                                                 sn=g/t;
997                                                 arrayS[j]=t;
998                                                 f=cs*e[j]+sn*arrayS[j+1];
999                                                 arrayS[j+1]=-sn*e[j]+cs*arrayS[j+1];
1000                                                g=sn*e[j+1];
1001                                                e[j+1]*=cs;
1002                                                if(j<Nm1) {
1003                                                        for(int i=0;i<N;i++) {
1004                                                                t=cs*arrayU[i][j]+sn*arrayU[i][j+1];
1005                                                                arrayU[i][j+1]=-sn*arrayU[i][j]+cs*arrayU[i][j+1];
1006                                                                arrayU[i][j]=t;
1007                                                        }
1008                                                }
1009                                        }
1010                                        e[p-2]=f;
1011                                        iter++;
1012                                        } break;
1013                                // convergence
1014 case 4: {
1015                                        // make the singular values positive
1016 if(arrayS[k]<=0.0) {
1017                                                arrayS[k]=-arrayS[k];
1018                                                for(int i=0;i<p;i++)
1019                                                        arrayV[i][k]=-arrayV[i][k];
1020                                        }
1021                                        // order the singular values
1022 while(k<p-1) {
1023                                                if(arrayS[k]>=arrayS[k+1])
1024                                                        break;
1025                                                double tmp=arrayS[k];
1026                                                arrayS[k]=arrayS[k+1];
1027                                                arrayS[k+1]=tmp;
1028                                                if(k<Nm1) {
1029                                                        for(int i=0;i<N;i++) {
1030                                                                tmp=arrayU[i][k+1];
1031                                                                arrayU[i][k+1]=arrayU[i][k];
1032                                                                arrayU[i][k]=tmp;
1033                                                                tmp=arrayV[i][k+1];
1034                                                                arrayV[i][k+1]=arrayV[i][k];
1035                                                                arrayV[i][k]=tmp;
1036                                                        }
1037                                                }
1038                                                k++;
1039                                        }
1040                                        iter=0;
1041                                        p--;
1042                                        } break;
1043                        }
1044                }
1045                final AbstractDoubleSquareMatrix svd[]=new AbstractDoubleSquareMatrix[3];
1046                svd[0]=new DoubleSquareMatrix(arrayU);
1047                svd[1]=new DoubleDiagonalMatrix(arrayS);
1048                svd[2]=new DoubleSquareMatrix(arrayV);
1049                return svd;
1050        }
1051
1052}
1053
Popular Tags