// ComplexMatrix.java
// solve, invert, etc. last argument is output
// void solve(Complex A[][], Complex Y[], Complex X[]);        X=A^-1*Y
// void invert(Complex A[][]);                                 A=A^-1
// Complex determinant(Complex A[]);                           d=det A
// void eigenvalues(Complex A[][], ComplexV[][], Complex Y[]); V,Y=eigen A
// void eigenCheck(Complex A[][], V[][], Complex Y[]);         printout
// void multiply(Complex A[][], Complex B[][], Complex C[][]); C=A*B
// void add(Complex A[][], Complex B[][], Complex C[][]);      C=A+B
// void subtract(Complex C[][], Complex A[][], complex B[][]); C=A-B
// double norm1(Complex A[][]);                                d=norm1 A
// double norm2(Complex A[][]);  sqrt largest eigenvalue A^T*A d=norm2 A
// double normFro(Complex A[][]); Frobenius                    d=normFro A
// double normInf(Complex A[][]);                              d=normInf A
// void copy(Complex A[][], Complex B[][]);                    B=A
// void boolean equals(Complex A[][], Complex B[][]);          A==B
// void fromDouble(double A[], double b[][], Complex C[][]);   C=(A,b)
// void fromDouble(double A[], Complex B[][]);                 B=A
// void identity(A[][]);                                       A=I
// void zero(A[][]);                                           A=0
// void print(Complex A[][]);                                  A
// void multiply(Complex A[][], Complex X[], Complex Y[]);     Y=A*X
// void add(Complex X[], Complex Y[], Complex Z[]);            Z=X+Y
// void subtract(Complex X[], Complex Y[], Complex Z[]);       Z=X-Y
// double norm1(Complex X[]);                                  d=norm1 X
// double norm2(Complex X[]);                                  d=norm2 X
// double normInf(Complex X[]);                                d=normInf X
// void copy(Complex X[], Complex Y[]);                        Y=X
// void boolean equals(Complex X[], Complex Y[]);              X==Y
// void fromDouble(double X[], double Y[], Complex Z[]);       Z=(X,Y)
// void fromDouble(double X[], Complex Y[]);                   Y=(X,0)
// void unitVector(X[],I);                                     X[I]=1, else 0
// void zero(X[]);                                             X=0
// void print(Complex X[]);                                    X
// void readSize(String, rowCol[])
// void read(String, Complex[][])
// void closeInput()
// void write(String, Complex[][])
// void closeOutput()

package myjava;

import myjava.*;
import java.io.*;

public strictfp class ComplexMatrix
{

  public static void solve(final Complex A[][], final Complex Y[],
                           Complex X[])
  {
    // solve complex linear equations for X where Y = A * X
    // method: Gauss-Jordan elimination using maximum pivot
    // usage:  ComplexMatrix.solve(A,Y,X);
    //    Translated to java by : Jon Squire , 4 April 2003
    //    First written by Jon Squire December 1959 for IBM 650, translated to
    //    other languages  e.g. Fortran converted to Ada converted to C
    //    converted to java
    int n=A.length;
    int m=n+1;
    Complex B[][]=new Complex[n][m];   // working matrix
    int row[]=new int[n];              // row interchange indicies
    int hold , I_pivot;                // pivot indicies
    Complex pivot;                     // pivot element value
    double abs_pivot;
    if(A[0].length!=n || Y.length!=n || X.length!=n)
    {
      System.out.println("Error in ComplexMatrix.solve,"+
                         " inconsistent array sizes.");
    }
    // build working data structure
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        B[i][j] = A[i][j];
      }
      B[i][n] = Y[i];
    }
    // set up row interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k] = k;
    }
    //  begin main reduction loop
    for(int k=0; k<n; k++)
    {
      // find largest element for pivot
      pivot = B[row[k]][k];
      abs_pivot = pivot.abs();
      I_pivot = k;
      for(int i=k+1; i<n; i++)
      {
        if(B[row[i]][k].abs() > abs_pivot)
        {
          I_pivot = i;
          pivot = B[row[i]][k];
          abs_pivot = pivot.abs();
        }
      }
      // have pivot, interchange row indicies
      hold = row[k];
      row[k] = row[I_pivot];
      row[I_pivot] = hold;
      // check for near singular
      if(abs_pivot < 1.0E-10)
      {
        for(int j=k+1; j<n+1; j++)
        {
          B[row[k]][j] = new Complex(0.0, 0.0);
        }
        System.out.println("redundant row (singular) "+row[k]);
      } // singular, delete row
      else
      {
        // reduce about pivot
        for(int j=k+1; j<n+1; j++)
        {
          B[row[k]][j] = B[row[k]][j].divide(B[row[k]][k]);
        }
        //  inner reduction loop
        for(int i=0; i<n; i++)
        {
          if(i != k)
          {
            for(int j=k+1; j<n+1; j++)
            {
              B[row[i]][j] = B[row[i]][j].subtract(
                             B[row[i]][k].multiply(B[row[k]][j]));
            }
          }
        }
      }
      //  finished inner reduction
    }
    //  end main reduction loop
    //  build  X  for return, unscrambling rows
    for(int i=0; i<n; i++)
    {
      X[i] = B[row[i]][n];
    }
  } // end solve

  public static final void invert(Complex A[][])
  {
    int n = A.length;
    int row[] = new int[n];
    int col[] = new int[n];
    Complex temp[] = new Complex[n];
    int hold , I_pivot , J_pivot;
    Complex pivot;
    double  abs_pivot;

    if(A[0].length!=n)
    {
      System.out.println("Error in Complex.Matrix.invert,"+
                         " matrix not square.");
    }
    // set up row and column interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k] = k ;
      col[k] = k ;
    }
    // begin main reduction loop
    for(int k=0; k<n; k++)
    {
      // find largest element for pivot
      pivot = A[row[k]][col[k]];
      I_pivot = k;
      J_pivot = k;
      for(int i=k; i<n; i++)
      {
        for(int j=k; j<n; j++)
        {
          abs_pivot = pivot.abs();
          if(A[row[i]][col[j]].abs() > abs_pivot)
          {
            I_pivot = i ;
            J_pivot = j ;
            pivot = A[row[i]][col[j]] ;
          }
        }
      }
      if(pivot.abs() < 1.0E-10)
      {
        System.out.println("ComplexMatrix is singular !");
        return;
      }
      hold = row[k];
      row[k]= row[I_pivot];
      row[I_pivot] = hold ;
      hold = col[k];
      col[k]= col[J_pivot];
      col[J_pivot] = hold ;
      // reduce about pivot
      A[row[k]][col[k]] = (new Complex(1.0,0.0)).divide(pivot) ;
      for(int j=0; j<n; j++)
      {
        if(j != k)
        {
          A[row[k]][col[j]] = A[row[k]][col[j]].multiply(A[row[k]][col[k]]);
        }
      }
      // inner reduction loop
      for(int i=0; i<n; i++)
      {
        if(k != i)
        {
          for(int j=0; j<n; j++)
          {
            if( k != j )
            {
              A[row[i]][col[j]] = A[row[i]][col[j]].subtract(
                                  A[row[i]][col[k]].multiply(A[row[k]][col[j]]));
            }
          }
          A[row[i]][col [k]] = A[row[i]][col[k]].multiply(
                               A[row[k]][col[k]]).negate();
        }
      }
    }
    // end main reduction loop

    // unscramble rows
    for(int j=0; j<n; j++)
    {
      for(int i=0; i<n; i++)
      {
        temp[col[i]] = A[row[i]][j];
      }
      for(int i=0; i<n; i++)
      {
        A[i][j] = temp[i] ;
      }
    }
    // unscramble columns
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        temp[row[j]] = A[i][col[j]] ;
      }
      for(int j=0; j<n; j++)
      {
        A[i][j] = temp[j] ;
      }
    }
  } // end invert

  public static final Complex determinant(final Complex A[][])
  {
    int n=A.length;
    Complex D = new Complex(1.0,0.0);  // determinant
    Complex B[][]=new Complex[n][n];   // working matrix
    int row[]=new int[n];              // row interchange indicies
    int hold , I_pivot;                // pivot indicies
    Complex pivot;                     // pivot element value
    double abs_pivot;

    if(A[0].length!=n)
    {
      System.out.println("Error in ComplexMatrix.determinant,"+
                         " inconsistent array sizes.");
    }
    // build working matrix
    for(int i=0; i<n; i++)
      for(int j=0; j<n; j++)
        B[i][j]=A[i][j];
    // set up row interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k]= k;
    }
    // begin main reduction loop
    for(int k=0; k<n-1; k++)
    {
      // find largest element for pivot
      pivot = B[row[k]][k];
      abs_pivot = pivot.abs();
      I_pivot = k;
      for(int i=k; i<n; i++)
      {
        if(B[row[i]][k].abs() > abs_pivot )
        {
          I_pivot = i;
          pivot = B[row[i]][k];
          abs_pivot = pivot.abs();
        }
      }
      // have pivot, interchange row indicies
      if(I_pivot != k)
      {
        hold = row[k];
        row[k] = row[I_pivot];
        row[I_pivot] = hold;
        D = D.negate();
      }
      // check for near singular
      if(abs_pivot < 1.0E-10)
      {
        return new Complex(0.0,0.0);
      }
      else
      {
        D = D.multiply(pivot);
        // reduce about pivot
        for(int j=k+1; j<n; j++)
        {
          B[row[k]][j] = B[row[k]][j].divide(B[row[k]][k]);
        }
        //  inner reduction loop
        for(int i=0; i<n; i++)
        {
          if(i != k)
          {
            for(int j=k+1; j<n; j++)
            {
              B[row[i]][j] = B[row[i]][j].subtract(
                             B[row[i]][k].multiply(B[row[k]][j]));
            }
          }
        }
      }
      //  finished inner reduction
    }
    // end of main reduction loop
    return D.multiply(B[row[n-1]][n-1]);
  } // end determinant

  public static final void eigenvalues(final Complex A[][], Complex V[][],
                                       Complex Y[])
  {
    // cyclic Jacobi iterative method of finding eigenvalues
    // advertized for symmetric real
    int n=A.length;
    Complex AA[][] = new Complex[n][n];
    double norm;
    Complex c= new Complex(1.0,0.0);
    Complex s= new Complex(0.0,0.0);
    if(A[0].length!=n || V.length!=n || V[0].length!=n || Y.length!=n)
    {
      System.out.println("Error in ComplexMatrix.eigenvalues,"+
                         " inconsistent array sizes.");
    }
    identity(V); // start V as identity matrix
    copy(A, AA);
    for(int k=0; k<n; k++)
    {
      norm=norm4(AA);
      for(int i=0; i<n-1; i++)
      {
        for(int j=i+1; j<n; j++)
        {
          schur2(AA, i, j, c, s);
          mat44(i, j, c, s, AA, V);
        }
      } // end one iteration
    }
    norm = norm4(AA); // final quality check if desired
    for(int i=0; i<n; i++) // copy eigenvalues back to caller
      Y[i] = AA[i][i];
  } // end eigenvalues

  public static final void eigenCheck(final Complex A[][],
                                      final Complex V[][], final Complex Y[])
  {
    if(A==null || V==null || Y==null) return;
    // check  A * X = lambda X  lambda=Y[i] X=V[i]
    // check  determinant(A- lambda I) = 0
    int n=A.length;
    Complex B[][] = new Complex[n][n];
    Complex C[][] = new Complex[n][n];
    Complex X[] = new Complex[n];
    Complex Z[] = new Complex[n];
    Complex T[] = new Complex[n];
    double norm=0.0;

    if(A[0].length!=n || V.length!=n || V[0].length!=n || Y.length!=n)
    {
      System.out.println("Error in ComplexMatrix.eigenCheck,"+
                         " inconsistent array sizes.");
    }
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        X[j]=V[j][i];
      }
      multiply(A, X, T);
      for(int j=0; j<n; j++)
      {
        Z[j]=T[j].subtract(Y[i].multiply(X[j]));
      }
      System.out.println("check for near zero norm of Z["+i+"]="+Z[i]);
    }
    norm = norm2(Z);
    System.out.println("norm ="+norm+" is eigen vector error indication 1.");
    System.out.println("det V = "+ComplexMatrix.determinant(V));
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        Z[j]=V[j][i];
      }
      System.out.println("check for 1.0 = "+norm2(Z));
    }
    
    for(int i=0; i<n; i++)
    {
      identity(B);
      multiply(B,Y[i],C);
      subtract(A,C,B);
      Z[i]=determinant(B);
    }
    norm = norm2(Z);
    System.out.println("norm ="+norm+" is eigen value error indication.");
  } // end eigenCheck

  static void schur2(final Complex A[][], final int p, final int q,
                     Complex c, Complex s)
  {
    Complex tau;
    Complex tau_tau_1;
    Complex t;

    if(A[0].length!=A.length)
    {
      System.out.println("Error in schur2 of Complex jacobi,"+
                         " inconsistent array sizes.");
    }
    if(A[p][q].abs()!=0.0)
    {
      tau=(A[q][q].subtract(A[p][p])).divide((A[p][q].multiply(2.0)));
      tau_tau_1=(tau.multiply(tau).add(1.0)).sqrt();

      if(tau.abs()>=0.0)
         t=tau.add(tau_tau_1).invert();
      else
         t=(tau_tau_1.subtract(tau)).invert().negate();
      c=(t.multiply(t)).add(1.0).sqrt().invert();
      s=t.multiply(c);
    }
    else
    {
      c=new Complex(1.0,0.0);
      s=new Complex(0.0,0.0);
    }
  } // end schur2

  static void mat22(final Complex c, final Complex s, final Complex A[][],
                    Complex B[][])
  {
    if(A.length!=2 || A[0].length!=2 || B.length!=2 || B[0].length!=2)
    {
      System.out.println("Error in mat22 of Jacobi, not both 2 by 2");
    }
    Complex T[][] = new Complex[2][2];

    T[0][0] = c.multiply(A[0][0]).subtract(s.multiply(A[0][1]));
    T[0][1] = s.multiply(A[0][0]).add     (c.multiply(A[0][1]));
    T[1][0] = c.multiply(A[1][0]).subtract(s.multiply(A[1][1]));
    T[1][1] = s.multiply(A[1][0]).add     (c.multiply(A[1][1]));

    B[0][0] = c.multiply(T[0][0]).subtract(s.multiply(T[1][0]));
    B[0][1] = c.multiply(T[0][1]).subtract(s.multiply(T[1][1]));
    B[1][0] = s.multiply(T[0][0]).add     (c.multiply(T[1][0]));
    B[1][1] = s.multiply(T[0][1]).add     (c.multiply(T[1][1]));
  } // end mat2

  static void mat44(final int p, final int q, final Complex c,
                    final Complex s, final Complex A[][], Complex V[][])
  {
    int n = A.length;
    Complex B[][] = new Complex[n][n];
    Complex J[][] = new Complex[n][n];

    if(A[0].length!=n || V.length!=n || V[0].length!=n)
    {
      System.out.println("Error in mat44 of Complex Jacobi,"+
                         " A or V not same and square");
    }
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        J[i][j]=new Complex(0.0,0.0);
      }
      J[i][i]=new Complex(1.0,0.0);
    }
    J[p][p]=c; /* J transpose */
    J[p][q]=s.negate();
    J[q][q]=c;
    J[q][p]=s;
    multiply(J, A, B);
    J[p][q]=s;
    J[q][p]=s.negate();
    multiply(B, J, A);
    multiply(V, J, B);
    copy(B, V);
  } // end mat44

  static double norm4(final Complex A[][]) // for Jacobi
  {
    int n=A.length;
    int nr=A[0].length;
    double nrm=0.0;
    if(n!=nr)
    {
      System.out.println("Error in Complex norm4, non square A["+
                         n+"]["+nr+"]");
    }
    for(int i=0; i<n-1; i++)
    {
      for(int j=i+1; j<n; j++)
      {
        nrm=nrm+A[i][j].abs()+A[j][i].abs();
      }
    }
    return nrm/(n*n-n);
  } // end norm4

  public static final void multiply(final Complex A[][], final Complex B[][],
                                    Complex C[][])
  {
    int ni = A.length;
    int nk = A[0].length;
    int nj = B[0].length;
    if(B.length!=nk || C.length!=ni || C[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.multiply,"+
                         " incompatible sizes");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
      {
        C[i][j] = new Complex(0.0, 0.0);
        for(int k=0; k<nk; k++)
          C[i][j] = C[i][j].add(A[i][k].multiply(B[k][j]));
      }
  } // end multiply

  public static final void multiply(final Complex A[][], final Complex B,
                                    Complex C[][])
  {
    int ni = A.length;
    int nj = A[0].length;
    if(C.length!=ni || C[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.multiply,"+
                         " incompatible sizes");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        C[i][j] = A[i][j].multiply(B);
  } // end multiply

  public static final void add(final Complex A[][], final Complex B[][],
                               Complex C[][])
  {
    int ni=A.length;
    int nj=A[0].length;
    if(B.length!=ni || C.length!=ni || B[0].length!=nj || C[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.add,"+
                         " incompatible sizes");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        C[i][j] = A[i][j].add(B[i][j]);
  } // end add

  public static final void add(final Complex A[][], final Complex B,
                               Complex C[][])
  {
    int ni=A.length;
    int nj=A[0].length;
    if(C.length!=ni || C[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.add,"+
                         " incompatible sizes");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        C[i][j] = A[i][j].add(B);
  } // end add

  public static final void subtract(final Complex A[][], final Complex B[][],
                                    Complex C[][])
  {
    int ni=A.length;
    int nj=A[0].length;
    if(B.length!=ni || C.length!=ni || B[0].length!=nj || C[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.subtract,"+
                         " incompatible sizes");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        C[i][j] = A[i][j].subtract(B[i][j]);
  } // end subtract

  public static final void subtract(final Complex A[][], final Complex B,
                                    Complex C[][])
  {
    int ni=A.length;
    int nj=A[0].length;
    if(C.length!=ni || C[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.subtract,"+
                         " incompatible sizes");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        C[i][j] = A[i][j].subtract(B);
  } // end subtract

  public static final double norm1(final Complex A[][])
  {
    double norm=0.0;
    double colSum;
    int ni=A.length;
    int nj=A[0].length;
    for(int j=0; j<nj; j++)
    {
      colSum = 0.0;
      for(int i=0; i<ni; i++)
        colSum = colSum + A[i][j].abs();
      norm = Math.max(norm, colSum);
    }
    return norm;
  } // end norm1

  public static final double normInf(final Complex A[][])
  {
    double norm=0.0;
    double rowSum;
    int ni=A.length;
    int nj=A[0].length;
    for(int i=0; i<ni; i++)
    {
      rowSum = 0.0;
      for(int j=0; j<nj; j++)
        rowSum = rowSum + A[i][j].abs();
      norm = Math.max(norm, rowSum);
    }
    return norm;
  } // end normInf

  public static final double normFro(final Complex A[][])
  {
    double norm=0.0;
    int n=A.length;
    for(int i=0; i<n; i++)
      for(int j=0; j<n; j++)
        norm = norm + A[i][j].abs()*A[i][j].abs();
    return Math.sqrt(norm);
  } // end normFro


  public static final double norm2(final Complex A[][])
  {
    double r=0.0; // largest eigenvalue
    int n=A.length;
    Complex B[][] = new Complex[n][n];
    Complex V[][] = new Complex[n][n];
    Complex BI[] = new Complex[n];
    if(A[0].length!=n)
    {
      System.out.println("Error in ComplexMatrix.norm2,"+
                         " matrix not square.");
    }
    for(int i=0; i<n; i++)  // B = A^T * A
    {
      for(int j=0; j<n; j++)
      {
        B[i][j]=new Complex(0.0,0.0);
        for(int k=0; k<n; k++)
          B[i][j] = B[i][j].add(A[k][i].multiply(A[k][j]));
      }
    }
    eigenvalues(B, V, BI);
    for(int i=0; i<n; i++) r=Math.max(r,BI[i].abs());
    return Math.sqrt(r);
  } // end subtract

  public static final void copy(final Complex A[][], Complex B[][])
  {
    int ni = A.length;
    int nj = A[0].length;
    if(B.length!=ni || B[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.copy,"+
                         " inconsistent sizes.");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        B[i][j] = A[i][j];
  } // end copy

  public static final boolean equals(final Complex A[][], final Complex B[][])
  {
    int ni = A.length;
    int nj = A[0].length;
    boolean same = true;
    if(B.length!=ni || B[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.equals,"+
                         " inconsistent sizes.");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        same = same && A[i][j].equals(B[i][j]);
    return same;
  } // end equals

  public static final void fromDouble(final double A[][],
                                      final double B[][], Complex C[][])
  {
    int ni = A.length;
    int nj = A[0].length;
    if(C.length!=ni || C[0].length!=nj || B.length!=ni || B[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.fromDouble,"+
                         " inconsistent sizes.");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        C[i][j] = new Complex(A[i][j],B[i][j]);
  } // end fromDouble

  public static final void fromDouble(final double A[][], Complex C[][])
  {
    int ni = A.length;
    int nj = A[0].length;
    if(C.length!=ni || C[0].length!=nj)
    {
      System.out.println("Error in ComplexMatrix.fromDouble,"+
                         " inconsistent sizes.");
    }
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        C[i][j] = new Complex(A[i][j]);
  } // end fromDouble

  public static final void identity(Complex A[][])
  {
    int n = A.length;
    if(n!=A[0].length)
    {
      System.out.println("Error in ComplexMatrix.identity,"+
                         " inconsistent sizes.");
    }
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
        A[i][j]=new Complex(0.0);
      A[i][i]=new Complex(1.0);
    }
  } // end identity

  public static final void zero(Complex A[][])
  {
    int ni = A.length;
    int nj = A[0].length;
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        A[i][j]=new Complex(0.0);
  } // end zero

  public static final void print(final Complex A[][])
  {
    int ni = A.length;
    int nj = A[0].length;
    for(int i=0; i<ni; i++)
      for(int j=0; j<nj; j++)
        System.out.println("A["+i+"]["+j+"]="+A[i][j]);
  } // end print

  public static final void multiply(final Complex A[][],
                                          final Complex B[], Complex C[])
  {
    int ni=A.length;
    int nj=A[0].length;
    if(B.length!=nj || C.length!=ni)
    {
 System.out.println("Error in ComplexMatrix.multiply,"+
                           " incompatible sizes.");
    }
    for(int i=0; i<ni; i++){
      C[i] = new Complex(0.0,0.0);
      for(int j=0; j<nj; j++){
        C[i] = C[i].add(A[i][j].multiply(B[j]));
      }
    }
  } // end multiply

  public static final void add(final Complex X[], final Complex Y[],
                                     Complex Z[])
  {
    int n=X.length;
    if(Y.length!=n || Z.length!=n)
    {
 System.out.println("Error in ComplexMatrix.add,"+
                           " incompatible sizes.");
    }
    for(int i=0; i<n; i++) Z[i] = X[i].add(Y[i]);
  } // end add

  public static final void subtract(final Complex X[], final Complex Y[],
                                          Complex Z[])
  {
    int n=X.length;
    if(Y.length!=n || Z.length!=n)
    {
 System.out.println("Error in ComplexMatrix.subtract,"+
                           " incompatible sizes.");
    }
    for(int i=0; i<n; i++) Z[i] = X[i].subtract(Y[i]);
  } // end subtract

  public static final double norm1(final Complex X[])
  {
    double norm=0.0;
    int n=X.length;
    for(int i=0; i<n; i++)
        norm = norm + X[i].abs();
    return norm;
  } // end norm1

  public static final double norm2(final Complex X[])
  {
    double norm=0.0;
    int n=X.length;
    for(int i=0; i<n; i++)
        norm = norm + X[i].abs()*X[i].abs();
    return StrictMath.sqrt(norm);
  } // end norm2

  public static final double normInf(final Complex X[])
  {
    double norm=0.0;
    int n=X.length;
    for(int i=0; i<n; i++)
      norm = Math.max(norm, X[i].abs());
    return norm;
  } // end normInf

  public static final void copy(final Complex X[], Complex Y[])
  {
    int n = X.length;
    if(Y.length!=n)
    {
      System.out.println("Error in ComplexMatrix.copy,"+
                           " incompatible sizes");
    }
    for(int i=0; i<n; i++) Y[i] = X[i];
  } // end copy

  public static final boolean equals(final Complex X[], final Complex Y[])
  {
    int n = X.length;
    boolean same = true;
    if(Y.length!=n)
    {
 System.out.println("Error in ComplexMatrix.equals,"+
                           " incompatible sizes");
    }
    for(int i=0; i<n; i++) same = same && X[i].equals(Y[i]);
    return same;
  } // end equals

  public static final void fromDouble(final double X[], Complex Z[])
  {
    int n = X.length;
    if(Z.length!=n)
    {
      System.out.println("Error in ComplexMatrix.fromDouble,"+
                         " incompatible sizes");
    }
    for(int i=0; i<n; i++) Z[i] = new Complex(X[i]);
  } // end fromDouble

  public static final void fromDouble(final double X[],
                                      final double Y[], Complex Z[])
  {
    int n = X.length;
    if(Z.length!=n || Y.length!=n)
    {
      System.out.println("Error in ComplexMatrix.fromDouble,"+
                         " incompatible sizes");
    }
    for(int i=0; i<n; i++) Z[i] = new Complex(X[i], Y[i]);
  } // end fromDouble

  public static void fromRoots(Complex X[], Complex Y[])
  {
    int n = X.length;
    if(Y.length!=n+1)
    {
      System.out.println("Error in ComplexMatrix.fromRoots,"+
                         " incompatible sizes");
    }
    Y[0] = X[0].negate();
    Y[1] = new Complex(1.0);
    if(n==1) return;
    for(int i=1; i<n; i++)
    {
      Y[i+1] = new Complex(0.0);
      for(int j=0; j<=i; j++)
      {
        Y[i+1-j] = Y[i-j].subtract(Y[i+1-j].multiply(X[i]));
      }
      Y[0] = Y[0].multiply(X[i]).negate();
    }
  }

  public static final void unitVector(Complex X[], int j)
  {
    int n = X.length;
    for(int i=0; i<n; i++)
      X[i]=new Complex(0.0);
    X[j]=new Complex(1.0);
  } // end unitVector

  public static final void zero(Complex X[])
  {
    int n = X.length;
    for(int i=0; i<n; i++)
      X[i]=new Complex(0.0);
  } // end zero

  public static final void print(final Complex X[])
  {
    int n = X.length;
    for(int i=0; i<n; i++)
      System.out.println("X["+i+"]="+X[i]);
  } // end print

  private static String input_file_name; // for read
  private static BufferedReader in; // for read
  private static String output_file_name; // for write
  private static BufferedWriter out; // for write
  private static PrintWriter file_out; // for write

  public static final void readSize(String file_name,  int rowCol[])
  {
    String input_line = new String("@"); // unread
    int len;       // input_line length
    int index;     // start of field
    int last;      // end of field
    String intStr; // for parseInt (does not ignore leading and 
                   // trailing white space)
    int ni;        // number of rows
    int nj;        // number of columns

    if(input_file_name==null || !file_name.equals(input_file_name))
    {
      input_file_name=file_name;
      try
      {
        in = new BufferedReader(new FileReader(file_name));
      }
      catch(Exception e)
      {
        System.out.println("ComplexMatrix.read unable to open file "+file_name);
        return;
      }
    }
    ni=0;
    nj=0;
    try
    {
      input_line = in.readLine(); // first read before  'while'
      while(input_line != null) // allow leading blank lines
      {
        input_line = input_line.trim();
        len = input_line.length();
        if(len==0)
        {
          input_line = in.readLine();
          continue; // skip blank lines
        }
        if(input_line.charAt(0)=='(')
        {
          System.out.println("ComplexMatrix.readSize unable to get size "+file_name);
          break;
        }
        index = 0;
        last = input_line.indexOf(' '); // first blank after number
        if(last==-1) last=len;
        intStr= input_line.substring(index,last);
        ni = Integer.parseInt(intStr);
        input_line = input_line.substring(last,len);
        input_line = input_line.trim();
        len = input_line.length();
        if(len==0)
        {
          nj=ni;
          break;
        }
        index = 0;
        last = input_line.indexOf(' '); // first blank after number
        if(last==-1) last=len;
        intStr= input_line.substring(index,last);
        nj = Integer.parseInt(intStr);
        break;
      } // end while
    }
    catch(Exception e)
    {
      System.out.println("ComplexMatrix.readSize unable to get size "+file_name);
    }
    rowCol[0]=ni;
    rowCol[1]=nj;
  } // end readSize

  public static final void read(String file_name, Complex A[][])
  {
    String input_line = new String("@"); // unread
    int len;       // input_line length
    int index;     // start of field
    int last;      // end of field
    String intStr; // for parseInt (does not ignore leading and 
                   // trailing white space)
    int i, ni;     // 0..ni-1
    int j, nj;     // 0..nj-1
    boolean have_line = false;

    if(input_file_name==null || !file_name.equals(input_file_name))
    {
      input_file_name=file_name;
      try
      {
        in = new BufferedReader(new FileReader(file_name));
      }
      catch(Exception e)
      {
        System.out.println("ComplexMatrix.read unable to open file "+file_name);
        return;
      }
    }
    ni=0;
    nj=0;
    try
    {
      input_line = in.readLine(); // first read before  'while'
      while(input_line != null) // allow leading blank lines
      {
        input_line = input_line.trim();
        len = input_line.length();
        if(len==0)
        {
          input_line = in.readLine();
          continue; // skip blank lines
        }
        if(input_line.charAt(0)=='(')
        {
          ni=A.length; // no size, just complex data
          nj=A[0].length;
          have_line=true;
          break;
        }
        index = 0;
        last = input_line.indexOf(' '); // first blank after number
        if(last==-1) last=len;
        intStr= input_line.substring(index,last);
        ni = Integer.parseInt(intStr);
        input_line = input_line.substring(last,len);
        input_line = input_line.trim();
        len = input_line.length();
        if(len==0)
        {
          nj=ni;
          break;
        }
        index = 0;
        last = input_line.indexOf(' '); // first blank after number
        if(last==-1) last=len;
        intStr= input_line.substring(index,last);
        nj = Integer.parseInt(intStr);
        break;
      } // end while
    }
    catch(Exception e)
    {
      System.out.println("ComplexMatrix.read unable to get size "+file_name);
    }
    // now read complex numbers
    i=0;
    j=0;
    if(A.length!=ni || A[0].length!=nj)
    {
      System.out.println("incompatible size in ComplexMatrix.read");
      return;
    }
    try
    {
      if(!have_line)
        input_line = in.readLine(); // first read before  'while'
      have_line=false;
      while(input_line != null)
      {
        input_line = input_line.trim();
        len = input_line.length();
        if(len==0)
        {
          input_line = in.readLine();
          continue; // skip blank lines
        }
        index = 0;
        last = input_line.indexOf(')'); // closing )
        if(last==-1)
        {
          input_line = in.readLine();
          continue;
        }
        intStr= input_line.substring(index,last+1);
        A[i][j] = Complex.parseComplex(intStr);
        j++;
        if(j==nj) { j=0; i++; }
        if(i==ni) break;
        input_line = input_line.substring(last+1);
      } // end while
    }
    catch(Exception e)
    {
      System.out.println("ComplexMatrix.read unable to read data "+file_name);
    }
  } // end read

  public static final void closeInput()
  {
    try
    {
      in.close();
    }
    catch(Exception e)
    {
      System.out.println("ComplexMatrix.closeInput not closed");
    }
    input_file_name = null;
  } // end closeInput


  public static final void write(String file_name, Complex A[][])
  {
    int ni = A.length;
    int nj = A[0].length;

    if(output_file_name==null || !file_name.equals(output_file_name))
    {
      output_file_name=file_name;
      try
      {
        out = new BufferedWriter(new FileWriter(file_name));
        file_out = new PrintWriter(out);
      }
      catch(Exception e)
      {
        System.out.println("ComplexMatrix.write unable to open file "+file_name);
        return;
      }
    }
    // write size
    if(ni == nj)
      file_out.println(ni);
    else
      file_out.println(ni+" "+nj);

    // now write complex numbers
    try
    {
      for(int i=0; i<ni; i++)
        for(int j=0; j<nj; j++)
          file_out.println(A[i][j].toString());
      file_out.println();
    }
    catch(Exception e)
    {
      System.out.println("ComplexMatrix.write unable to write data "+file_name);
    }
  } // end write

  public static final void closeOutput()
  {
    file_out.close();
    output_file_name = null;
  } // end closeOutput

} // end class ComplexMatrix