# TITLE:  LU Decomposition
# AUTHOR: Roy F.A. Maclean 
# EMAIL:  rfamgm at gmail
# WEB:    http://www.spiderpixel.co.uk/caspro        
# DATE:   25Apr2006
# MAKE:   CASIO
# MODEL:  9850G

# NOTES:  To solve Ax=B by LU Decomposition:  
#         Before running the program, put the square matrix
#         to be factored, into Mat A, and a column matrix into Mat B.
#         The matrices L and U are displayed.
#         The system becomes LUx=B
#         Ux is displayed, and then x is displayed.

# Example Put the following into Mat A
# 1 2 1
# 2 1 1
# 1 2 2
# and next bit into Mat B:
# 6
# 7
# 8
# Run program LU1, to get:
# L
# 1  0  0
# 2  1  0
# 1 -1  1
# U
# 1  1  1
# 0 -1 -1
# 0  0 -1

# Run program LU2, to get:
# UX
# 6
# -5
# -3
# X
# 1
# 2
# 3
# 

@@ Program "LU1"

Dim Mat A
List Ans[1]->N
Identity N->Mat L
Mat L->Mat U
For 1->K To N
Mat A[1,K]->Mat U[1,K]:Next
Mat A[1,1]->A
For 2->J To N
Mat A[J,1]%A->Mat L[J,1]:Next
For 2->J To N
For 2->K To N
If K<J:Then 0->T:For 1->S To K-1
T+Mat L[J,S]Mat U[S,K]->T:Next
(Mat A[J,K]-T)%Mat U[K,K]->Mat L[J,K]
Else 0->T:For 1->S To J-1
T+Mat L[J,S]Mat U[S,K]->T:Next
Mat A[J,K]-T->Mat U[J,K]:IfEnd
Next
Next
"L"_
Mat L_
"U"_
Mat U

@@ Program "LU2"

Mat B->Mat Y
For 2->J To N
0->D:For 1->K To J-1
Mat Y[K,1]Mat L[J,K]+D->D:Next
Mat B[J,1]-D->Mat Y[J,1]:Next
"UX"_
Mat Y_
(1%Mat U[N,N])Mat Y->Mat X
For N-1->J To 1 Step (-)1
0->D:For J+1->K To N
Mat X[K,1]Mat U[J,K]+D->D
Next
(Mat Y[J,1]-D)%Mat U[J,J]->Mat X[J,1]
Next
"X"_
Mat X