' ******************************************************************
' Solution of a system of linear equations by the SVD method
' ******************************************************************
' The SVD algorithm is applied to a linear system A * X = B, such
' that A(n x m) is the system matrix (with n >= m), B(n) is the
' constant vector and X(m) is the solution vector.
'
' The matrix A is factored into a product U * S * V' where U is a
' (n x m) column-orthogonal matrix, S is a (m x m) diagonal matrix
' with elements >= 0 (the singular values) and V is a (m x m)
' orthogonal matrix.
'
' If n > m, the algorithm gives the least squares solution.
'
' Function MathErr can be called to check the error code:
'   * MatOk      = no error
'   * MatErrDim  = non-compatible dimensions (n < m)
'   * MatNonConv = non-convergence of the iterative procedure
' ******************************************************************

#INCLUDE "fbmath.bi"

SUB PrintMatrix(Title AS STRING, A() AS DOUBLE)

  DIM AS INTEGER I, J, Lb, Ub1, Ub2

  PRINT : PRINT Title : PRINT

  Lb  = LBOUND(A, 1)
  Ub1 = UBOUND(A, 1)
  Ub2 = UBOUND(A, 2)

  FOR I = Lb TO Ub1
    FOR J = Lb TO Ub2
      PRINT USING "#####.######"; A(I, J);
    NEXT J
    PRINT
  NEXT I
END SUB

SUB PrintVector(Title AS STRING, B() AS DOUBLE)

  DIM AS INTEGER I, Lb, Ub

  PRINT : PRINT Title : PRINT

  Lb = LBOUND(B)
  Ub = UBOUND(B)

  FOR I = Lb TO Ub
    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 S(1 TO N)                    ' Singular values
DIM AS DOUBLE V(1 TO N, 1 TO N)            ' V factor
DIM AS DOUBLE X(1 TO N)                    ' Solution vector

CLS

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

SV_Decomp A(), S(), V()

SELECT CASE MathErr()
  CASE MatOk
    SV_Solve A(), S(), V(), B(), X()
    PrintVector "Solution vector:", X()
  CASE MatErrDim
    PRINT "Non-compatible dimensions!"
  CASE MatNonConv
    PRINT "Non-convergence!"
END SELECT

SLEEP
