' ******************************************************************
' Solution of a system of linear equations with
' single constant vector by Gauss-Jordan method
' ******************************************************************
' The algorithm solves the linear system A * X = B
'
' On input:
'   * A contains the system matrix
'   * B contains the constant vector
'
' On output:
'   * A contains the inverse matrix
'   * B contains the solution vector
'   * The determinant of the system matrix is returned in Det
'
' Function MathErr can be called to check the error code:
'   * MatOk     = no error
'   * MatErrDim = non-compatible dimensions
'   * MatSing   = quasi-singular matrix (Det ~ 0)
' ******************************************************************

#INCLUDE "fbmath.bi"

SUB PrintMatrix(Title AS STRING, A() AS DOUBLE)
' ------------------------------------------------------------------
' Print matrix A
' ------------------------------------------------------------------

  DIM AS INTEGER I, J

  PRINT : PRINT Title : PRINT

  FOR I = 1 TO UBOUND(A, 1)
    FOR J = 1 TO UBOUND(A, 2)
      PRINT USING "#####.######"; A(I, J);
    NEXT J
    PRINT
  NEXT I
END SUB

SUB PrintVector(Title AS STRING, B() AS DOUBLE)
' ------------------------------------------------------------------
' Print vector B
' ------------------------------------------------------------------

  DIM AS INTEGER I

  PRINT : PRINT Title : PRINT

  FOR I = 1 TO UBOUND(B)
    PRINT USING "#####.######"; B(I)
  NEXT I
END SUB

CONST N = 4                                ' Size of matrix

DIM AS DOUBLE A(1 TO N, 1 TO N) = _        ' System matrix
{{ 2,  1,   5, -8}, _
 { 7,  6,   2,  2}, _
 {-1, -3, -10,  4}, _
 { 2,  2,   2,  1}}

DIM AS DOUBLE B(1 TO N) = {0, 17, -10, 7}  ' Constant vector

DIM AS DOUBLE Det                          ' Determinant

CLS

PrintMatrix "System matrix:", A()
PrintVector "Constant vector:", B()

LinEq A(), B(), Det

SELECT CASE MathErr()
  CASE MatOk
    PrintMatrix "Inverse matrix:", A()
    PrintVector "Solution vector:", B()
    PRINT : PRINT USING "Determinant = #####.######"; Det
  CASE MatSing
    PRINT "Quasi-singular matrix!"
  CASE MatErrDim
    PRINT "Non-compatible dimensions!"
END SELECT

SLEEP
