This shows that numerical analysis algorithms can be implemented elegantly in Prolog.
check_determ
domains
scalar = real
vector = scalar*
matrix = vector*
predicates
gauss_jordan_elimination(matrix,matrix)
gauss_jordan_elimination_1(matrix,matrix,matrix)
gaussian_elimination_with_backsubstitution(matrix,vector)
gaussian_elimination(matrix,matrix,matrix)
backsubstitution(matrix,vector,vector)
backsubstitution_1(vector,vector,scalar,scalar)
pivot_row(matrix,vector,matrix)
normalized_pivot_row(vector,vector)
normalized_pivot_row_1(vector,scalar,vector)
normalized_other_rows(matrix,vector,matrix)
normalized_other_row(vector,scalar,vector,vector)
reverse(matrix,matrix,matrix)
clauses
/*-------------------------------------------------------------------------*/
/* Gauss-Jordan Elimination with Partial Pivoting */
/* It handles any number of right-hand sides */
/* It can be used to invert a matrix */
/* */
/* e.g. gauss_jordan_elimination([[2,1,3, 13, 7], Gives [[1,1], */
/* [3,2,1, 10, 8], [2,2], */
/* [2,4,1, 13,11]], X). [3,1]] */
/*-------------------------------------------------------------------------*/
gauss_jordan_elimination(Ass, Xss):-
gauss_jordan_elimination_1(Ass, [], Xss).
gauss_jordan_elimination_1([], ReversedXss, Xss):-
!,
reverse(ReversedXss, [], Xss).
gauss_jordan_elimination_1(Lower, Upper, Xss):-
pivot_row(Lower, PivotRow, OtherRows),
normalized_pivot_row(PivotRow, NormalizedPivotRow),
normalized_other_rows(OtherRows, NormalizedPivotRow, NewLower),
normalized_other_rows(Upper, NormalizedPivotRow, NewUpper),
gauss_jordan_elimination_1(NewLower, [NormalizedPivotRow|NewUpper], Xss).
pivot_row([PivotRow], PivotRow, []).
pivot_row([[X1|Row1],[X2|Row2]|Rows], PivotRow, [[X2|Row2]|Rest]):-
abs(X1) > abs(X2),
!,
pivot_row([[X1|Row1]|Rows], PivotRow, Rest).
pivot_row([Row1,Row2|Rows], PivotRow, [Row1|Rest]):-
pivot_row([Row2|Rows], PivotRow, Rest).
normalized_pivot_row([0|_], _):-
!,
write("No non-zero pivot\n"),
fail.
normalized_pivot_row([Pivot|As], Bs):-
normalized_pivot_row_1(As, Pivot, Bs).
normalized_pivot_row_1([], _, []).
normalized_pivot_row_1([A|As], Pivot, [B|Bs]):-
B = A / Pivot,
normalized_pivot_row_1(As, Pivot, Bs).
normalized_other_rows([], _, []).
normalized_other_rows([[A|As]|Ass], PivotRow, [Bs|Bss]):-
normalized_other_row(As, A, PivotRow, Bs),
normalized_other_rows(Ass, PivotRow, Bss).
normalized_other_row([], _, [], []).
normalized_other_row([A|As], X, [P|Ps], [B|Bs]):-
B = A - X * P,
normalized_other_row(As, X, Ps, Bs).
/* reverse(Xs, Ys, Zs) is true if Zs is the result of concatenating the */
/* reverse of the list Xs to the front of the list Ys. */
reverse([], Ys, Ys).
reverse([X|Xs], Ys, Zs):-reverse(Xs, [X|Ys], Zs).
/*-------------------------------------------------------------------------*/
/* Gaussian Elimination with Partial Pivoting and Backsubstitution */
/* It handles only one right-hand side */
/* */
/* e.g. gaussian_elimination_with_backsubstitution([[2,1,3, 13], gives [1, */
/* [3,2,1, 10], 2, */
/* [2,4,1, 13]], X). 3] */
/*-------------------------------------------------------------------------*/
gaussian_elimination_with_backsubstitution(Ass, Xs):-
gaussian_elimination(Ass, [], ReversedReducedAss),
backsubstitution(ReversedReducedAss, [], Xs).
gaussian_elimination([], Xss, Xss):-!.
gaussian_elimination(Lower, Upper, Xss):-
pivot_row(Lower, PivotRow, OtherRows),
normalized_pivot_row(PivotRow, NormalizedPivotRow),
normalized_other_rows(OtherRows, NormalizedPivotRow, NewLower),
gaussian_elimination(NewLower, [NormalizedPivotRow|Upper], Xss).
backsubstitution([], Xs, Xs).
backsubstitution([As|Ass], Ys, Xs):-
backsubstitution_1(Ys, As, 0, Y),
backsubstitution(Ass, [Y|Ys], Xs).
backsubstitution_1([], [B], Acc0, Acc):-
Acc = B - Acc0.
backsubstitution_1([X|Xs], [A|As], Acc0, Acc):-
Acc1 = Acc0 + X * A,
backsubstitution_1(Xs, As, Acc1, Acc).