DECLARE SUB rowk (n%, k%, r#())
DECLARE SUB dple (n%, b#(), c#(), ierr%)
'***********************************************************************
'* Solve AX = B using a partial pivoting algorithm and reduced storage *
'* ------------------------------------------------------------------- *
'* SAMPLE RUN:                                                         *
'*                                                                     *
'* System to solve:                                                    *
'*   2.0000  -1.0000   1.0000   7.0000 -12.5400    5.0000              *
'*   1.0000   5.0000  -2.0000  -8.0000 100.0000    1.0000              *
'*   3.0000  -2.0000   3.0000  45.0000  27.3333    3.0000              *
'*  11.0000   0.5500  -2.0000  -4.0000   1.0000    4.0000              *
'*  33.0000   2.0000  -3.0000   5.0000   7.3333  -10.0000              *
'*                                                                     *
'* Solution is:                                                        *
'*  2.111496                                                           *
'* -25.829027                                                          *
'*  8.174232                                                           *
'* -2.521467                                                           *
'*  1.242104                                                           *
'*                                                                     *
'* ------------------------------------------------------------------- *
'* Ref.: "Wassyng, A. - Solving Ax = b: A method with reduced storage  *
'*        requirements, SIAM J. Numerical Analysis, vol.19 (1982),     *
'*        pp. 197-204".                                                *
'*                                                                     *
'*                                Basic Release By J-P Moreau, Paris.  *
'*                                        (www.jpmoreau.fr)            *
'***********************************************************************
'PROGRAM Test_DPLE
DEFDBL A-H, O-Z
DEFINT I-N

DIM a(25), b(25), soln(25)
DIM row AS INTEGER

' Define the right-hand side (n=5)
b(1) = 5#: b(2) = 1#: b(3) = 3#: b(4) = 4#: b(5) = -10#

CLS
PRINT
PRINT " System to solve:"
FOR i = 1 TO 5
  rowk 5, i, soln()
  FOR j = 1 TO 5
    PRINT USING " ###.####"; soln(j);
  NEXT j
  PRINT USING "  ###.####"; b(i)
NEXT i
PRINT

dple 5, b(), soln(), ierr

IF ierr = 0 THEN
  PRINT " Solution is:"
  FOR i = 1 TO 5
    PRINT USING "####.######"; soln(i)
  NEXT i
ELSE
  PRINT " Error = "; ierr
END IF

END 'of main program

SUB dple (n, b(), c(), ierr)
 
'******************************************************************
'       SOLUTION OF LINEAR EQUATIONS WITH REDUCED STORAGE
'******************************************************************

' Uses the Henderson-Wassyng partial pivot algorithm.
' Wassyng, A. "Solving Ax = b: A method with reduced storage requirements",
' SIAM J. Numerical Analysis, vol.19 (1982), pp. 197-204.

' The user must provide a routine ROWK to return the requested row of the
' matrix A.

  zero = 0#
DIM wk(n * n / 4 + n + 3)
DIM iwk(25)

'     SET THE NECESSARY CONSTANTS

ierr = 0
maxwk = n * n / 4 + n + 3
np1 = n + 1
k = 1
iflag = -1

'     GET THE FIRST COLUMN OF THE TRANSPOSED SYSTEM

rowk n, 1, c()
bk = b(1)

IF n <= 1 THEN
  IF c(1) = zero THEN GOTO 130
  c(1) = bk / c(1)
  EXIT SUB
END IF

'     FIND THE PIVOT FOR COLUMN 1

m = 1
FOR i = 2 TO n
  IF ABS(c(m)) < ABS(c(i)) THEN m = i
NEXT i

iwk(1) = m
c1 = c(m)
c(m) = c(1)
c(1) = c1
IF c(1) <> zero THEN
  
'     FIND THE FIRST ELEMENTARY MATRIX AND STORE IT IN D
  
  FOR i = 2 TO n
    wk(i - 1) = -c(i) / c(1)
  NEXT i
  wk(n) = bk / c(1)
  
'     K LOOP - EACH K FOR A NEW COLUMN OF THE TRANSPOSED SYSTEM
  
  FOR k = 2 TO n
    kp1 = k + 1
    km1 = k - 1
    
'       GET COLUMN K
    
    rowk n, k, c()
    FOR j = 1 TO km1
      m = iwk(j)
      cj = c(j)
      c(j) = c(m)
      c(m) = cj
    NEXT j
    bk = b(k)
    
    iflag = -iflag
    lcol = np1 - k
    lcolp1 = lcol + 1
    lastm1 = 1
    last = maxwk - n + k
    IF k <> 2 THEN
      lastm1 = maxwk - n + km1
      IF iflag < 0 THEN last = last - n + k - 2
      IF iflag > 0 THEN lastm1 = lastm1 - n + k - 3
    END IF
    
'     J LOOP - EFFECT OF COLUMNS 1 TO K-1 OF L-INVERSE
    
    FOR j = 1 TO km1
      cj = c(j)
      ij = (j - 1) * lcolp1
      IF j = km1 THEN ij = lastm1 - 1
      
'     I LOOP - EFFECT OF L-INVERSE ON ROWS K TO N+1
      
      FOR i = k TO n
	ij = ij + 1
	c(i) = c(i) + wk(ij) * cj
      NEXT i
      bk = bk - wk(ij + 1) * cj
    NEXT j
    
'       K=N CASE
    
    m = k
    IF k >= n THEN
      IF c(k) = zero THEN GOTO 130
      wk(last) = bk / c(k)
    ELSE
      
'     FIND THE PIVOT
      
      FOR i = kp1 TO n
	IF ABS(c(m)) < ABS(c(i)) THEN m = i
      NEXT i
      
      iwk(k) = m
      ck = c(m)
      c(m) = c(k)
      c(k) = ck
      IF c(k) = zero THEN GOTO 130
      
'     FIND THE K-TH ELEMENTARY MATRIX
      
      ik = last
      FOR i = kp1 TO n
	wk(ik) = -c(i) / c(k)
	ik = ik + 1
      NEXT i
      wk(ik) = bk / c(k)
    END IF
    
'     FORM THE PRODUCT OF THE ELEMENTARY MATRICES
    
    FOR j = 1 TO km1
      kjold = j * lcolp1 + k - np1
      mjold = kjold + m - k
      ij = (j - 1) * lcol
      ijold = ij + j
      IF j = km1 THEN
	kjold = lastm1
	mjold = lastm1 + m - k
	ijold = lastm1
      END IF
      
      ik = last - 1
      dkj = wk(mjold)
      wk(mjold) = wk(kjold)
      FOR i = kp1 TO np1
	ij = ij + 1
	ijold = ijold + 1
	ik = ik + 1
	wk(ij) = wk(ijold) + wk(ik) * dkj
      NEXT i
    NEXT j
  NEXT k
  
  last = maxwk
  IF iflag < 0 THEN last = maxwk - 2
  wk(n) = wk(last)
  
'     INSERT THE SOLUTION IN C
  
  FOR i = 1 TO n
    c(i) = wk(i)
  NEXT i
  
  nm1 = n - 1
  FOR i = 1 TO nm1
    k = n - i
    m = iwk(k)
    ck = c(k)
    c(k) = c(m)
    c(m) = ck
  NEXT i
  EXIT SUB
END IF

' THE SYSTEM IS SINGULAR

130 ierr = k

END SUB 'dple

SUB rowk (n, k, r())

IF k = 1 THEN
  r(1) = 2#: r(2) = -1#: r(3) = 1#: r(4) = 7#: r(5) = -12.54
ELSEIF k = 2 THEN
  r(1) = 1#: r(2) = 5#: r(3) = -2#: r(4) = -8#: r(5) = 100#
ELSEIF k = 3 THEN
  r(1) = 3#: r(2) = -2#: r(3) = 3#: r(4) = 45#: r(5) = 27.3333
ELSEIF k = 4 THEN
  r(1) = 11#: r(2) = .55: r(3) = -2#: r(4) = -4#: r(5) = 1#
ELSEIF k = 5 THEN
  r(1) = 33#: r(2) = 2#: r(3) = -3#: r(4) = 5#: r(5) = 7.3333
END IF

END SUB  'rowk