package JSci.maths;

/**
* The ComplexTridiagonalMatrix class provides an object for encapsulating tridiagonal matrices containing complex numbers.
* Uses compressed diagonal storage.
* @version 2.2
* @author Mark Hale
*/
public class ComplexTridiagonalMatrix extends ComplexSquareMatrix implements TridiagonalMatrix {
//
// example: element storage for a 3x3 tridiagonal matrix
//
// matrix[1][0] matrix[2][0]
// matrix[0][1] matrix[1][1] matrix[2][1]
//              matrix[0][2] matrix[1][2]
//
        /**
        * Storage format identifier.
        */
        protected final static int TRIDIAGONAL=3;
        /**
        * Tridiagonal data.
        */
        protected double ldiagRe[],ldiagIm[];
        protected double diagRe[],diagIm[];
        protected double udiagRe[],udiagIm[];
        /**
        * Constructs a matrix.
        */
        protected ComplexTridiagonalMatrix(final int size,final int storeID) {
                super(size,storeID);
        }
        /**
        * Constructs an empty matrix.
        * @param size the number of rows/columns
        */
        public ComplexTridiagonalMatrix(final int size) {
                this(size,TRIDIAGONAL);
                matrixRe=new double[3][size];
                matrixIm=new double[3][size];
                ldiagRe=matrixRe[0];
                ldiagIm=matrixIm[0];
                diagRe=matrixRe[1];
                diagIm=matrixIm[1];
                udiagRe=matrixRe[2];
                udiagIm=matrixIm[2];
        }
        /**
        * Constructs a matrix from an array.
        * Any non-tridiagonal elements in the array are ignored.
        * @param array an assigned value
        * @exception MatrixDimensionException If the array is not square.
        */
        public ComplexTridiagonalMatrix(final Complex array[][]) {
                this(array.length);
                if(array.length==array[0].length) {
                        diagRe[0]=array[0][0].real();
                        diagIm[0]=array[0][0].imag();
                        udiagRe[0]=array[0][1].real();
                        udiagIm[0]=array[0][1].imag();
                        int i=1;
                        for(;i<array.length-1;i++) {
                                ldiagRe[i]=array[i][i-1].real();
                                ldiagIm[i]=array[i][i-1].imag();
                                diagRe[i]=array[i][i].real();
                                diagIm[i]=array[i][i].imag();
                                udiagRe[i]=array[i][i+1].real();
                                udiagIm[i]=array[i][i+1].imag();
                        }
                        ldiagRe[i]=array[i][i-1].real();
                        ldiagIm[i]=array[i][i-1].imag();
                        diagRe[i]=array[i][i].real();
                        diagIm[i]=array[i][i].imag();
                } else {
                        matrixRe=null;
                        matrixIm=null;
                        throw new MatrixDimensionException("The array is not square.");
                }
        }
        /**
        * Compares two complex tridiagonal matrices for equality.
        * @param m a complex tridiagonal matrix
        */
        public boolean equals(Object m) {
                if((m instanceof ComplexTridiagonalMatrix) &&
                numRows==((ComplexTridiagonalMatrix)m).rows()) {
                        final ComplexTridiagonalMatrix ctm=(ComplexTridiagonalMatrix)m;
                        if(!ctm.getElement(0,0).equals(diagRe[0],diagIm[0]))
                                return false;
                        if(!ctm.getElement(0,1).equals(udiagRe[0],udiagIm[0]))
                                return false;
                        int i=1;
                        for(;i<numRows-1;i++) {
                                if(!ctm.getElement(i,i-1).equals(ldiagRe[i],ldiagIm[i]))
                                        return false;
                                if(!ctm.getElement(i,i).equals(diagRe[i],diagIm[i]))
                                        return false;
                                if(!ctm.getElement(i,i+1).equals(udiagRe[i],udiagIm[i]))
                                        return false;
                        }
                        if(!ctm.getElement(i,i-1).equals(ldiagRe[i],ldiagIm[i]))
                                return false;
                        if(!ctm.getElement(i,i).equals(diagRe[i],diagIm[i]))
                                return false;
                        return true;
                } else
                        return false;
        }
        /**
        * Returns a string representing this matrix.
        */
        public String toString() {
                final StringBuffer buf=new StringBuffer(5*rows()*columns());
                for(int j,i=0;i<rows();i++) {
                        for(j=0;j<columns();j++) {
                                buf.append(getElement(i,j).toString());
                                buf.append(' ');
                        }
                        buf.append('\n');
                }
                return buf.toString();
        }
        /**
        * Returns the real part of this complex matrix.
        * @return a double tridiagonal matrix
        */
        public DoubleMatrix real() {
                final DoubleTridiagonalMatrix m=new DoubleTridiagonalMatrix(numRows);
                m.matrix[1][0]=diagRe[0];
                m.matrix[2][0]=udiagRe[0];
                int i=1;
                for(;i<numRows-1;i++) {
                        m.matrix[0][i]=ldiagRe[i];
                        m.matrix[1][i]=diagRe[i];
                        m.matrix[2][i]=udiagRe[i];
                }
                m.matrix[0][i]=ldiagRe[i];
                m.matrix[1][i]=diagRe[i];
                return m;
        }
        /**
        * Returns the imaginary part of this complex matrix.
        * @return a double tridiagonal matrix
        */
        public DoubleMatrix imag() {
                final DoubleTridiagonalMatrix m=new DoubleTridiagonalMatrix(numRows);
                m.matrix[1][0]=diagIm[0];
                m.matrix[2][0]=udiagIm[0];
                int i=1;
                for(;i<numRows-1;i++) {
                        m.matrix[0][i]=ldiagIm[i];
                        m.matrix[1][i]=diagIm[i];
                        m.matrix[2][i]=udiagIm[i];
                }
                m.matrix[0][i]=ldiagIm[i];
                m.matrix[1][i]=diagIm[i];
                return m;
        }
        /**
        * Returns an element of the matrix.
        * @param i row index of the element
        * @param j column index of the element
        * @exception MatrixDimensionException If attempting to access an invalid element.
        */
        public Complex getElement(final int i, final int j) {
                if(i>=0 && i<numRows && j>=0 && j<numCols) {
                        if(j>=i-1 && j<=i+1)
                                return new Complex(matrixRe[j-i+1][i],matrixIm[j-i+1][i]);
                        else
                                return Complex.ZERO;
                } else
                        throw new MatrixDimensionException(getInvalidElementMsg(i,j));
        }
        /**
        * Sets the value of an element of the matrix.
        * Should only be used to initialise this matrix.
        * @param i row index of the element
        * @param j column index of the element
        * @param z a complex number
        * @exception MatrixDimensionException If attempting to access an invalid element.
        */
        public void setElement(final int i, final int j, final Complex z) {
                if(i>=0 && i<numRows && j>=0 && j<numCols) {
                        if(j>=i-1 && j<=i+1) {
                                matrixRe[j-i+1][i]=z.real();
                                matrixIm[j-i+1][i]=z.imag();
                        } else
                                throw new MatrixDimensionException(getInvalidElementMsg(i,j));
                } else
                        throw new MatrixDimensionException(getInvalidElementMsg(i,j));
        }
        /**
        * Sets the value of an element of the matrix.
        * Should only be used to initialise this matrix.
        * @param i row index of the element
        * @param j column index of the element
        * @param x the real part of a complex number
        * @param y the imaginary part of a complex number
        * @exception MatrixDimensionException If attempting to access an invalid element.
        */
        public void setElement(final int i, final int j, final double x, final double y) {
                if(i>=0 && i<numRows && j>=0 && j<numCols) {
                        if(j>=i-1 && j<=i+1) {
                                matrixRe[j-i+1][i]=x;
                                matrixIm[j-i+1][i]=y;
                        } else
                                throw new MatrixDimensionException(getInvalidElementMsg(i,j));
                } else
                        throw new MatrixDimensionException(getInvalidElementMsg(i,j));
        }
        /**
        * Returns the trace.
        */
        public Complex trace() {
                double trRe=diagRe[0];
                double trIm=diagIm[0];
                for(int i=1;i<numRows;i++) {
                        trRe+=diagRe[i];
                        trIm+=diagIm[i];
                }
                return new Complex(trRe,trIm);
        }
        /**
        * Returns the l<sup><img border=0 alt="infinity" src="doc-files/infinity.gif"></sup>-norm.
        * @author Taber Smith
        */
        public double infNorm() {
                double result=#modA("diag" "[0]")+#modA("udiag" "[0]");
                double tmpResult;
                int i=1;
                for(;i<numRows-1;i++) {
                        tmpResult=#modA("ldiag" "[i]")+#modA("diag" "[i]")+#modA("udiag" "[i]");
                        if(tmpResult>result)
                                result=tmpResult;
                }
                tmpResult=#modA("ldiag" "[i]")+#modA("diag" "[i]");
                if(tmpResult>result)
                        result=tmpResult;
                return result;
        }
        /**
        * Returns the Frobenius (l<sup>2</sup>) norm.
        * @author Taber Smith
        */
        public double frobeniusNorm() {
                double result=diagRe[0]*diagRe[0]+diagIm[0]*diagIm[0]+
                        udiagRe[0]*udiagRe[0]+udiagIm[0]*udiagIm[0];
                int i=1;
                for(;i<numRows-1;i++) {
                        result+=ldiagRe[i]*ldiagRe[i]+ldiagIm[i]*ldiagIm[i]+
                                diagRe[i]*diagRe[i]+diagIm[i]*diagIm[i]+
                                udiagRe[i]*udiagRe[i]+udiagIm[i]*udiagIm[i];
                }
                result+=ldiagRe[i]*ldiagRe[i]+ldiagIm[i]*ldiagIm[i]+
                        diagRe[i]*diagRe[i]+diagIm[i]*diagIm[i];
                return Math.sqrt(result);
        }
        /**
        * Returns the operator norm.
        * @exception MaximumIterationsExceededException If it takes more than 50 iterations to determine an eigenvalue.
        */
        public double operatorNorm() throws MaximumIterationsExceededException {
                return Math.sqrt(ArrayMath.max(LinearMath.eigenvalueSolveHermitian((ComplexTridiagonalMatrix)(this.hermitianAdjoint().multiply(this)))));
        }

//============
// OPERATIONS
//============

// ADDITION

        /**
        * Returns the addition of this matrix and another.
        * @param m a complex matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexMatrix add(final ComplexMatrix m) {
                switch(m.storageFormat) {
                        case ARRAY_2D: return rawAdd(m);
                        case TRIDIAGONAL: return rawAddTridiagonal(m);
                        default: 
                                if(numRows==m.rows() && numCols==m.columns()) {
                                        final double arrayRe[][]=new double[numRows][numCols];
                                        final double arrayIm[][]=new double[numRows][numCols];
                                        Complex elem;
                                        for(int j,i=0;i<numRows;i++) {
                                                elem=getElement(i,0).add(m.getElement(i,0));
                                                arrayRe[i][0]=elem.real();
                                                arrayIm[i][0]=elem.imag();
                                                for(j=1;j<numCols;j++) {
                                                        elem=getElement(i,j).add(m.getElement(i,j));
                                                        arrayRe[i][j]=elem.real();
                                                        arrayIm[i][j]=elem.imag();
                                                }
                                        }
                                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                                } else
                                        throw new MatrixDimensionException("Matrices are different sizes.");
                }
        }
        private ComplexSquareMatrix rawAdd(final ComplexMatrix m) {
                if(numRows==m.numRows && numCols==m.numCols) {
                        final double arrayRe[][]=new double[numRows][numRows];
                        final double arrayIm[][]=new double[numRows][numRows];
                        for(int i=0;i<numRows;i++) {
                                System.arraycopy(m.matrixRe[i],0,arrayRe[i],0,numRows);
                                System.arraycopy(m.matrixIm[i],0,arrayIm[i],0,numRows);
                        }
                        arrayRe[0][0]+=diagRe[0];
                        arrayIm[0][0]+=diagIm[0];
                        arrayRe[0][1]+=udiagRe[0];
                        arrayIm[0][1]+=udiagIm[0];
                        int n=numCols-1;
                        for(int i=1;i<n;i++) {
                                arrayRe[i][i-1]+=ldiagRe[i];
                                arrayIm[i][i-1]+=ldiagIm[i];
                                arrayRe[i][i]+=diagRe[i];
                                arrayIm[i][i]+=diagIm[i];
                                arrayRe[i][i+1]+=udiagRe[i];
                                arrayIm[i][i+1]+=udiagIm[i];
                        }
                        arrayRe[n][n-1]+=ldiagRe[n];
                        arrayIm[n][n-1]+=ldiagIm[n];
                        arrayRe[n][n]+=diagRe[n];
                        arrayIm[n][n]+=diagIm[n];
                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                } else
                        throw new MatrixDimensionException("Matrices are different sizes.");
        }
        /**
        * Returns the addition of this matrix and another.
        * @param m a complex square matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexSquareMatrix add(final ComplexSquareMatrix m) {
                switch(m.storageFormat) {
                        case ARRAY_2D: return rawAdd(m);
                        case TRIDIAGONAL: return rawAddTridiagonal(m);
                        default: 
                                if(numRows==m.rows()) {
                                        final double arrayRe[][]=new double[numRows][numCols];
                                        final double arrayIm[][]=new double[numRows][numCols];
                                        Complex elem;
                                        for(int j,i=0;i<numRows;i++) {
                                                elem=getElement(i,0).add(m.getElement(i,0));
                                                arrayRe[i][0]=elem.real();
                                                arrayIm[i][0]=elem.imag();
                                                for(j=1;j<numCols;j++) {
                                                        elem=getElement(i,j).add(m.getElement(i,j));
                                                        arrayRe[i][j]=elem.real();
                                                        arrayIm[i][j]=elem.imag();
                                                }
                                        }
                                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                                } else
                                        throw new MatrixDimensionException("Matrices are different sizes.");
                }
        }
        /**
        * Returns the addition of this matrix and another.
        * @param m a complex tridiagonal matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexTridiagonalMatrix add(final ComplexTridiagonalMatrix m) {
                switch(m.storageFormat) {
                        case TRIDIAGONAL: return rawAddTridiagonal(m);
                        default: 
                                int mRow=numRows;
                                if(mRow==m.rows()) {
                                        final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                                        Complex elem=m.getElement(0,0);
                                        ans.diagRe[0]=diagRe[0]+elem.real();
                                        ans.diagIm[0]=diagIm[0]+elem.imag();
                                        elem=m.getElement(0,1);
                                        ans.udiagRe[0]=udiagRe[0]+elem.real();
                                        ans.udiagIm[0]=udiagIm[0]+elem.imag();
                                        mRow--;
                                        for(int i=1;i<mRow;i++) {
                                                elem=m.getElement(i,i-1);
                                                ans.ldiagRe[i]=ldiagRe[i]+elem.real();
                                                ans.ldiagIm[i]=ldiagIm[i]+elem.imag();
                                                elem=m.getElement(i,i);
                                                ans.diagRe[i]=diagRe[i]+elem.real();
                                                ans.diagIm[i]=diagIm[i]+elem.imag();
                                                elem=m.getElement(i,i+1);
                                                ans.udiagRe[i]=udiagRe[i]+elem.real();
                                                ans.udiagIm[i]=udiagIm[i]+elem.imag();
                                        }
                                        elem=m.getElement(mRow,mRow-1);
                                        ans.ldiagRe[mRow]=ldiagRe[mRow]+elem.real();
                                        ans.ldiagIm[mRow]=ldiagIm[mRow]+elem.imag();
                                        elem=m.getElement(mRow,mRow);
                                        ans.diagRe[mRow]=diagRe[mRow]+elem.real();
                                        ans.diagIm[mRow]=diagIm[mRow]+elem.imag();
                                        return ans;
                                } else
                                        throw new MatrixDimensionException("Matrices are different sizes.");
                }
        }
        private ComplexTridiagonalMatrix rawAddTridiagonal(final ComplexMatrix m) {
                int mRow=numRows;
                if(mRow==m.numRows) {
                        final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                        ans.diagRe[0]=diagRe[0]+m.matrixRe[1][0];
                        ans.diagIm[0]=diagIm[0]+m.matrixIm[1][0];
                        ans.udiagRe[0]=udiagRe[0]+m.matrixRe[2][0];
                        ans.udiagIm[0]=udiagIm[0]+m.matrixIm[2][0];
                        mRow--;
                        for(int i=1;i<mRow;i++) {
                                ans.ldiagRe[i]=ldiagRe[i]+m.matrixRe[0][i];
                                ans.ldiagIm[i]=ldiagIm[i]+m.matrixIm[0][i];
                                ans.diagRe[i]=diagRe[i]+m.matrixRe[1][i];
                                ans.diagIm[i]=diagIm[i]+m.matrixIm[1][i];
                                ans.udiagRe[i]=udiagRe[i]+m.matrixRe[2][i];
                                ans.udiagIm[i]=udiagIm[i]+m.matrixIm[2][i];
                        }
                        ans.ldiagRe[mRow]=ldiagRe[mRow]+m.matrixRe[0][mRow];
                        ans.ldiagIm[mRow]=ldiagIm[mRow]+m.matrixIm[0][mRow];
                        ans.diagRe[mRow]=diagRe[mRow]+m.matrixRe[1][mRow];
                        ans.diagIm[mRow]=diagIm[mRow]+m.matrixIm[1][mRow];
                        return ans;
                } else
                        throw new MatrixDimensionException("Matrices are different sizes.");
        }

// SUBTRACTION

        /**
        * Returns the subtraction of this matrix and another.
        * @param m a complex matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexMatrix subtract(final ComplexMatrix m) {
                switch(m.storageFormat) {
                        case ARRAY_2D: return rawSubtract(m);
                        case TRIDIAGONAL: return rawSubtractTridiagonal(m);
                        default: 
                                if(numRows==m.rows() && numRows==m.columns()) {
                                        final double arrayRe[][]=new double[numRows][numCols];
                                        final double arrayIm[][]=new double[numRows][numCols];
                                        Complex elem;
                                        for(int j,i=0;i<numRows;i++) {
                                                elem=getElement(i,0).subtract(m.getElement(i,0));
                                                arrayRe[i][0]=elem.real();
                                                arrayIm[i][0]=elem.imag();
                                                for(j=1;j<numCols;j++) {
                                                        elem=getElement(i,j).subtract(m.getElement(i,j));
                                                        arrayRe[i][j]=elem.real();
                                                        arrayIm[i][j]=elem.imag();
                                                }
                                        }
                                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                                } else
                                        throw new MatrixDimensionException("Matrices are different sizes.");
                }
        }
        private ComplexSquareMatrix rawSubtract(final ComplexMatrix m) {
                if(numRows==m.numRows && numCols==m.numCols) {
                        final double arrayRe[][]=new double[numRows][numCols];
                        final double arrayIm[][]=new double[numRows][numCols];
                        Complex elem;
                        for(int j,i=0;i<numRows;i++) {
                                elem=getElement(i,0);
                                arrayRe[i][0]=elem.real()-m.matrixRe[i][0];
                                arrayIm[i][0]=elem.imag()-m.matrixIm[i][0];
                                for(j=1;j<numCols;j++) {
                                        elem=getElement(i,j);
                                        arrayRe[i][j]=elem.real()-m.matrixRe[i][j];
                                        arrayIm[i][j]=elem.imag()-m.matrixIm[i][j];
                                }
                        }
                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                } else
                        throw new MatrixDimensionException("Matrices are different sizes.");
        }
        /**
        * Returns the subtraction of this matrix by another.
        * @param m a complex square matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexSquareMatrix subtract(final ComplexSquareMatrix m) {
                switch(m.storageFormat) {
                        case ARRAY_2D: return rawSubtract(m);
                        case TRIDIAGONAL: return rawSubtractTridiagonal(m);
                        default: 
                                if(numRows==m.rows()) {
                                        final double arrayRe[][]=new double[numRows][numCols];
                                        final double arrayIm[][]=new double[numRows][numCols];
                                        Complex elem;
                                        for(int j,i=0;i<numRows;i++) {
                                                elem=getElement(i,0).subtract(m.getElement(i,0));
                                                arrayRe[i][0]=elem.real();
                                                arrayIm[i][0]=elem.imag();
                                                for(j=1;j<numCols;j++) {
                                                        elem=getElement(i,j).subtract(m.getElement(i,j));
                                                        arrayRe[i][j]=elem.real();
                                                        arrayIm[i][j]=elem.imag();
                                                }
                                        }
                                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                                } else
                                        throw new MatrixDimensionException("Matrices are different sizes.");
                }
        }
        /**
        * Returns the subtraction of this matrix and another.
        * @param m a complex tridiagonal matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexTridiagonalMatrix subtract(final ComplexTridiagonalMatrix m) {
                switch(m.storageFormat) {
                        case TRIDIAGONAL: return rawSubtractTridiagonal(m);
                        default: 
                                int mRow=numRows;
                                if(mRow==m.rows()) {
                                        final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                                        Complex elem=m.getElement(0,0);
                                        ans.diagRe[0]=diagRe[0]-elem.real();
                                        ans.diagIm[0]=diagIm[0]-elem.imag();
                                        elem=m.getElement(0,1);
                                        ans.udiagRe[0]=udiagRe[0]-elem.real();
                                        ans.udiagIm[0]=udiagIm[0]-elem.imag();
                                        mRow--;
                                        for(int i=1;i<mRow;i++) {
                                                elem=m.getElement(i,i-1);
                                                ans.ldiagRe[i]=ldiagRe[i]-elem.real();
                                                ans.ldiagIm[i]=ldiagIm[i]-elem.imag();
                                                elem=m.getElement(i,i);
                                                ans.diagRe[i]=diagRe[i]-elem.real();
                                                ans.diagIm[i]=diagIm[i]-elem.imag();
                                                elem=m.getElement(i,i+1);
                                                ans.udiagRe[i]=udiagRe[i]-elem.real();
                                                ans.udiagIm[i]=udiagIm[i]-elem.imag();
                                        }
                                        elem=m.getElement(mRow,mRow-1);
                                        ans.ldiagRe[mRow]=ldiagRe[mRow]-elem.real();
                                        ans.ldiagIm[mRow]=ldiagIm[mRow]-elem.imag();
                                        elem=m.getElement(mRow,mRow);
                                        ans.diagRe[mRow]=diagRe[mRow]-elem.real();
                                        ans.diagIm[mRow]=diagIm[mRow]-elem.imag();
                                        return ans;
                                } else
                                        throw new MatrixDimensionException("Matrices are different sizes.");
                }
        }
        private ComplexTridiagonalMatrix rawSubtractTridiagonal(final ComplexMatrix m) {
                int mRow=numRows;
                if(mRow==m.numRows) {
                        final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                        ans.diagRe[0]=diagRe[0]-m.matrixRe[1][0];
                        ans.diagIm[0]=diagIm[0]-m.matrixIm[1][0];
                        ans.udiagRe[0]=udiagRe[0]-m.matrixRe[2][0];
                        ans.udiagIm[0]=udiagIm[0]-m.matrixIm[2][0];
                        mRow--;
                        for(int i=1;i<mRow;i++) {
                                ans.ldiagRe[i]=ldiagRe[i]-m.matrixRe[0][i];
                                ans.ldiagIm[i]=ldiagIm[i]-m.matrixIm[0][i];
                                ans.diagRe[i]=diagRe[i]-m.matrixRe[1][i];
                                ans.diagIm[i]=diagIm[i]-m.matrixIm[1][i];
                                ans.udiagRe[i]=udiagRe[i]-m.matrixRe[2][i];
                                ans.udiagIm[i]=udiagIm[i]-m.matrixIm[2][i];
                        }
                        ans.ldiagRe[mRow]=ldiagRe[mRow]-m.matrixRe[0][mRow];
                        ans.ldiagIm[mRow]=ldiagIm[mRow]-m.matrixIm[0][mRow];
                        ans.diagRe[mRow]=diagRe[mRow]-m.matrixRe[1][mRow];
                        ans.diagIm[mRow]=diagIm[mRow]-m.matrixIm[1][mRow];
                        return ans;
                } else
                        throw new MatrixDimensionException("Matrices are different sizes.");
        }

// SCALAR MULTIPLICATION

        /**
        * Returns the multiplication of this matrix by a scalar.
        * @param z a complex number
        * @return a complex tridiagonal matrix
        */
        public ComplexMatrix scalarMultiply(final Complex z) {
                final double real=z.real();
                final double imag=z.imag();
                int mRow=numRows;
                final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                ans.diagRe[0]=real*diagRe[0]-imag*diagIm[0];
                ans.diagIm[0]=imag*diagRe[0]+real*diagIm[0];
                ans.udiagRe[0]=real*udiagRe[0]-imag*udiagIm[0];
                ans.udiagIm[0]=imag*udiagRe[0]+real*udiagIm[0];
                mRow--;
                for(int i=1;i<mRow;i++) {
                        ans.ldiagRe[i]=real*ldiagRe[i]-imag*ldiagIm[i];
                        ans.ldiagIm[i]=imag*ldiagRe[i]+real*ldiagIm[i];
                        ans.diagRe[i]=real*diagRe[i]-imag*diagIm[i];
                        ans.diagIm[i]=imag*diagRe[i]+real*diagIm[i];
                        ans.udiagRe[i]=real*udiagRe[i]-imag*udiagIm[i];
                        ans.udiagIm[i]=imag*udiagRe[i]+real*udiagIm[i];
                }
                ans.ldiagRe[mRow]=real*ldiagRe[mRow]-imag*ldiagIm[mRow];
                ans.ldiagIm[mRow]=imag*ldiagRe[mRow]+real*ldiagIm[mRow];
                ans.diagRe[mRow]=real*diagRe[mRow]-imag*diagIm[mRow];
                ans.diagIm[mRow]=imag*diagRe[mRow]+real*diagIm[mRow];
                return ans;
        }
        /**
        * Returns the multiplication of this matrix by a scalar.
        * @param x a double
        * @return a complex tridiagonal matrix
        */
        public ComplexMatrix scalarMultiply(final double x) {
                int mRow=numRows;
                final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                ans.diagRe[0]=x*diagRe[0];
                ans.diagIm[0]=x*diagIm[0];
                ans.udiagRe[0]=x*udiagRe[0];
                ans.udiagIm[0]=x*udiagIm[0];
                mRow--;
                for(int i=1;i<mRow;i++) {
                        ans.ldiagRe[i]=x*ldiagRe[i];
                        ans.ldiagIm[i]=x*ldiagIm[i];
                        ans.diagRe[i]=x*diagRe[i];
                        ans.diagIm[i]=x*diagIm[i];
                        ans.udiagRe[i]=x*udiagRe[i];
                        ans.udiagIm[i]=x*udiagIm[i];
                }
                ans.ldiagRe[mRow]=x*ldiagRe[mRow];
                ans.ldiagIm[mRow]=x*ldiagIm[mRow];
                ans.diagRe[mRow]=x*diagRe[mRow];
                ans.diagIm[mRow]=x*diagIm[mRow];
                return ans;
        }

// MATRIX MULTIPLICATION

        /**
        * Returns the multiplication of a vector by this matrix.
        * @param v a complex vector
        * @exception DimensionException If the matrix and vector are incompatible.
        */
        public ComplexVector multiply(final ComplexVector v) {
                int mRow=numRows;
                if(mRow==v.dimension()) {
                        final double arrayRe[]=new double[mRow];
                        final double arrayIm[]=new double[mRow];
                        Complex comp;
                        comp=v.getComponent(0);
                        arrayRe[0]=#multAZRe("diag" "[0]" "comp");
                        arrayIm[0]=#multAZIm("diag" "[0]" "comp");
                        comp=v.getComponent(1);
                        arrayRe[0]+=#multAZRe("udiag" "[0]" "comp");
                        arrayIm[0]+=#multAZIm("udiag" "[0]" "comp");
                        mRow--;
                        for(int i=1;i<mRow;i++) {
                                comp=v.getComponent(i-1);
                                arrayRe[i]=#multAZRe("ldiag" "[i]" "comp");
                                arrayIm[i]=#multAZIm("ldiag" "[i]" "comp");
                                comp=v.getComponent(i);
                                arrayRe[i]+=#multAZRe("diag" "[i]" "comp");
                                arrayIm[i]+=#multAZIm("diag" "[i]" "comp");
                                comp=v.getComponent(i+1);
                                arrayRe[i]+=#multAZRe("udiag" "[i]" "comp");
                                arrayIm[i]+=#multAZIm("udiag" "[i]" "comp");
                        }
                        comp=v.getComponent(mRow-1);
                        arrayRe[mRow]=#multAZRe("ldiag" "[mRow]" "comp");
                        arrayIm[mRow]=#multAZIm("ldiag" "[mRow]" "comp");
                        comp=v.getComponent(mRow);
                        arrayRe[mRow]+=#multAZRe("diag" "[mRow]" "comp");
                        arrayIm[mRow]+=#multAZIm("diag" "[mRow]" "comp");
                        return new ComplexVector(arrayRe,arrayIm);
                } else
                        throw new DimensionException("Matrix and vector are incompatible.");
        }
        /**
        * Returns the multiplication of this matrix and another.
        * @param m a complex matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexMatrix multiply(final ComplexMatrix m) {
                switch(m.storageFormat) {
                        case ARRAY_2D: return rawMultiply(m);
                        case TRIDIAGONAL: return rawMultiplyTridiagonal(m);
                        default: 
                                if(numCols==m.rows()) {
                                        int n,k;
                                        final double arrayRe[][]=new double[numRows][m.columns()];
                                        final double arrayIm[][]=new double[numRows][m.columns()];
                                        Complex elem;
                                        for(int j=0;j<numRows;j++) {
                                                for(k=0;k<m.columns();k++) {
                                                        elem=getElement(j,0).multiply(m.getElement(0,k));
                                                        arrayRe[j][k]=elem.real();
                                                        arrayIm[j][k]=elem.imag();
                                                        for(n=1;n<numCols;n++) {
                                                                elem=getElement(j,n).multiply(m.getElement(n,k));
                                                                arrayRe[j][k]+=elem.real();
                                                                arrayIm[j][k]+=elem.imag();
                                                        }
                                                }
                                        }
                                        return new ComplexMatrix(arrayRe,arrayIm);
                                } else
                                        throw new MatrixDimensionException("Incompatible matrices.");
                }
        }
        private ComplexMatrix rawMultiply(final ComplexMatrix m) {
                if(numCols==m.numRows) {
                        int n,k;
                        final double arrayRe[][]=new double[numRows][m.numCols];
                        final double arrayIm[][]=new double[numRows][m.numCols];
                        double real,imag;
                        Complex elem;
                        for(int j=0;j<numRows;j++) {
                                for(k=0;k<m.numCols;k++) {
                                        elem=getElement(j,0);
                                        real=#multAZRe("m.matrix" "[0][k]" "elem");
                                        imag=#multAZIm("m.matrix" "[0][k]" "elem");
                                        for(n=1;n<numCols;n++) {
                                                elem=getElement(j,n);
                                                real+=#multAZRe("m.matrix" "[n][k]" "elem");
                                                imag+=#multAZIm("m.matrix" "[n][k]" "elem");
                                        }
                                        arrayRe[j][k]=real;
                                        arrayIm[j][k]=imag;
                                }
                        }
                        return new ComplexMatrix(arrayRe,arrayIm);
                } else
                        throw new MatrixDimensionException("Incompatible matrices.");
        }
        /**
        * Returns the multiplication of this matrix and another.
        * @param m a complex square matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexSquareMatrix multiply(final ComplexSquareMatrix m) {
                switch(m.storageFormat) {
                        case ARRAY_2D: return rawMultiply(m);
                        case TRIDIAGONAL: return rawMultiplyTridiagonal(m);
                        default: 
                                if(numCols==m.rows()) {
                                        int n,k;
                                        final double arrayRe[][]=new double[numRows][numCols];
                                        final double arrayIm[][]=new double[numRows][numCols];
                                        Complex elem;
                                        for(int j=0;j<numRows;j++) {
                                                for(k=0;k<numCols;k++) {
                                                        elem=getElement(j,0).multiply(m.getElement(0,k));
                                                        arrayRe[j][k]=elem.real();
                                                        arrayIm[j][k]=elem.imag();
                                                        for(n=1;n<numCols;n++) {
                                                                elem=getElement(j,n).multiply(m.getElement(n,k));
                                                                arrayRe[j][k]+=elem.real();
                                                                arrayIm[j][k]+=elem.imag();
                                                        }
                                                }
                                        }
                                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                                } else
                                        throw new MatrixDimensionException("Incompatible matrices.");
                }
        }
        private ComplexSquareMatrix rawMultiply(final ComplexSquareMatrix m) {
                if(numCols==m.numRows) {
                        int n,k;
                        final double arrayRe[][]=new double[numRows][numCols];
                        final double arrayIm[][]=new double[numRows][numCols];
                        double real,imag;
                        Complex elem;
                        for(int j=0;j<numRows;j++) {
                                for(k=0;k<numCols;k++) {
                                        elem=getElement(j,0);
                                        real=#multAZRe("m.matrix" "[0][k]" "elem");
                                        imag=#multAZIm("m.matrix" "[0][k]" "elem");
                                        for(n=1;n<numCols;n++) {
                                                elem=getElement(j,n);
                                                real+=#multAZRe("m.matrix" "[n][k]" "elem");
                                                imag+=#multAZIm("m.matrix" "[n][k]" "elem");
                                        }
                                        arrayRe[j][k]=real;
                                        arrayIm[j][k]=imag;
                                }
                        }
                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                } else
                        throw new MatrixDimensionException("Incompatible matrices.");
        }
        /**
        * Returns the multiplication of this matrix and another.
        * @param m a complex tridiagonal matrix
        * @exception MatrixDimensionException If the matrices are different sizes.
        */
        public ComplexSquareMatrix multiply(final ComplexTridiagonalMatrix m) {
                switch(m.storageFormat) {
                        case TRIDIAGONAL: return rawMultiplyTridiagonal(m);
                        default: 
                                int mRow=numRows;
                                if(numCols==m.rows()) {
                                        final double arrayRe[][]=new double[mRow][mRow];
                                        final double arrayIm[][]=new double[mRow][mRow];
                                        Complex elem1,elem2,elem3;
                                        elem1=m.getElement(0,0);elem2=m.getElement(1,0);
                                        arrayRe[0][0]=#multAZRe("diag" "[0]" "elem1")+#multAZRe("udiag" "[0]" "elem2");
                                        arrayIm[0][0]=#multAZIm("diag" "[0]" "elem1")+#multAZIm("udiag" "[0]" "elem2");
                                        elem1=m.getElement(0,1);elem2=m.getElement(1,1);
                                        arrayRe[0][1]=#multAZRe("diag" "[0]" "elem1")+#multAZRe("udiag" "[0]" "elem2");
                                        arrayIm[0][1]=#multAZIm("diag" "[0]" "elem1")+#multAZIm("udiag" "[0]" "elem2");
                                        elem1=m.getElement(1,2);
                                        arrayRe[0][2]=#multAZRe("udiag" "[0]" "elem1");
                                        arrayIm[0][2]=#multAZIm("udiag" "[0]" "elem1");
                                        if(mRow>3) {
                                                elem1=m.getElement(0,0);elem2=m.getElement(1,0);
                                                arrayRe[1][0]=#multAZRe("ldiag" "[1]" "elem1")+#multAZRe("diag" "[1]" "elem2");
                                                arrayIm[1][0]=#multAZIm("ldiag" "[1]" "elem1")+#multAZIm("diag" "[1]" "elem2");
                                                elem1=m.getElement(0,1);elem2=m.getElement(1,1);elem3=m.getElement(2,1);
                                                arrayRe[1][1]=#multAZRe("ldiag" "[1]" "elem1")+#multAZRe("diag" "[1]" "elem2")+#multAZRe("udiag" "[1]" "elem3");
                                                arrayIm[1][1]=#multAZIm("ldiag" "[1]" "elem1")+#multAZIm("diag" "[1]" "elem2")+#multAZIm("udiag" "[1]" "elem3");
                                                elem1=m.getElement(1,2);elem2=m.getElement(2,2);
                                                arrayRe[1][2]=#multAZRe("diag" "[1]" "elem1")+#multAZRe("udiag" "[1]" "elem2");
                                                arrayIm[1][2]=#multAZIm("diag" "[1]" "elem1")+#multAZIm("udiag" "[1]" "elem2");
                                                elem1=m.getElement(2,3);
                                                arrayRe[1][3]=#multAZRe("udiag" "[1]" "elem1");
                                                arrayIm[1][3]=#multAZIm("udiag" "[1]" "elem1");
                                        }
                                        if(mRow==3) {
                                                elem1=m.getElement(0,0);elem2=m.getElement(1,0);
                                                arrayRe[1][0]=#multAZRe("ldiag" "[1]" "elem1")+#multAZRe("diag" "[1]" "elem2");
                                                arrayIm[1][0]=#multAZIm("ldiag" "[1]" "elem1")+#multAZIm("diag" "[1]" "elem2");
                                                elem1=m.getElement(0,1);elem2=m.getElement(1,1);elem3=m.getElement(2,1);
                                                arrayRe[1][1]=#multAZRe("ldiag" "[1]" "elem1")+#multAZRe("diag" "[1]" "elem2")+#multAZRe("udiag" "[1]" "elem3");
                                                arrayIm[1][1]=#multAZIm("ldiag" "[1]" "elem1")+#multAZIm("diag" "[1]" "elem2")+#multAZIm("udiag" "[1]" "elem3");
                                                elem1=m.getElement(1,2);elem2=m.getElement(2,2);
                                                arrayRe[1][2]=#multAZRe("diag" "[1]" "elem1")+#multAZRe("udiag" "[1]" "elem2");
                                                arrayIm[1][2]=#multAZIm("diag" "[1]" "elem1")+#multAZIm("udiag" "[1]" "elem2");
                                        } else if(mRow>4) {
                                                for(int i=2;i<mRow-2;i++) {
                                                        elem1=m.getElement(i-1,i-2);
                                                        arrayRe[i][i-2]=#multAZRe("ldiag" "[i]" "elem1");
                                                        arrayIm[i][i-2]=#multAZIm("ldiag" "[i]" "elem1");
                                                        elem1=m.getElement(i-1,i-1);elem2=m.getElement(i,i-1);
                                                        arrayRe[i][i-1]=#multAZRe("ldiag" "[i]" "elem1")+#multAZRe("diag" "[i]" "elem2");
                                                        arrayIm[i][i-1]=#multAZIm("ldiag" "[i]" "elem1")+#multAZIm("diag" "[i]" "elem2");
                                                        elem1=m.getElement(i-1,i);elem2=m.getElement(i,i);elem3=m.getElement(i+1,i);
                                                        arrayRe[i][i]=#multAZRe("ldiag" "[i]" "elem1")+#multAZRe("diag" "[i]" "elem2")+#multAZRe("udiag" "[i]" "elem3");
                                                        arrayIm[i][i]=#multAZIm("ldiag" "[i]" "elem1")+#multAZIm("diag" "[i]" "elem2")+#multAZIm("udiag" "[i]" "elem3");
                                                        elem1=m.getElement(i,i+1);elem2=m.getElement(i+1,i+1);
                                                        arrayRe[i][i+1]=#multAZRe("diag" "[i]" "elem1")+#multAZRe("udiag" "[i]" "elem2");
                                                        arrayIm[i][i+1]=#multAZIm("diag" "[i]" "elem1")+#multAZIm("udiag" "[i]" "elem2");
                                                        elem1=m.getElement(i+1,i+2);
                                                        arrayRe[i][i+2]=#multAZRe("udiag" "[i]" "elem1");
                                                        arrayIm[i][i+2]=#multAZIm("udiag" "[i]" "elem1");
                                                }
                                        }
                                        if(mRow>3) {
                                                elem1=m.getElement(mRow-3,mRow-4);
                                                arrayRe[mRow-2][mRow-4]=#multAZRe("ldiag" "[mRow-2]" "elem1");
                                                arrayIm[mRow-2][mRow-4]=#multAZIm("ldiag" "[mRow-2]" "elem1");
                                                elem1=m.getElement(mRow-3,mRow-3);elem2=m.getElement(mRow-2,mRow-3);
                                                arrayRe[mRow-2][mRow-3]=#multAZRe("ldiag" "[mRow-2]" "elem1")+#multAZRe("diag" "[mRow-2]" "elem2");
                                                arrayIm[mRow-2][mRow-3]=#multAZIm("ldiag" "[mRow-2]" "elem1")+#multAZIm("diag" "[mRow-2]" "elem2");
                                                elem1=m.getElement(mRow-3,mRow-2);elem2=m.getElement(mRow-2,mRow-2);elem3=m.getElement(mRow-1,mRow-2);
                                                arrayRe[mRow-2][mRow-2]=#multAZRe("ldiag" "[mRow-2]" "elem1")+#multAZRe("diag" "[mRow-2]" "elem2")+#multAZRe("udiag" "[mRow-2]" "elem3");
                                                arrayIm[mRow-2][mRow-2]=#multAZIm("ldiag" "[mRow-2]" "elem1")+#multAZIm("diag" "[mRow-2]" "elem2")+#multAZIm("udiag" "[mRow-2]" "elem3");
                                                elem1=m.getElement(mRow-2,mRow-1);elem2=m.getElement(mRow-1,mRow-1);
                                                arrayRe[mRow-2][mRow-1]=#multAZRe("diag" "[mRow-2]" "elem1")+#multAZRe("udiag" "[mRow-2]" "elem2");
                                                arrayIm[mRow-2][mRow-1]=#multAZIm("diag" "[mRow-2]" "elem1")+#multAZIm("udiag" "[mRow-2]" "elem2");
                                        }
                                        mRow--;
                                        elem1=m.getElement(mRow-1,mRow-2);
                                        arrayRe[mRow][mRow-2]=#multAZRe("ldiag" "[mRow]" "elem1");
                                        arrayIm[mRow][mRow-2]=#multAZIm("ldiag" "[mRow]" "elem1");
                                        elem1=m.getElement(mRow-1,mRow-1);elem2=m.getElement(mRow,mRow-1);
                                        arrayRe[mRow][mRow-1]=#multAZRe("ldiag" "[mRow]" "elem1")+#multAZRe("diag" "[mRow]" "elem2");
                                        arrayIm[mRow][mRow-1]=#multAZIm("ldiag" "[mRow]" "elem1")+#multAZIm("diag" "[mRow]" "elem2");
                                        elem1=m.getElement(mRow-1,mRow);elem2=m.getElement(mRow,mRow);
                                        arrayRe[mRow][mRow]=#multAZRe("ldiag" "[mRow]" "elem1")+#multAZRe("diag" "[mRow]" "elem2");
                                        arrayIm[mRow][mRow]=#multAZIm("ldiag" "[mRow]" "elem1")+#multAZIm("diag" "[mRow]" "elem2");
                                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                                } else
                                        throw new MatrixDimensionException("Incompatible matrices.");
                }
        }
        private ComplexSquareMatrix rawMultiplyTridiagonal(final ComplexMatrix m) {
                int mRow=numRows;
                if(numCols==m.numRows) {
                        final double arrayRe[][]=new double[mRow][mRow];
                        final double arrayIm[][]=new double[mRow][mRow];
                        arrayRe[0][0]=#multAARe("diag" "[0]" "m.matrix" "[1][0]")+#multAARe("udiag" "[0]" "m.matrix" "[0][1]");
                        arrayIm[0][0]=#multAAIm("diag" "[0]" "m.matrix" "[1][0]")+#multAAIm("udiag" "[0]" "m.matrix" "[0][1]");
                        arrayRe[0][1]=#multAARe("diag" "[0]" "m.matrix" "[2][0]")+#multAARe("udiag" "[0]" "m.matrix" "[1][1]");
                        arrayIm[0][1]=#multAAIm("diag" "[0]" "m.matrix" "[2][0]")+#multAAIm("udiag" "[0]" "m.matrix" "[1][1]");
                        arrayRe[0][2]=#multAARe("udiag" "[0]" "m.matrix" "[2][1]");
                        arrayIm[0][2]=#multAAIm("udiag" "[0]" "m.matrix" "[2][1]");
                        if(mRow>3) {
                                arrayRe[1][0]=#multAARe("ldiag" "[1]" "m.matrix" "[1][0]")+#multAARe("diag" "[1]" "m.matrix" "[0][1]");
                                arrayIm[1][0]=#multAAIm("ldiag" "[1]" "m.matrix" "[1][0]")+#multAAIm("diag" "[1]" "m.matrix" "[0][1]");
                                arrayRe[1][1]=#multAARe("ldiag" "[1]" "m.matrix" "[2][0]")+#multAARe("diag" "[1]" "m.matrix" "[1][1]")+#multAARe("udiag" "[1]" "m.matrix" "[0][2]");
                                arrayIm[1][1]=#multAAIm("ldiag" "[1]" "m.matrix" "[2][0]")+#multAAIm("diag" "[1]" "m.matrix" "[1][1]")+#multAAIm("udiag" "[1]" "m.matrix" "[0][2]");
                                arrayRe[1][2]=#multAARe("diag" "[1]" "m.matrix" "[2][1]")+#multAARe("udiag" "[1]" "m.matrix" "[1][2]");
                                arrayIm[1][2]=#multAAIm("diag" "[1]" "m.matrix" "[2][1]")+#multAAIm("udiag" "[1]" "m.matrix" "[1][2]");
                                arrayRe[1][3]=#multAARe("udiag" "[1]" "m.matrix" "[2][2]");
                                arrayIm[1][3]=#multAAIm("udiag" "[1]" "m.matrix" "[2][2]");
                        }
                        if(mRow==3) {
                                arrayRe[1][0]=#multAARe("ldiag" "[1]" "m.matrix" "[1][0]")+#multAARe("diag" "[1]" "m.matrix" "[0][1]");
                                arrayIm[1][0]=#multAAIm("ldiag" "[1]" "m.matrix" "[1][0]")+#multAAIm("diag" "[1]" "m.matrix" "[0][1]");
                                arrayRe[1][1]=#multAARe("ldiag" "[1]" "m.matrix" "[2][0]")+#multAARe("diag" "[1]" "m.matrix" "[1][1]")+#multAARe("udiag" "[1]" "m.matrix" "[0][2]");
                                arrayIm[1][1]=#multAAIm("ldiag" "[1]" "m.matrix" "[2][0]")+#multAAIm("diag" "[1]" "m.matrix" "[1][1]")+#multAAIm("udiag" "[1]" "m.matrix" "[0][2]");
                                arrayRe[1][2]=#multAARe("diag" "[1]" "m.matrix" "[2][1]")+#multAARe("udiag" "[1]" "m.matrix" "[1][2]");
                                arrayIm[1][2]=#multAAIm("diag" "[1]" "m.matrix" "[2][1]")+#multAAIm("udiag" "[1]" "m.matrix" "[1][2]");
                        } else if(mRow>4) {
                                for(int i=2;i<mRow-2;i++) {
                                        arrayRe[i][i-2]=#multAARe("ldiag" "[i]" "m.matrix" "[0][i-1]");
                                        arrayIm[i][i-2]=#multAAIm("ldiag" "[i]" "m.matrix" "[0][i-1]");
                                        arrayRe[i][i-1]=#multAARe("ldiag" "[i]" "m.matrix" "[1][i-1]")+#multAARe("diag" "[i]" "m.matrix" "[0][i]");
                                        arrayIm[i][i-1]=#multAAIm("ldiag" "[i]" "m.matrix" "[1][i-1]")+#multAAIm("diag" "[i]" "m.matrix" "[0][i]");
                                        arrayRe[i][i]=#multAARe("ldiag" "[i]" "m.matrix" "[2][i-1]")+#multAARe("diag" "[i]" "m.matrix" "[1][i]")+#multAARe("udiag" "[i]" "m.matrix" "[0][i+1]");
                                        arrayIm[i][i]=#multAAIm("ldiag" "[i]" "m.matrix" "[2][i-1]")+#multAAIm("diag" "[i]" "m.matrix" "[1][i]")+#multAAIm("udiag" "[i]" "m.matrix" "[0][i+1]");
                                        arrayRe[i][i+1]=#multAARe("diag" "[i]" "m.matrix" "[2][i]")+#multAARe("udiag" "[i]" "m.matrix" "[1][i+1]");
                                        arrayIm[i][i+1]=#multAAIm("diag" "[i]" "m.matrix" "[2][i]")+#multAAIm("udiag" "[i]" "m.matrix" "[1][i+1]");
                                        arrayRe[i][i+2]=#multAARe("udiag" "[i]" "m.matrix" "[2][i+1]");
                                        arrayIm[i][i+2]=#multAAIm("udiag" "[i]" "m.matrix" "[2][i+1]");
                                }
                        }
                        if(mRow>3) {
                                arrayRe[mRow-2][mRow-4]=#multAARe("ldiag" "[mRow-2]" "m.matrix" "[0][mRow-3]");
                                arrayIm[mRow-2][mRow-4]=#multAAIm("ldiag" "[mRow-2]" "m.matrix" "[0][mRow-3]");
                                arrayRe[mRow-2][mRow-3]=#multAARe("ldiag" "[mRow-2]" "m.matrix" "[1][mRow-3]")+#multAARe("diag" "[mRow-2]" "m.matrix" "[0][mRow-2]");
                                arrayIm[mRow-2][mRow-3]=#multAAIm("ldiag" "[mRow-2]" "m.matrix" "[1][mRow-3]")+#multAAIm("diag" "[mRow-2]" "m.matrix" "[0][mRow-2]");
                                arrayRe[mRow-2][mRow-2]=#multAARe("ldiag" "[mRow-2]" "m.matrix" "[2][mRow-3]")+#multAARe("diag" "[mRow-2]" "m.matrix" "[1][mRow-2]")+#multAARe("udiag" "[mRow-2]" "m.matrix" "[0][mRow-1]");
                                arrayIm[mRow-2][mRow-2]=#multAAIm("ldiag" "[mRow-2]" "m.matrix" "[2][mRow-3]")+#multAAIm("diag" "[mRow-2]" "m.matrix" "[1][mRow-2]")+#multAAIm("udiag" "[mRow-2]" "m.matrix" "[0][mRow-1]");
                                arrayRe[mRow-2][mRow-1]=#multAARe("diag" "[mRow-2]" "m.matrix" "[2][mRow-2]")+#multAARe("udiag" "[mRow-2]" "m.matrix" "[1][mRow-1]");
                                arrayIm[mRow-2][mRow-1]=#multAAIm("diag" "[mRow-2]" "m.matrix" "[2][mRow-2]")+#multAAIm("udiag" "[mRow-2]" "m.matrix" "[1][mRow-1]");
                        }
                        mRow--;
                        arrayRe[mRow][mRow-2]=#multAARe("ldiag" "[mRow]" "m.matrix" "[0][mRow-1]");
                        arrayIm[mRow][mRow-2]=#multAAIm("ldiag" "[mRow]" "m.matrix" "[0][mRow-1]");
                        arrayRe[mRow][mRow-1]=#multAARe("ldiag" "[mRow]" "m.matrix" "[1][mRow-1]")+#multAARe("diag" "[mRow]" "m.matrix" "[0][mRow]");
                        arrayIm[mRow][mRow-1]=#multAAIm("ldiag" "[mRow]" "m.matrix" "[1][mRow-1]")+#multAAIm("diag" "[mRow]" "m.matrix" "[0][mRow]");
                        arrayRe[mRow][mRow]=#multAARe("ldiag" "[mRow]" "m.matrix" "[2][mRow-1]")+#multAARe("diag" "[mRow]" "m.matrix" "[1][mRow]");
                        arrayIm[mRow][mRow]=#multAAIm("ldiag" "[mRow]" "m.matrix" "[2][mRow-1]")+#multAAIm("diag" "[mRow]" "m.matrix" "[1][mRow]");
                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                } else
                        throw new MatrixDimensionException("Incompatible matrices.");
        }

// HERMITIAN ADJOINT

        /**
        * Returns the hermitian adjoint of this matrix.
        * @return a complex tridiagonal matrix
        */
        public ComplexMatrix hermitianAdjoint() {
                int mRow=numRows;
                final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                System.arraycopy(ldiagRe,1,ans.udiagRe,0,ldiagRe.length-1);
                System.arraycopy(diagRe,0,ans.diagRe,0,diagRe.length);
                System.arraycopy(udiagRe,0,ans.ldiagRe,1,udiagRe.length-1);
                ans.diagIm[0]=-diagIm[0];
                ans.ldiagIm[1]=-udiagIm[0];
                mRow--;
                for(int i=1;i<mRow;i++) {
                        ans.udiagIm[i-1]=-ldiagIm[i];
                        ans.diagIm[i]=-diagIm[i];
                        ans.ldiagIm[i+1]=-udiagIm[i];
                }
                ans.udiagIm[mRow-1]=-ldiagIm[mRow];
                ans.diagIm[mRow]=-diagIm[mRow];
                return ans;
        }

// CONJUGATE

        /**
        * Returns the complex conjugate of this matrix.
        * @return a complex tridiagonal matrix
        */
        public ComplexMatrix conjugate() {
                int mRow=numRows;
                final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                System.arraycopy(ldiagRe,1,ans.ldiagRe,0,ldiagRe.length-1);
                System.arraycopy(diagRe,0,ans.diagRe,0,diagRe.length);
                System.arraycopy(udiagRe,0,ans.udiagRe,1,udiagRe.length-1);
                ans.diagIm[0]=-diagIm[0];
                ans.udiagIm[0]=-udiagIm[0];
                mRow--;
                for(int i=1;i<mRow;i++) {
                        ans.ldiagIm[i]=-ldiagIm[i];
                        ans.diagIm[i]=-diagIm[i];
                        ans.udiagIm[i]=-udiagIm[i];
                }
                ans.ldiagIm[mRow]=-ldiagIm[mRow];
                ans.diagIm[mRow]=-diagIm[mRow];
                return ans;
        }

// TRANSPOSE

        /**
        * Returns the transpose of this matrix.
        * @return a complex tridiagonal matrix
        */
        public Matrix transpose() {
                final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(numRows);
                System.arraycopy(ldiagRe,1,ans.udiagRe,0,ldiagRe.length-1);
                System.arraycopy(ldiagIm,1,ans.udiagIm,0,ldiagIm.length-1);
                System.arraycopy(diagRe,0,ans.diagRe,0,diagRe.length);
                System.arraycopy(diagIm,0,ans.diagIm,0,diagIm.length);
                System.arraycopy(udiagRe,0,ans.ldiagRe,1,udiagRe.length-1);
                System.arraycopy(udiagIm,0,ans.ldiagIm,1,udiagIm.length-1);
                return ans;
        }

// LU DECOMPOSITION

        /**
        * Returns the LU decomposition of this matrix.
        * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
        */
        public ComplexSquareMatrix[] luDecompose(int pivot[]) {
                if(LU!=null) {
                        if(pivot!=null)
                                System.arraycopy(LUpivot,0,pivot,0,pivot.length);
                        return LU;
                }
                int i,j,k,pivotrow;
                final int N=numRows;
                final double arrayLRe[][]=new double[N][N];
                final double arrayLIm[][]=new double[N][N];
                final double arrayURe[][]=new double[N][N];
                final double arrayUIm[][]=new double[N][N];
                final double buf[]=new double[N];
                double tmp,tmpRe,tmpIm;
                double max;
                Complex elem;
                if(pivot==null)
                        pivot=new int[N+1];
                for(i=0;i<N;i++)
                        pivot[i]=i;
                pivot[N]=1;
        // LU decomposition to arrayU
                for(j=0;j<N;j++) {
                        for(i=0;i<j;i++) {
                                elem=getElement(pivot[i],j);
                                tmpRe=elem.real();
                                tmpIm=elem.imag();
                                for(k=0;k<i;k++) {
                                        tmpRe -= #multAARe("arrayU" "[i][k]" "arrayU" "[k][j]");
                                        tmpIm -= #multAAIm("arrayU" "[i][k]" "arrayU" "[k][j]");
                                }
                                arrayURe[i][j]=tmpRe;
                                arrayUIm[i][j]=tmpIm;
                        }
                        max=0.0;
                        pivotrow=j;
                        for(i=j;i<N;i++) {
                                elem=getElement(pivot[i],j);
                                tmpRe=elem.real();
                                tmpIm=elem.imag();
                                for(k=0;k<j;k++) {
                                        tmpRe -= #multAARe("arrayU" "[i][k]" "arrayU" "[k][j]");
                                        tmpIm -= #multAAIm("arrayU" "[i][k]" "arrayU" "[k][j]");
                                }
                                arrayURe[i][j]=tmpRe;
                                arrayUIm[i][j]=tmpIm;
                        // while we're here search for a pivot for arrayU[j][j]
                                tmp=tmpRe*tmpRe+tmpIm*tmpIm;
                                if(tmp>max) {
                                        max=tmp;
                                        pivotrow=i;
                                }
                        }
                // swap row j with pivotrow
                        if(pivotrow!=j) {
                                System.arraycopy(arrayURe[j],0,buf,0,j+1);
                                System.arraycopy(arrayURe[pivotrow],0,arrayURe[j],0,j+1);
                                System.arraycopy(buf,0,arrayURe[pivotrow],0,j+1);
                                System.arraycopy(arrayUIm[j],0,buf,0,j+1);
                                System.arraycopy(arrayUIm[pivotrow],0,arrayUIm[j],0,j+1);
                                System.arraycopy(buf,0,arrayUIm[pivotrow],0,j+1);
                                k=pivot[j];
                                pivot[j]=pivot[pivotrow];
                                pivot[pivotrow]=k;
                                // update parity
                                pivot[N]=-pivot[N];
                        }
                // divide by pivot
                        tmpRe=arrayURe[j][j];
                        tmpIm=arrayUIm[j][j];
                        double a,denom;
                        if(Math.abs(tmpRe)<Math.abs(tmpIm)) {
                                a=tmpRe/tmpIm;
                                denom=tmpRe*a+tmpIm;
                                for(i=j+1;i<N;i++) {
                                        tmp=(arrayURe[i][j]*a+arrayUIm[i][j])/denom;
                                        arrayUIm[i][j]=(arrayUIm[i][j]*a-arrayURe[i][j])/denom;
                                        arrayURe[i][j]=tmp;
                                }
                        } else {
                                a=tmpIm/tmpRe;
                                denom=tmpRe+tmpIm*a;
                                for(i=j+1;i<N;i++) {
                                        tmp=(arrayURe[i][j]+arrayUIm[i][j]*a)/denom;
                                        arrayUIm[i][j]=(arrayUIm[i][j]-arrayURe[i][j]*a)/denom;
                                        arrayURe[i][j]=tmp;
                                }
                        }
                }
                // move lower triangular part to arrayL
                for(j=0;j<N;j++) {
                        arrayLRe[j][j]=1.0;
                        for(i=j+1;i<N;i++) {
                                arrayLRe[i][j]=arrayURe[i][j];
                                arrayLIm[i][j]=arrayUIm[i][j];
                                arrayURe[i][j]=0.0;
                                arrayUIm[i][j]=0.0;
                        }
                }
                LU=new ComplexSquareMatrix[2];
                LU[0]=new ComplexSquareMatrix(arrayLRe,arrayLIm);
                LU[1]=new ComplexSquareMatrix(arrayURe,arrayUIm);
                LUpivot=new int[pivot.length];
                System.arraycopy(pivot,0,LUpivot,0,pivot.length);
                return LU;
        }

// MAP ELEMENTS

        /**
        * Applies a function on all the matrix elements.
        * @param f a user-defined function
        * @return a complex tridiagonal matrix
        */
        public ComplexMatrix mapElements(final ComplexMapping f) {
                int mRow=numRows;
                final ComplexTridiagonalMatrix ans=new ComplexTridiagonalMatrix(mRow);
                ans.setElement(0,0,f.map(diagRe[0],diagIm[0]));
                ans.setElement(0,1,f.map(udiagRe[0],udiagIm[0]));
                mRow--;
                for(int i=1;i<mRow;i++) {
                        ans.setElement(i,i-1,f.map(ldiagRe[i],ldiagIm[i]));
                        ans.setElement(i,i,f.map(diagRe[i],diagIm[i]));
                        ans.setElement(i,i+1,f.map(udiagRe[i],udiagIm[i]));
                }
                ans.setElement(mRow,mRow-1,f.map(ldiagRe[mRow],ldiagIm[mRow]));
                ans.setElement(mRow,mRow,f.map(diagRe[mRow],diagIm[mRow]));
                return ans;
        }
}
