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

package edu.rice.linpack.Matrix.FMatrix;
import edu.rice.linpack.util.*;
import edu.rice.linpack.Vector.*;

public class FTDiag extends FMatrix {

  public FTDiag() {
    super();
  }
  public FTDiag(int i, int j) {
    cons(3,j);
  }
  public FTDiag(float[][] f) {
    cols = f.length;
    rows = 3;
    int cm = cols - 1;
    Mat = new float[3][cols];
    Mat[1][0] = f[0][0];
    Mat[0][0] = f[0][1]; 
    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];
  }
  public FTDiag(float[] a, float[] b, float[] c) {
    cols = a.length;
    rows = 3;
    Mat = new float[3][cols];
    Mat[1][0] = b[0];
    Mat[0][0] = c[0];
    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];
  }
  public FTDiag(FTDiag F) {
    cons(F);
  }

  public void solve(Vector Be) 
       throws WrongDataTypeException, SingularMatrixException
  {
    float[] B = Be.getFloatArray();

    this.Mat[2][0] = this.Mat[1][0];

    int nm = cols-1;
    if(cols > 0) {
      this.Mat[1][0] = this.Mat[0][0];
      this.Mat[0][0] = 0;
      this.Mat[0][nm] = 0;  

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

	if(Math.abs(this.Mat[2][kp]) >= Math.abs(this.Mat[2][k])) {
	  this.swap(3,1,0,kp,this,1,0,k);
	  FUtil.swapElems(B,kp,k);
	}
	if(this.Mat[2][k] == 0) 
	  throw new SingularMatrixException(k+1);
	else {
	  float T = -this.Mat[2][kp]/this.Mat[2][k];
	  this.Mat[2][kp] = this.Mat[1][kp] + T*this.Mat[1][k];
	  this.Mat[1][kp] = this.Mat[0][kp] + T*this.Mat[0][k];
	  this.Mat[0][kp] = 0;
	  B[kp] += T*B[k];
	}
      }
    }
    if(this.Mat[2][nm] == 0) 
      throw new SingularMatrixException(cols);
    else {
      int nm2 = nm - 1;
      B[nm] = B[nm]/this.Mat[2][nm];
      if(cols > 0) {
	B[nm2] = (B[nm2] - this.Mat[1][nm2]*B[nm])/this.Mat[2][nm2];
	if(nm2 > 0) {
	  for(int k=nm2-1;k>=0;k--) {
	    B[k] = (B[k] - this.Mat[1][k]*B[k+1] - 
		    this.Mat[0][k]*B[k+2])/this.Mat[2][k];
	  }
	}
      }
    }
  }
}