(*CROUT FACTORIZATION FOR TRIDIAGONAL LINEAR SYSTEMS ALGORITHM 6.7
*
*   To solve the n by n linear system
*
* E1: A[1,1] X[1] + A[1,2] X[2] + ... + A[1,n] X[n] = A[1,n+1] 
* E2: A[2,1] X[1] + A[2,2] X[2] + ... + A[2,n] X[n] = A[2,n+1] 
*  :
*  .
* EN: A[n,1] X[1] + A[n,2] X[2] + ... + A[n,n] X[n] = A[n,n+1] 
*
*   Input:  the dimension n; the entries of A.
*
*  Output:  solution x(1), x(2), ..., x(n). 
*)
Print["\n"];
Print["All diagonal entries, all lower sub-diagonal entries,\n
 all upper sub-diagonal entries, inhomogeneous term\n"]; 
Print["\n"];
Print["Place as many entries as desired on each line, but \n"];
Print["separate entries with at least one blank\n"];
Print["\n"];
Print["\n"];
AA = InputString["This is Crout Method for tridiagonal linear systems\n
     The array will be input from a text file in the order:(see screen)\n
     Has the input file been created?\n
     Enter 'yes' or 'no'\n"];
OK=0;
If[AA == "yes" || AA == "y" || AA == "Y",
   NAME = InputString["Input the file name in the form - \n
      drive:\ name.ext      for example\n
         A:\\DATA.DTA\n"];
   INP = OpenRead[NAME];
   OK = 0;
   While[OK == 0,
      n = Input["Input the number of equations - an integer\n"];
      If[n > 0,
	 For[i = 1,i <= n,i++,
	    A[i-1] = Read[INP,Number];
   	 ];
	 For[i = 2,i <= n,i++,
	    B[i-1] = Read[INP,Number];
   	 ];
	 NN = n-1;
	 For[i = 1,i <= NN,i++,
	    c[i-1]=Read[INP,Number];
   	 ];
	 For[i = 1,i <= n,i++,
	    BB[i-1]=Read[INP,Number];
   	 ];         
	 OK = 1;
	 Close[INP],
	 Input["Number must be a positive integer\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ],
   Input["This program will end so the input file\n
   can be created.\n
   \n
   Press 1 [enter] to continue\n"];
];
If[OK == 1,
   (* Steps 1-3 - Set up and solve LZ=B *)
   (* Step 1 *)
   c[0] = c[0]/A[0];
   Z[0] = BB[0]/A[0];
   (* Step 2 *)
   For[i = 2, i <= NN, i++,
      A[i-1] = A[i-1]-B[i-1]*c[i-2];
      c[i-1] = c[i-1]/A[i-1];
      Z[i-1] = (BB[i-1]-B[i-1]*Z[i-2])/A[i-1];
   ];
   (* Step 3 *)
   A[n-1] = A[n-1]-B[n-1]*c[n-2];
   Z[n-1] = (BB[n-1]-B[n-1]*Z[n-2])/A[n-1];
   (* Steps 4 & 5 - Solve UX=Z *)
   (* Step 4 *)
   X[n-1] = Z[n-1];
   (* Step 5 *)
   For[II = 1, II <= NN, II++,
      i = NN-II+1;
      X[i-1] = Z[i-1]-c[i-1]*X[i];
   ];
];
(* Step 6 *)
If[OK == 1,
   FLAG = Input["Select output destination\n
      1. Screen\n
      2. Text file\n
      Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
      For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout"
   ];
   Write[OUP,"CROUT METHOD FOR TRIDIAGONAL LINEAR SYSTEMS\n"];
   Write[OUP,"\n"];
   Write[OUP,"The solution is\n"];
   For[i = 1, i <= n, i++,
      Write[OUP,SetPrecision[X[i-1],9]];
   ];
   Write[OUP,"\n"];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];