/*
 ##  class FTDiag
 ##
 ##  This class holds the routines for a tridiagonal matrix.  The storage
 ##  is the same as banded matrices.
 ##
 */

package edu.rice.linpack.Matrix.NMatrix;

import edu.rice.linpack.util.*;
import edu.rice.linpack.LNumber.*;

public class NTDiag extends NMatrix {
  
  public NTDiag() {
    super();
  }
  public NTDiag(int i, int j) {
    cons(3,j);
  }
  public NTDiag(LNumber[][] f) {
    cols = f.length;
    rows = 3;
    int cm = cols - 1;
    Mat = new LNumber[3][cols];
    Mat[1][0] = f[0][0];
    Mat[0][0] = f[0][1]; 
    Mat[2][0] = Mat[0][0].Clone();
    Mat[2][0].setZero();
    for(int i=1;i<cm;i++) {
      Mat[2][i] = f[i][i-1];
      Mat[1][i] = f[i][i];
      Mat[0][i] = f[i][i+1];
    }
    Mat[1][cm] = f[cm][cm];
    Mat[2][cm] = f[cm][cm-1];
    Mat[0][cm] = Mat[2][cm].Clone();
    Mat[0][cm].setZero();
  }
  public NTDiag(LNumber[] a, LNumber[] b, LNumber[] c) {
    cols = a.length;
    rows = 3;
    Mat = new LNumber[3][cols];
    Mat[1][0] = b[0];
    Mat[0][0] = c[0];
    Mat[2][0] = Mat[0][0].Clone();
    Mat[2][0].setZero();
    for(int i=1;i<cols-1;i++) {
      Mat[2][i] = a[i];
      Mat[1][i] = b[i];
      Mat[0][i] = c[i];
    }
    Mat[2][cols-1] = a[cols-1];
    Mat[1][cols-1] = b[cols-1];
    Mat[0][cols-1] = Mat[2][cols-1].Clone();
    Mat[0][cols-1].setZero();
  }
  public NTDiag(NTDiag F) {
    cons(F);
  }

  public void solve(LNumber[] B) 
       throws SingularMatrixException
  {
    Mat[2][0] = Mat[1][0].Clone();
    
    int nm = cols-1;
    if(cols > 0) {
      Mat[1][0] = Mat[0][0].Clone();
      Mat[0][0].setZero();
      Mat[0][nm].setZero();  

      for(int k=0;k<nm;k++) {
	int kp = k+1;

	if((Mat[2][kp].abs()).greaterOrEqual(Mat[2][k].abs())) {
	  this.swap(3,1,0,kp,this,1,0,k);
	  NUtil.swapElems(B,kp,k);
	}
	if(Mat[2][k].equals(0)) 
	  throw new SingularMatrixException(k+1);
	else {
	  LNumber T = (Mat[2][kp].div(Mat[2][k])).negateTo();
	  Mat[2][kp] = Mat[1][kp].add(T.mult(Mat[1][k]));
	  Mat[1][kp] = Mat[0][kp].add(T.mult(Mat[0][k]));
	  Mat[0][kp].setZero();
	  B[kp].addTo(T.mult(B[k]));
	}
      }
    }
    if(Mat[2][nm].equals(0)) 
      throw new SingularMatrixException(cols);
    else {
      int nm2 = nm - 1;
      B[nm].divTo(Mat[2][nm]);
      if(cols > 0) {
	B[nm2] = (B[nm2].sub(Mat[1][nm2].mult(B[nm]))).div(Mat[2][nm2]);
	if(nm2 > 0) {
	  for(int k=nm2-1;k>=0;k--) {
	    B[k] = ((B[k].sub(Mat[1][k].mult(B[k+1])))
		      .sub(Mat[0][k].mult(B[k+2]))).div(Mat[2][k]);
	  }
	}
      }
    }
  }
}