// sparse.java    Sparse storage of a matrix
//                specifically designed for solving PDE
//                simultaneous equations.
//                zero based subscripting
//                column [nrow] is right hand side

import java.io.*;

public class sparse
{
  cell A[];

  public sparse(int nrow) // constructor
  {
    A = new cell[nrow];
    for(int i=0; i<nrow; i++)
    {
      A[i] = new cell(); // every row must have one value
                         // else singular, no unique solution
    }
  } // end constructor sparse

  public sparse(double M[][]) // constructor
  {
    int n = M.length;
    if(M[0].length!=n+1)
    {
      System.out.println("sparse constructor error n="+n+" rows, yet columns="
                         +(M[0].length)+" rather than "+(n+1));
    }
    A = new cell[n];
    for(int i=0; i<n; i++)
    {
      A[i] = new cell(); // every row must have one value
                         // else singular, no unique solution
    }
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<=n; j++)
      {
	if(M[i][j]!=0.0) A[i].setval(j,M[i][j]);
      }
    }
  } // end constructor sparse


  public class cell // basic element in linked list, i is row index, j is col
  {
    int j;
    double val;
    cell next;

    public cell()
    {
      j = -1;
      next = null;
      val = 0.0;
    } // end default constructor

    public cell(int aj, double aval)
    {
      j = aj;
      next = null;
      val = aval;
    } // end j,val constructor

    public cell(int aj, cell anext, double aval)
    {
      j = aj;
      next = anext;
      val = aval;
    } // end j,next,val constructor

    public int getj()
    {
      return j;
    } // end get j of this cell

    public void setj(int aj)
    {
      j = aj;
    } // end set j of this cell

    public cell getnext()
    {
      return next;
    } // end get next of this cell

    public void setnext(cell anext)
    {
      next = anext;
    } // end set next of this cell

    public double getval()
    {
      return val;
    } // end get val of this cell

    public void setval(double aval)
    {
      val = aval;
    } // end set val of this cell

    public void addval(double aval)
    {
      val += aval;
    } // end add a val to this cell val

    public double getval(int aj)
    {
      cell temp = this;
      while(true)
      {
	if(aj==temp.j) return temp.val;
	if(aj<temp.j) // we passed it
	{
	  // System.out.println("getval("+aj+") no such j"); // not error
	  return 0.0;
	}
        if(temp.next!=null) // aj>temp.j
	{
	  temp = temp.next;
	}
	else
	{
	  // System.out.println("getval("+aj+") no such j"); // not error
	  return 0.0;
	}
      }
    } // end get val from cell j in this row, zero if not in list

    public void setval(int aj, double aval) // creates j if not in list
    {
      if(j<0 || j==aj)
      {
	j = aj;
	val = aval;
	return;
      }
      if(aj<j) // must insert before
      {
	next = new cell(j, next, val);
        j = aj;
        val = aval;
	return;
      }
      if(next==null) // must insert at end
      {
	next = new cell(aj, aval);
	return;
      }
      cell p = this; // will be previous, need its next
      cell q;
      while(true)
      {
        if(p.next!=null)
	{
	  q = p.next;
	}
        else // add to end of list
	{
	  p.next = new cell(aj, aval);
          return;
	}
	if(aj==q.j)
	{
	  p.next.setval(aval); // not q!  q is just a temporary
          return;
	}
	if(aj<q.j) // we passed it, insert before
	{
	  p.next.next = new cell(p.next.j, p.next.next, p.next.val);
          p.next.j = aj;
          p.next.val = aval;
	  return ;
	}
        if(q.next!=null) // move along list
	{
	  p = q;
	  q = q.next; 
	}
	else
	{
	  p.next.next = new cell(aj, aval);
	  return;
	}
      }
    } // end set val to j cell

    public void addval(int aj, double aval)
    {
      if(j<0 || j==aj)
      {
	j = aj;
	val += aval;
	return;
      }
      if(aj<j) // must insert before
      {
	next = new cell(j, next, val);
        j = aj;
        val = aval;
	return;
      }
      if(next==null) // must insert at end
      {
	next = new cell(aj, aval);
	return;
      }
      cell p = this; // will be previous, need its next
      cell q;
      while(true)
      {
        if(p.next!=null)
	{
	  q = p.next;
	}
        else // add to end of list
	{
	  p.next = new cell(aj,aval);
          return;
	}
	if(aj==q.j)
	{
	  p.next.addval(aval); // not q!  q is just a temporary
          return;
	}
	if(aj<q.j) // we passed it, insert before
	{
	  p.next.next = new cell(p.next.j, p.next.next, p.next.val);
          p.next.j = aj;
          p.next.val = aval;
	  return ;
	}
        if(q.next!=null) // move along list
	{
	  p = q;
	  q = q.next; 
	}
	else
	{
	  p.next.next = new cell(aj, aval);
	  return;
	}
      }
    } // end add a val to j cell val
  } // end class cell, end basic element in linked list


  public double get(int i, int j) // slow, improve for special cases
  {
    if(i>A.length || i<0)
    {
      System.out.println("sparse.get("+i+","+j+")  bad index");
      return 0.0;
    }
    return A[i].getval(j);
  } // end get

  public void put(int i, int j, double val)
  {
    if(i>A.length || i<0)
    {
      System.out.println("sparse.put("+i+","+j+")  bad index");
      return;
    }
    A[i].setval(j,val);
  } // end put

  public double[] getRHS()
  {
    cell B;
    int thisj;
    int n = A.length;
    double Y[] = new double[n];
    Boolean more;

    for(int i=0; i<n; i++)
    {
      if(A[i].getj() < 0)
      {
        Y[i] = 0.0;
        continue;
      }
      B = A[i];
      more = true;
      while(more)
      {
        thisj = B.getj();
        if(n==thisj)
        {
          Y[i] = B.getval();
          more = false;
        }
	else if(B.getnext()!=null)
	{
          B = B.next;
        }
	else
	{
          Y[i] = 0.0;
          more = false;
	  System.out.println("getRHS not found ("+i+","+n+")");
        } // end if  n==
      } // end while   
    } // end i loop
    return Y;
  } // end getRHS

  public void setRHS(double Y[])
  {
    int n = A.length;

    if(n!=Y.length)
    {
      System.out.println("ERROR sparse.setRHS wrong length of Y");
      System.out.println("Y.length="+Y.length+" yet should be A.length="+n);
      return;
    }
    for(int i=0; i<n; i++)
    {
      A[i].setval(n,Y[i]);
    }
  } // end setRHS

  public void add(int i, int j, double val)
  {
    if(i>=A.length || i<0)
    {
      System.out.println("sparse.add("+i+","+j+")  bad index");
      return;
    }
    A[i].addval(j, val);
  } // end add

  public void divide(int kr, int k) // works on a row j=k+1 .. n
  {
    // B[row[k]][j] = B[row[k]][j] / B[row[k]][k];
    cell Akk, Akj; // j > k traverse list to end
    int thisk;
    double valk, valj;

    Akk = A[kr];
    while(true)
    {
      thisk = Akk.getj();
      if(thisk>k)
      {
	System.out.println("divide("+kr+","+k+") ERROR, no k in this row!");
        return;
      } 
      else if(thisk==k)
      {
	  break;
      }
      else
      {
	if(Akk.getnext()!=null)
	{
          Akk = Akk.next;
	}
	else
	{
	  System.out.println("divide("+kr+","+k+") ERROR, no k this big");
          return;
	}
      } // end else
    } // end first while, positioned at col index k
    valk = Akk.getval();
    Akj = Akk;
    while(Akj.getnext()!=null)
    {
      valj = Akj.next.getval();
      Akj.next.setval(valj/valk);
      Akj = Akj.next;
    } // end second while
  } // end divide

  public void reduce(int ir, int k, int kr)  // works on complete row  j=k+1..n
  {
    // B[row[i]][j] = B[row[i]][j] - B[row[i]][k] * B[row[k]][j];
    cell Aij, Aij_prev, Aik, Akj;
    int thisij, thisik, thiskj; // j > k  traverse all j
    double valij, valik, valkj; 

    // find Aik and valik  constant for  while loop on j
    Aik = A[ir];
    thisik = Aik.getj();
    valik = 0.0;
    valij = 0.0;
    valkj = 0.0;
    if(thisik==k)
    {
      valik = Aik.getval();
    }
    else
    {
      while(thisik<=k)
      {
        if(thisik==k)
        {
	  valik = Aik.getval();
	  break;
        }
        else
	{
	  if(Aik.getnext()!=null)
	  {
	    Aik = Aik.next;
            thisik = Aik.getj();
	  }
          else
	  {
	    break; // valik is zero
	  } // end if !=
	} // end if ==k
      } // end while
    } // end ==k
    if(valik == 0.0) return; // no change possible
    Aij = Aik; // go from here.next
    Aij_prev = Aik;
    Akj = A[kr];

    // find valij and valkj (both move in j, no change if either not exists)
    while(true) // walk j, must have a pair with neither zero
    {
      thisij = Aij.getj();
      thiskj = Akj.getj();
      valij = Aij.getval();
      valkj = Akj.getval();
      if(thiskj<=k) // kj must be bigger than k
      {
	if(Akj.getnext()!=null)
	{
	  Akj = Akj.next;
          thiskj = Akj.getj();
          valkj = Akj.getval();
          continue;
	}
        else
	{
          return;
	}
      }
      if(thisij<thiskj) // ij must be bigger that k, equal to thiskj
      {
	if(Aij.getnext()!=null)
	{
	  Aij_prev = Aij;
	  Aij = Aij.next;
	  continue;
	}
        else
	{
          return;
	}
      }
      if(thisij>thiskj) // must insert a node
      {
	  Aij_prev.next = new cell(thiskj,Aij_prev.next,0.0);
          Aij = Aij_prev.next;
          thisij = Aij.getj();
          valij = Aij.getval();
      }
      // must be equal
      Aij_prev.next.setval(valij-valik*valkj);
      if(Aij.getnext()!=null)
      {
	Aij_prev = Aij;
        Aij = Aij.next;
      }
      else
      {
        return;
      }
      if(Akj.getnext()!=null)
      {
        Akj = Akj.next;
      }
      else
      {
        return;
      }
    } // end while walking j
  } // end reduce

  public void clear()
  {
    for(int i=0; i<A.length; i++)
    {
      A[i].j = -1;
      A[i].val = 0.0;
      A[i].next = null;
    }
  } // end clear

  public void write_all()
  {
    int j, jprev;
    double val;
    cell B;

    System.out.println("write_all");
    for(int i=0; i<A.length; i++)
    {
      jprev = -3;
      B = A[i];
      j = B.getj();
      while(j>=0)
      {
	val = B.getval();
	if(j<=jprev)
	{
	  System.out.println("i="+i+",j="+j+", val="+val+" ERROR in matrix");
	}
        else
        {
	  System.out.println("i="+i+",j="+j+", val="+val);
	}
        if(B.next != null)
	{
          B = B.next;
          jprev = j;
          j = B.getj();
	}
        else j = -2;
      } // end  while
    } // end  i
  } // end write_all

  // simeq      solve simultaneous equations AX=Y 
  // solve real linear equations for X where Y = A * X
  // method: Gauss-Jordan elimination using maximum pivot
  //    Translated to java by : Jon Squire , 26 March 2003
  //    First written by Jon Squire December 1959 for IBM 650, translated to
  //    other languages  e.g. Fortran converted to Ada converted to C
  //    then converted to java, then to sparse
  public void simeq(double X[]) // X returned solution
  {
    int hold , I_pivot;             // pivot indicies
    double pivot;                   // pivot element value
    double abs_pivot;
    int n = A.length;
    int row[] = new int[n];         // row interchange indicies

    // 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 = get(row[k],k); // B[row[k]][k];
      abs_pivot = Math.abs(pivot);
      I_pivot = k;
      for(int i=k+1; i<n; i++)
      {
        pivot = get(row[i],k); // B[row[i]][k])
	if(Math.abs(pivot) > abs_pivot) // B[row[i]][k])
        {
          I_pivot = i;
          abs_pivot = Math.abs(pivot);
        }
      }
      // 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; j++)
        {
          X[j] = 0.0;
        }
        System.out.println("redundant row (singular) "+row[k]);
        return;
      }
      else // can divide by pivot
      {
        // reduce about pivot
        // for(int j=k+1; j<n+1; j++)
        // {
	// B[row[k]][j] = B[row[k]][j] / B[row[k]][k];
	  divide(row[k],k); // j loop in here  should use pivot
	// }
        //  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] - B[row[i]][k] * B[row[k]][j];
              reduce(row[i],k,row[k]); // j loop in here on both
	                               // row[i] and row[k]
	      // }
          }
        }
      }
      //  finished inner reduction
    }
    //  end main reduction loop
    //  build  X  for return, unscrambling rows
    for(int i=0; i<n; i++)
    {
	X[i] = A[row[i]].getval(n); // this could be direct if matrix shrunk
    }
  } // end simeq

  public double[] multiply(double X[])
  {
    cell B;
    int thisj;
    int n = A.length;
    double Y[] = new double[n];

    for(int i=0; i<n; i++)
    {
      Y[i] = 0.0; // redundant for Java, yet multi languages
      B = A[i];
      while(B.getnext()!=null)
      {
	thisj = B.getj();
	Y[i] += X[thisj]*B.getval();
        B = B.next;
      } // end while
    } // end i loop
    return Y;
  }

  public double[][] export() // import uses a constructor
  {
    int n = A.length;
    double M[][] = new double[n][n+1]; // simultaneous equations
    int j, jprev;
    double val;
    cell B;

    for(int i=0; i<n; i++)
    {
      for(int jj=0; jj<=n; jj++) M[i][jj] = 0.0; 
      jprev = -3;
      B = A[i];
      j = B.getj();
      while(j>=0)
      {
	val = B.getval();
	if(j<=jprev)
	{
	  System.out.println("i="+i+",j="+j+", val="+val+
                             " ERROR in sparse matrix");
          return M;
	}
        else
        {
	  M[i][j] = val;
	}
        if(B.next != null)
	{
          B = B.next;
          jprev = j;
          j = B.getj();
	}
        else j = -2;
      } // end  while
    } // end  i
    return M;
  } // end export

} // end sparse.java