(*:Mathematica Version: 2.0 *)

(* :Copyright: Copyright 1992-1999, Wolfram Research, Inc.*)

(*:Name: LinearAlgebra`Tridiagonal` *)

(*:Context: LinearAlgebra`Tridiagonal` *)

(*:Title: Tridiagonal Matrix Routines *)

(*:Author:
	Jerry B. Keiper (Wolfram Research, Inc.)
*)

(*:Keywords:
	tridiagonal matrix, diagonals
*)

(*:Requirements: Length[a] >= Length[r] - 1,
                 Length[b] >= Length[r],
                 Length[c] >= Length[r] - 1. *)

(*:Warnings: Infinity introduced if pivot becomes zero. *)

(*:Sources:

*)

(*:Summary:
TridiagonalSolve solves A . x == r for x, where A is an n x n
tridiagonal matrix and r is a vector of length n.
*)

(*:Discussion:
TridiagonalSolve solves A . x == r for x, where A is an n x n
tridiagonal matrix and r is a vector of length n.
The three main diagonals of A are given by

     b[[1]] c[[1]]
     a[[1]] b[[2]] c[[2]]
            a[[2]]    .     .
                      .     .     .

                                     c[[n-1]]
                            a[[n-1]] b[[n]]  ,

Regular Gaussian elimination is employed, but no pivoting is done.
*)

BeginPackage["LinearAlgebra`Tridiagonal`"]

TridiagonalSolve::usage =
"TridiagonalSolve[ a, b, c, r ] solves A . x == r for x,
where A is a tridiagonal matrix.  The three diagonals of A are given
by: \n

     a11 = b[[1]], a12 = c[[1]],\n
      a21 = a[[1]], a22 = b[[2]], a23 = c[[2]], ...\n

No pivoting is done, so the algorithm may fail even though a solution exists.
This cannot happen for a symmetric positive definite matrix."

Unprotect[ TridiagonalSolve ];

Begin["LinearAlgebra`Tridiagonal`Private`"]

TridiagonalSolve[a_List, b_List, c_List, r_List]:=
    Module[{len=Length[r], solution=Array[0, Length[r]], aux, 
			aux1=Array[0, Length[r]], a1 = Prepend[a, 0], iter}, 
	aux = 1/b[[1]];
	solution[[1]] = r[[1]] aux;
	Do[aux1[[iter]] = c[[iter-1]] aux;
	    aux = 1/(b[[iter]]-a1[[iter]] aux1[[iter]]);
	    solution[[iter]] = (r[[iter]]-a1[[iter]] solution[[iter-1]]) aux,
	    {iter, 2, len}]; 			
	Do[solution[[iter]] -= aux1[[iter+1]] solution[[iter+1]],
	    {iter, len-1, 1, -1}];
	solution
	] /; Length[a] + 1 == Length[b] == Length[c] + 1 == Length[r]


End[]  (* LinearAlgebra`Tridiagonal`Private` *)

Protect[ TridiagonalSolve ];

EndPackage[] (* LinearAlgebra`Tridiagonal` *)

(*:Limitations: No pivoting done. *)


(*:Examples:

TridiagonalSolve[ Table[-1,{i,9}],Table[2,{i,10}],
    Table[-1,{i,9}],{1,0,0,0,0,0,0,0,0,0}]

LinearSolve[Table[Switch[Abs[i-j], 0, 2, 1, -1, _, 0], {i, 10}, {j, 10}],
	{1,0,0,0,0,0,0,0,0,0}]

TridiagonalSolve[{2,3,4,5}, {0,1,2,3,4}, {9,8,7,6}, {9,9,9,9,9}]

LinearSolve[ { { 0,9,0,0,0},{2,1,8,0,0},{0,3,2,7,0},{0,0,4,3,6},{0,0,0,5,4}},
     {9,9,9,9,9}]

*)
