/* AUTO-GENERATED */
package jsci.maths.matrices;

import jsci.maths.ArrayMath;
import jsci.maths.Complex;
import jsci.maths.ComplexMapping;
import jsci.maths.DimensionException;
import jsci.maths.LinearMath;
import jsci.maths.MaximumIterationsExceededException;
import jsci.maths.vectors.AbstractComplexVector;
import jsci.maths.vectors.ComplexVector;

/**
* 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 AbstractComplexSquareMatrix implements TridiagonalMatrix {
        /**
        * Tridiagonal data.
        */
        protected final double ldiagRe[],ldiagIm[];
        protected final double diagRe[],diagIm[];
        protected final double udiagRe[],udiagIm[];
        /**
        * Constructs an empty matrix.
        * @param size the number of rows/columns
        */
        public ComplexTridiagonalMatrix(final int size) {
                super(size);
                ldiagRe=new double[size];
                ldiagIm=new double[size];
                diagRe=new double[size];
                diagIm=new double[size];
                udiagRe=new double[size];
                udiagIm=new double[size];
        }
        /**
        * 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(!ArrayMath.isSquare(array))
                        throw new MatrixDimensionException("Array is not square.");
                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=0 && i=0 && j=0 && i=0 && j=0 && i=0 && j=0 && i=0 && j=0 && i=0 && j-norm.
        * @author Taber Smith
        */
        public double infNorm() {
                double result = Math.sqrt((diagRe[0]*diagRe[0] + diagIm[0]*diagIm[0])) + Math.sqrt((udiagRe[0]*udiagRe[0] + udiagIm[0]*udiagIm[0]));
                double tmpResult;
                int i=1;
                for(;iresult)
                                result=tmpResult;
                }
                tmpResult = Math.sqrt((ldiagRe[i]*ldiagRe[i] + ldiagIm[i]*ldiagIm[i])) + Math.sqrt((diagRe[i]*diagRe[i] + diagIm[i]*diagIm[i]));
                if(tmpResult>result)
                        result=tmpResult;
                return result;
        }
        /**
        * Returns the Frobenius (l2) 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(;i3) {
                                arrayRe[1][0]=(ldiagRe[1]*m.diagRe[0] - ldiagIm[1]*m.diagIm[0])+(diagRe[1]*m.ldiagRe[1] - diagIm[1]*m.ldiagIm[1]);
                                arrayIm[1][0]=(ldiagIm[1]*m.diagRe[0] + ldiagRe[1]*m.diagIm[0])+(diagIm[1]*m.ldiagRe[1] + diagRe[1]*m.ldiagIm[1]);
                                arrayRe[1][1]=(ldiagRe[1]*m.udiagRe[0] - ldiagIm[1]*m.udiagIm[0])+(diagRe[1]*m.diagRe[1] - diagIm[1]*m.diagIm[1])+(udiagRe[1]*m.ldiagRe[2] - udiagIm[1]*m.ldiagIm[2]);
                                arrayIm[1][1]=(ldiagIm[1]*m.udiagRe[0] + ldiagRe[1]*m.udiagIm[0])+(diagIm[1]*m.diagRe[1] + diagRe[1]*m.diagIm[1])+(udiagIm[1]*m.ldiagRe[2] + udiagRe[1]*m.ldiagIm[2]);
                                arrayRe[1][2]=(diagRe[1]*m.udiagRe[1] - diagIm[1]*m.udiagIm[1])+(udiagRe[1]*m.diagRe[2] - udiagIm[1]*m.diagIm[2]);
                                arrayIm[1][2]=(diagIm[1]*m.udiagRe[1] + diagRe[1]*m.udiagIm[1])+(udiagIm[1]*m.diagRe[2] + udiagRe[1]*m.diagIm[2]);
                                arrayRe[1][3]=(udiagRe[1]*m.udiagRe[2] - udiagIm[1]*m.udiagIm[2]);
                                arrayIm[1][3]=(udiagIm[1]*m.udiagRe[2] + udiagRe[1]*m.udiagIm[2]);
                        }
                        if(mRow==3) {
                                arrayRe[1][0]=(ldiagRe[1]*m.diagRe[0] - ldiagIm[1]*m.diagIm[0])+(diagRe[1]*m.ldiagRe[1] - diagIm[1]*m.ldiagIm[1]);
                                arrayIm[1][0]=(ldiagIm[1]*m.diagRe[0] + ldiagRe[1]*m.diagIm[0])+(diagIm[1]*m.ldiagRe[1] + diagRe[1]*m.ldiagIm[1]);
                                arrayRe[1][1]=(ldiagRe[1]*m.udiagRe[0] - ldiagIm[1]*m.udiagIm[0])+(diagRe[1]*m.diagRe[1] - diagIm[1]*m.diagIm[1])+(udiagRe[1]*m.ldiagRe[2] - udiagIm[1]*m.ldiagIm[2]);
                                arrayIm[1][1]=(ldiagIm[1]*m.udiagRe[0] + ldiagRe[1]*m.udiagIm[0])+(diagIm[1]*m.diagRe[1] + diagRe[1]*m.diagIm[1])+(udiagIm[1]*m.ldiagRe[2] + udiagRe[1]*m.ldiagIm[2]);
                                arrayRe[1][2]=(diagRe[1]*m.udiagRe[1] - diagIm[1]*m.udiagIm[1])+(udiagRe[1]*m.diagRe[2] - udiagIm[1]*m.diagIm[2]);
                                arrayIm[1][2]=(diagIm[1]*m.udiagRe[1] + diagRe[1]*m.udiagIm[1])+(udiagIm[1]*m.diagRe[2] + udiagRe[1]*m.diagIm[2]);
                        } else if(mRow>4) {
                                for(int i=2;i3) {
                                arrayRe[mRow-2][mRow-4]=(ldiagRe[mRow-2]*m.ldiagRe[mRow-3] - ldiagIm[mRow-2]*m.ldiagIm[mRow-3]);
                                arrayIm[mRow-2][mRow-4]=(ldiagIm[mRow-2]*m.ldiagRe[mRow-3] + ldiagRe[mRow-2]*m.ldiagIm[mRow-3]);
                                arrayRe[mRow-2][mRow-3]=(ldiagRe[mRow-2]*m.diagRe[mRow-3] - ldiagIm[mRow-2]*m.diagIm[mRow-3])+(diagRe[mRow-2]*m.ldiagRe[mRow-2] - diagIm[mRow-2]*m.ldiagIm[mRow-2]);
                                arrayIm[mRow-2][mRow-3]=(ldiagIm[mRow-2]*m.diagRe[mRow-3] + ldiagRe[mRow-2]*m.diagIm[mRow-3])+(diagIm[mRow-2]*m.ldiagRe[mRow-2] + diagRe[mRow-2]*m.ldiagIm[mRow-2]);
                                arrayRe[mRow-2][mRow-2]=(ldiagRe[mRow-2]*m.udiagRe[mRow-3] - ldiagIm[mRow-2]*m.udiagIm[mRow-3])+(diagRe[mRow-2]*m.diagRe[mRow-2] - diagIm[mRow-2]*m.diagIm[mRow-2])+(udiagRe[mRow-2]*m.ldiagRe[mRow-1] - udiagIm[mRow-2]*m.ldiagIm[mRow-1]);
                                arrayIm[mRow-2][mRow-2]=(ldiagIm[mRow-2]*m.udiagRe[mRow-3] + ldiagRe[mRow-2]*m.udiagIm[mRow-3])+(diagIm[mRow-2]*m.diagRe[mRow-2] + diagRe[mRow-2]*m.diagIm[mRow-2])+(udiagIm[mRow-2]*m.ldiagRe[mRow-1] + udiagRe[mRow-2]*m.ldiagIm[mRow-1]);
                                arrayRe[mRow-2][mRow-1]=(diagRe[mRow-2]*m.udiagRe[mRow-2] - diagIm[mRow-2]*m.udiagIm[mRow-2])+(udiagRe[mRow-2]*m.diagRe[mRow-1] - udiagIm[mRow-2]*m.diagIm[mRow-1]);
                                arrayIm[mRow-2][mRow-1]=(diagIm[mRow-2]*m.udiagRe[mRow-2] + diagRe[mRow-2]*m.udiagIm[mRow-2])+(udiagIm[mRow-2]*m.diagRe[mRow-1] + udiagRe[mRow-2]*m.diagIm[mRow-1]);
                        }
                        mRow--;
                        arrayRe[mRow][mRow-2]=(ldiagRe[mRow]*m.ldiagRe[mRow-1] - ldiagIm[mRow]*m.ldiagIm[mRow-1]);
                        arrayIm[mRow][mRow-2]=(ldiagIm[mRow]*m.ldiagRe[mRow-1] + ldiagRe[mRow]*m.ldiagIm[mRow-1]);
                        arrayRe[mRow][mRow-1]=(ldiagRe[mRow]*m.diagRe[mRow-1] - ldiagIm[mRow]*m.diagIm[mRow-1])+(diagRe[mRow]*m.ldiagRe[mRow] - diagIm[mRow]*m.ldiagIm[mRow]);
                        arrayIm[mRow][mRow-1]=(ldiagIm[mRow]*m.diagRe[mRow-1] + ldiagRe[mRow]*m.diagIm[mRow-1])+(diagIm[mRow]*m.ldiagRe[mRow] + diagRe[mRow]*m.ldiagIm[mRow]);
                        arrayRe[mRow][mRow]=(ldiagRe[mRow]*m.udiagRe[mRow-1] - ldiagIm[mRow]*m.udiagIm[mRow-1])+(diagRe[mRow]*m.diagRe[mRow] - diagIm[mRow]*m.diagIm[mRow]);
                        arrayIm[mRow][mRow]=(ldiagIm[mRow]*m.udiagRe[mRow-1] + ldiagRe[mRow]*m.udiagIm[mRow-1])+(diagIm[mRow]*m.diagRe[mRow] + diagRe[mRow]*m.diagIm[mRow]);
                        return new ComplexSquareMatrix(arrayRe,arrayIm);
                } else
                        throw new MatrixDimensionException("Incompatible matrices.");
        }
        private ComplexSquareMatrix multiplyTridiagonal(final AbstractComplexSquareMatrix m) {
                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]=(diagRe[0]*elem1.real() - diagIm[0]*elem1.imag())+(udiagRe[0]*elem2.real() - udiagIm[0]*elem2.imag());
                        arrayIm[0][0]=(diagIm[0]*elem1.real() + diagRe[0]*elem1.imag())+(udiagIm[0]*elem2.real() + udiagRe[0]*elem2.imag());
                        elem1=m.getElement(0,1);elem2=m.getElement(1,1);
                        arrayRe[0][1]=(diagRe[0]*elem1.real() - diagIm[0]*elem1.imag())+(udiagRe[0]*elem2.real() - udiagIm[0]*elem2.imag());
                        arrayIm[0][1]=(diagIm[0]*elem1.real() + diagRe[0]*elem1.imag())+(udiagIm[0]*elem2.real() + udiagRe[0]*elem2.imag());
                        elem1=m.getElement(1,2);
                        arrayRe[0][2]=(udiagRe[0]*elem1.real() - udiagIm[0]*elem1.imag());
                        arrayIm[0][2]=(udiagIm[0]*elem1.real() + udiagRe[0]*elem1.imag());
                        if(mRow>3) {
                                elem1=m.getElement(0,0);elem2=m.getElement(1,0);
                                arrayRe[1][0]=(ldiagRe[1]*elem1.real() - ldiagIm[1]*elem1.imag())+(diagRe[1]*elem2.real() - diagIm[1]*elem2.imag());
                                arrayIm[1][0]=(ldiagIm[1]*elem1.real() + ldiagRe[1]*elem1.imag())+(diagIm[1]*elem2.real() + diagRe[1]*elem2.imag());
                                elem1=m.getElement(0,1);elem2=m.getElement(1,1);elem3=m.getElement(2,1);
                                arrayRe[1][1]=(ldiagRe[1]*elem1.real() - ldiagIm[1]*elem1.imag())+(diagRe[1]*elem2.real() - diagIm[1]*elem2.imag())+(udiagRe[1]*elem3.real() - udiagIm[1]*elem3.imag());
                                arrayIm[1][1]=(ldiagIm[1]*elem1.real() + ldiagRe[1]*elem1.imag())+(diagIm[1]*elem2.real() + diagRe[1]*elem2.imag())+(udiagIm[1]*elem3.real() + udiagRe[1]*elem3.imag());
                                elem1=m.getElement(1,2);elem2=m.getElement(2,2);
                                arrayRe[1][2]=(diagRe[1]*elem1.real() - diagIm[1]*elem1.imag())+(udiagRe[1]*elem2.real() - udiagIm[1]*elem2.imag());
                                arrayIm[1][2]=(diagIm[1]*elem1.real() + diagRe[1]*elem1.imag())+(udiagIm[1]*elem2.real() + udiagRe[1]*elem2.imag());
                                elem1=m.getElement(2,3);
                                arrayRe[1][3]=(udiagRe[1]*elem1.real() - udiagIm[1]*elem1.imag());
                                arrayIm[1][3]=(udiagIm[1]*elem1.real() + udiagRe[1]*elem1.imag());
                        }
                        if(mRow==3) {
                                elem1=m.getElement(0,0);elem2=m.getElement(1,0);
                                arrayRe[1][0]=(ldiagRe[1]*elem1.real() - ldiagIm[1]*elem1.imag())+(diagRe[1]*elem2.real() - diagIm[1]*elem2.imag());
                                arrayIm[1][0]=(ldiagIm[1]*elem1.real() + ldiagRe[1]*elem1.imag())+(diagIm[1]*elem2.real() + diagRe[1]*elem2.imag());
                                elem1=m.getElement(0,1);elem2=m.getElement(1,1);elem3=m.getElement(2,1);
                                arrayRe[1][1]=(ldiagRe[1]*elem1.real() - ldiagIm[1]*elem1.imag())+(diagRe[1]*elem2.real() - diagIm[1]*elem2.imag())+(udiagRe[1]*elem3.real() - udiagIm[1]*elem3.imag());
                                arrayIm[1][1]=(ldiagIm[1]*elem1.real() + ldiagRe[1]*elem1.imag())+(diagIm[1]*elem2.real() + diagRe[1]*elem2.imag())+(udiagIm[1]*elem3.real() + udiagRe[1]*elem3.imag());
                                elem1=m.getElement(1,2);elem2=m.getElement(2,2);
                                arrayRe[1][2]=(diagRe[1]*elem1.real() - diagIm[1]*elem1.imag())+(udiagRe[1]*elem2.real() - udiagIm[1]*elem2.imag());
                                arrayIm[1][2]=(diagIm[1]*elem1.real() + diagRe[1]*elem1.imag())+(udiagIm[1]*elem2.real() + udiagRe[1]*elem2.imag());
                        } else if(mRow>4) {
                                for(int i=2;i3) {
                                elem1=m.getElement(mRow-3,mRow-4);
                                arrayRe[mRow-2][mRow-4]=(ldiagRe[mRow-2]*elem1.real() - ldiagIm[mRow-2]*elem1.imag());
                                arrayIm[mRow-2][mRow-4]=(ldiagIm[mRow-2]*elem1.real() + ldiagRe[mRow-2]*elem1.imag());
                                elem1=m.getElement(mRow-3,mRow-3);elem2=m.getElement(mRow-2,mRow-3);
                                arrayRe[mRow-2][mRow-3]=(ldiagRe[mRow-2]*elem1.real() - ldiagIm[mRow-2]*elem1.imag())+(diagRe[mRow-2]*elem2.real() - diagIm[mRow-2]*elem2.imag());
                                arrayIm[mRow-2][mRow-3]=(ldiagIm[mRow-2]*elem1.real() + ldiagRe[mRow-2]*elem1.imag())+(diagIm[mRow-2]*elem2.real() + diagRe[mRow-2]*elem2.imag());
                                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]=(ldiagRe[mRow-2]*elem1.real() - ldiagIm[mRow-2]*elem1.imag())+(diagRe[mRow-2]*elem2.real() - diagIm[mRow-2]*elem2.imag())+(udiagRe[mRow-2]*elem3.real() - udiagIm[mRow-2]*elem3.imag());
                                arrayIm[mRow-2][mRow-2]=(ldiagIm[mRow-2]*elem1.real() + ldiagRe[mRow-2]*elem1.imag())+(diagIm[mRow-2]*elem2.real() + diagRe[mRow-2]*elem2.imag())+(udiagIm[mRow-2]*elem3.real() + udiagRe[mRow-2]*elem3.imag());
                                elem1=m.getElement(mRow-2,mRow-1);elem2=m.getElement(mRow-1,mRow-1);
                                arrayRe[mRow-2][mRow-1]=(diagRe[mRow-2]*elem1.real() - diagIm[mRow-2]*elem1.imag())+(udiagRe[mRow-2]*elem2.real() - udiagIm[mRow-2]*elem2.imag());
                                arrayIm[mRow-2][mRow-1]=(diagIm[mRow-2]*elem1.real() + diagRe[mRow-2]*elem1.imag())+(udiagIm[mRow-2]*elem2.real() + udiagRe[mRow-2]*elem2.imag());
                        }
                        mRow--;
                        elem1=m.getElement(mRow-1,mRow-2);
                        arrayRe[mRow][mRow-2]=(ldiagRe[mRow]*elem1.real() - ldiagIm[mRow]*elem1.imag());
                        arrayIm[mRow][mRow-2]=(ldiagIm[mRow]*elem1.real() + ldiagRe[mRow]*elem1.imag());
                        elem1=m.getElement(mRow-1,mRow-1);elem2=m.getElement(mRow,mRow-1);
                        arrayRe[mRow][mRow-1]=(ldiagRe[mRow]*elem1.real() - ldiagIm[mRow]*elem1.imag())+(diagRe[mRow]*elem2.real() - diagIm[mRow]*elem2.imag());
                        arrayIm[mRow][mRow-1]=(ldiagIm[mRow]*elem1.real() + ldiagRe[mRow]*elem1.imag())+(diagIm[mRow]*elem2.real() + diagRe[mRow]*elem2.imag());
                        elem1=m.getElement(mRow-1,mRow);elem2=m.getElement(mRow,mRow);
                        arrayRe[mRow][mRow]=(ldiagRe[mRow]*elem1.real() - ldiagIm[mRow]*elem1.imag())+(diagRe[mRow]*elem2.real() - diagIm[mRow]*elem2.imag());
                        arrayIm[mRow][mRow]=(ldiagIm[mRow]*elem1.real() + ldiagRe[mRow]*elem1.imag())+(diagIm[mRow]*elem2.real() + diagRe[mRow]*elem2.imag());
                        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 AbstractComplexMatrix 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