'*************************************************************************
'* This program solves a tridiagonal linear set of equations M * U = R   *
'* --------------------------------------------------------------------- *
'* SAMPLE RUN:                                                           *
'*                                                                       *
'* Input data file (tridiag.dat):                                        *
'*                                                                       *
'*  5                                                                    *
'*  0.0   1.0   -2.0                                                     *
'* 10.0   5.0   -7.5                                                     *
'*  2.22 -3.25   0.00456                                                 *
'*  2.0   4.0   -3.3                                                     *
'* -1.0   3.0    0.0                                                     *
'*                                                                       *
'* (The first line contains the size of system                           *
'* Then the three columns give the three main diagonals of matrix M, the *
'* first column begins with a zero and the third column ends with zero)  *
'*                                                                       *
'* Output file (tridiag.lst):                                            *
'*                                                                       *
'* --------------------------------------------------------------------  *
'*   Tridiagonal system of equations                                     *
'* --------------------------------------------------------------------  *
'*  upper subdiagonal:                                                   *
'* -2.000000 -7.500000  0.004560 -3.300000  0.000000                     *
'*  main diagonal:                                                       *
'*  1.000000  5.000000 -3.250000  4.000000  3.000000                     *
'*  lower subdiagonal:                                                   *
'*  0.000000  10.00000  2.220000  2.000000 -1.000000                     *
'*  right side member:                                                   *
'*  1.000000  1.000000  1.000000  1.000000  1.000000                     *
'* --------------------------------------------------------------------  *
'*   System solution:                                                    *
'*  -0.136497 -0.568249 -0.694162  1.202870  0.734290                    *
'* --------------------------------------------------------------------  *
'*                                                                       *
'* Reference: "Numerical Recipes By W.H. Press, B.P. Flannery, S.A. Teu- *
'*             kolsky, W.T. Vetterling, Cambridge University Press, 1987 *
'*             [BIBLI 08].                                               *
'*                                                                       *
'*                       Basic version from FORTRAN by J-P Moreau, Paris *
'*                                     (www.jpmoreau.fr)                 *
'*************************************************************************
defint i-n
defdbl a-h,o-z

dim a(25), b(25), c(25), r(25), u(25)
'irc : INTEGER

  cls
  print
  print " Input datafile name (without .dat): "; : input nom$

  open nom$+".dat" for input as #1
  open nom$+".lst" for output as #2

  input #1, n

  for i=1 to n
    input #1, a(i),b(i),c(i)
    r(i)=1#
  next i

  close #1

  print #2,
  print #2,"-------------------------------------------------------"
  print #2,"    Tridiagonal system of equations"
  print #2,"-------------------------------------------------------"
  print #2,"    upper subdiagonal:"
  for i=1 to n
    print #2,using " ##.######"; c(i);
  next i
  print #2,
  print #2,"    main diagonal:"
  for i=1 to n
    print #2,using " ##.######"; b(i);
  next i
  print #2,
  print #2,"    lower subdiagonal:"
  for i=1 to n
    print #2,using " ##.######"; a(i);
  next i
  print #2,
  print #2,"    right side member:"
  for i=1 to n
    print #2,using " ##.######"; r(i);
  next i
  print #2,
  print #2,"-------------------------------------------------------"

  'TRIDAG(a,b,c,r,u,n,irc)
  gosub 1000

  if irc<>0 then
    print #2, " ERROR IN TRIDAG : RETURN CODE="; irc
  else
    print #2, "    System solution:"
  for i=1 to n
    print #2,using " ##.######"; u(i);
  next i
  end if
  print #2,
  print #2,"-------------------------------------------------------"

  print
  print " Results in ";nom$;".lst."
  close #2

END 'of main program


1000 'SUBROUTINE TRIDAG
'********************************************************************
'* Solves for a vector U of length N the tridiagonal linear set     *
'* M U = R, where A, B and C are the three main diagonals of matrix *
'* M(N,N), the other terms are 0. R is the right side vector.       *
'********************************************************************
  dim GAM(25)

  IF B(1)=0 THEN
    irc=1
    return
  end if

  BET=B(1)
  U(1)=R(1)/BET
  for J=2 to N                  'Decomposition and forward substitution
    GAM(J)=C(J-1)/BET
    BET=B(J)-A(J)*GAM(J)
    IF BET=0 THEN                   'Algorithm fails
      irc=2
      return
    end if
    U(J)=(R(J)-A(J)*U(J-1))/BET
  next J

  for J=N-1 to 1 step -1            'Back substitution
    U(J)=U(J)-GAM(J+1)*U(J+1)
  next J

  irc=0

Return

'End of file tridiag.bas