\ Linear equations with tridiagonal matrices by LU method
\   reference: Press, et al., "Numerical Recipes" (Cambridge
               U. Press, 1986)

\ ---------------------------------------------------
\     (c) Copyright 1999  Julian V. Noble.          \
\       Permission is granted by the author to      \
\       use this software for any application pro-  \
\       vided this copyright notice is preserved.   \
\ ---------------------------------------------------

\ This is an ANS Forth program requiring the
\   FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependences:
\       Assumes independent floating point stack

FALSE [IF]
Algorithm:

    Ax = r, A is tridiagonal

    b(k), k=1,...,n are the diagonal elements
    a(k), k=2,...,n are the lower subdiagonal elements
    c(k), k=1,...,n-1 are the upper subdiagonal elements

    Write A = LU where L is lower-triangular and U upper triangular

    If A were a general matrix this would be possible also but the
    decomposition would require O(n^3) steps. For the tridiagonal
    case, however, the decomposition requires only O(n) steps.

    Can assume L and U are bi-diagonal. In the case of L the lower
    subdiagonal is just a(k). In the case of U we can choose the
    diagonal elements to be 1 (unity). Thus we need to determine
    the diagonal elements of L and the upper subdiagonal of U.

    Call them L(k) and U(k) respectively. Then

        L(1) = b(1)
        U(k) = c(k)/L(k), k=1,...,n-1
        L(k) = b(k) - a(k)*U(k-1), k=2,...,n

    Finally let Ux = y and solve first

        Ly = r

    via

        y(1) = r(1)/L(1)
        y(2) = [r(2) - a(2)*y(1)] / L(2) etc.

    then

        x(n) = y(n)
        x(n-1) = y(n-1) - x(n)*U(n-1)  etc.

Usage:

    Say  a{ b{ c{ n }triangularize  r{ x{ n }backsolve

[THEN]


MARKER -tridiag

BL PARSE undefined DUP PAD C! PAD CHAR+ SWAP CHARS MOVE PAD FIND NIP 0=
[IF]  : undefined    BL WORD  FIND  NIP  0=  ;   [THEN]

include arrays.f
include ftran111.f

20 VALUE Nmax

Nmax long 1 FLOATS 1array a{    \ input array in vector form
Nmax long 1 FLOATS 1array b{
Nmax long 1 FLOATS 1array c{

0 VALUE aa{     0 VALUE bb{     0 VALUE cc{     0 VALUE NN

Nmax long 1 FLOATS 1array r{    \ inhomogeneous term


Nmax long 1 FLOATS 1array L{    \ diagonal of lower-triangular matrix
Nmax long 1 FLOATS 1array U{    \ subdiagonal of upper-triangular matrix

Nmax long 1 FLOATS 1array x{    \ solution vector

: }triangularize    ( a{ b{ c{ n --)
    TO NN   TO cc{  TO bb{  TO aa{
    f" bb{ 0 }"   FDUP  F0=  ABORT" Reduce # of equations by 1"
    L{ 0 }  F!
    f" U{ 0 } = cc{ 0 } / L{ 0 }"
    NN 1-  0  DO    f" U{ I } = cc{ I } / L{ I }"
                    f" L{I_1+} = bb{I_1+} - aa{I_1+} * U{I}"
    LOOP
;

: }backsolve    ( r{ x{ n --)
    TO NN   TO aa{  TO bb{
    f" bb{0} = bb{0} / L{0}"
    NN 1 DO     f" bb{ I } = (bb{ I } - a{I}*bb{I_1-}) / L{I}"
    LOOP
    f" aa{ NN_1- } = bb{ NN_1- }"
    0 NN 2 -  DO    f" aa{I} = bb{I} - U{I}*aa{I_1+}"
    LOOP
;