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.219 * @author Mark Hale20 */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/columns31         */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 value42         * @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=0 && i=0 && j=0 && i=0 && j-norm.200         * @author Taber Smith201         */202         public double infNorm() {203                 double result=Math.abs(diag[0])+Math.abs(udiag[0]);204                 double tmpResult;205                 int i=1;206                 for(;iresult)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 (l2) norm.218         * @author Taber Smith219         */220         public double frobeniusNorm() {221                 double result=diag[0]*diag[0]+udiag[0]*udiag[0];222                 int i=1;223                 for(;i3) {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;i3) {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                 } else603                         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;i3) {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 // TRANSPOSE648 649         /**650         * Returns the transpose of this matrix.651         * @return a double matrix652         */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 DECOMPOSITION662 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;jJAMA (public domain).704         * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.705         * @jsci.planetmath QRDecomposition706         */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 matrix713 array[0][0]=diag[0];714                 array[0][1]=udiag[0];715                 for(int i=1;i=0; k--) {747                         arrayQ[k][k] = 1.0;748                         for(int j=k; jJAMA (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 matrix786 array[0][0]=diag[0];787                 array[0][1]=udiag[0];788                 for(int i=1;i=0;k--) {862                         if(arrayS[k]!=0.0) {863                                 for(int j=k+1;j=0;k--) {883                         if(k0) {900                         int k, action;901                         // action = 1 if arrayS[p] and e[k-1] are negligible and 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=arrayS[k+1])1043                                                        break;1044                                                double tmp=arrayS[k];1045                                                arrayS[k]=arrayS[k+1];1046                                                arrayS[k+1]=tmp;1047                                                if(k