#### Solution of Linear Algebraic Equations

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).
```