(*
From: levy@math.Princeton.EDU
Date: Thu, 21 Dec 89 14:09:05 -0500
To: dbae@math.Princeton.EDU, gunn@math.Princeton.EDU, wpt@math.Princeton.EDU,
        WRI!wri-tech@uunet.UU.NET
*)
BeginPackage["Tridiagonal`"]

Tridiagonal::usage="\n
  Tridiagonal[A,b] gives the solution x of the equation A'x=b, where\n
  A' is the n-by-n tridiagonal matrix with coefficients given in the\n
  n-by-3 matrix A.  (The entries A[[1,1]] and A[[-1,-1]] are ignored.)"

CyclicTridiagonal::usage="\n
  CyclicTridiagonal[A,b] gives the solution x of the equation A'x=b, where\n
  A' is the n-by-n ``cyclic'' tridiagonal matrix with coefficients\n
  given in the n-by-3 matrix A, with indices taken modulo n (that is,\n
  A'[[1,-1]]=A[[1,1]] and A'[[-1,1]]=A[[-1,-1]])."

Begin["Tridiagonal`private`"]

Tridiagonal[A_?MatrixQ,b_?VectorQ] :=
MapAt[(Drop[#,2]/#[[2]])&,Transpose[Append[Transpose[A],b]],{1}] //
  MapAt[0&,#,{-1,-2}]& // Accumulate[ForwardSubstitute,#]& // Reverse // 
  MapAt[Last,#,{1}]& // Accumulate[BackSubstitute,#]& // Reverse

ForwardSubstitute[{gamma_,delta_},{a_,b_,c_,d_}] := {c,d-a delta}/(b-a gamma)
ForwardSubstitute[{alpha_,gamma_,delta_},{a_,b_,c_,d_}] :=
	 (m-=k alpha; n-=k delta; k*= -gamma; {-a alpha,c,d-a delta}/(b-a gamma))

BackSubstitute[x_,{gamma_,delta_}] := delta-gamma x
BackSubstitute[x_,{alpha_,gamma_,delta_}] := delta-gamma x-alpha ult

CyclicTridiagonal[A_?MatrixQ,b_?VectorQ] :=
Block[{k,l,m,n,penult,ult},
  MapAt[(Drop[#,{2}]/#[[2]])&,Transpose[Append[Transpose[A],b]],{1}] //
    ({l,m,k,n}=Last[#];Drop[#,-1])& //
    Accumulate[ForwardSubstitute,#]& // Reverse //
    ({penult,ult}=LinearSolve[{{1,#[[1,1]]+#[[1,2]]},{k+l,m}},{#[[1,3]],n}];
     MapAt[penult&,#,{1}])& //
    Accumulate[BackSubstitute,#]& // Reverse // Append[#,ult]&
]

End[]
EndPackage[]
