SUBROUTINE dple(rowk, n, b, c, ierr)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-06-16  Time: 12:26:32
 
!     ******************************************************************
!     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.

! N.B. Arguments D and IP have been removed.

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)

INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(IN)   :: b(n)
REAL (dp), INTENT(OUT)  :: c(n)
INTEGER, INTENT(OUT)    :: ierr

! EXTERNAL rowk
INTERFACE
  SUBROUTINE rowk(n, k, r)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
    INTEGER, INTENT(IN)     :: n, k
    REAL (dp), INTENT(OUT)  :: r(:)
  END SUBROUTINE rowk
END INTERFACE

! Local variables

REAL (dp)  :: bk, cj, ck, c1, dkj
REAL (dp), PARAMETER  :: zero = 0.0_dp
REAL (dp)  :: wk(n*n/4 + n + 3)
INTEGER    :: i, iflag, ij, ijold, ik, iwk(n), j, k, kjold, km1, kp1,   &
              last, lastm1, lcol, lcolp1, m, maxwk, mjold, nm1, np1

!     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

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

IF (n <= 1) THEN
  IF (c(1) == zero) GO TO 130
  c(1) = bk / c(1)
  RETURN
END IF

!     FIND THE PIVOT FOR COLUMN 1

m = 1
DO  i = 2, n
  IF (ABS(c(m)) < ABS(c(i))) m = i
END DO

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
  
  DO  i = 2, n
    wk(i-1) = -c(i) / c(1)
  END DO
  wk(n) = bk / c(1)
  
!     K LOOP - EACH K FOR A NEW COLUMN OF THE TRANSPOSED SYSTEM
  
  DO  k = 2, n
    kp1 = k + 1
    km1 = k - 1
    
!       GET COLUMN K
    
    CALL rowk(n, k, c)
    DO  j = 1, km1
      m = iwk(j)
      cj = c(j)
      c(j) = c(m)
      c(m) = cj
    END DO
    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) last = last - n + k - 2
      IF (iflag > 0) lastm1 = lastm1 - n + k - 3
    END IF
    
!     J LOOP - EFFECT OF COLUMNS 1 TO K-1 OF L-INVERSE
    
    DO  j = 1, km1
      cj = c(j)
      ij = (j-1) * lcolp1
      IF (j == km1) ij = lastm1 - 1
      
!     I LOOP - EFFECT OF L-INVERSE ON ROWS K TO N+1
      
      DO  i = k, n
        ij = ij + 1
        c(i) = c(i) + wk(ij) * cj
      END DO
      bk = bk - wk(ij+1) * cj
    END DO
    
!       K=N CASE
    
    m = k
    IF (k >= n) THEN
      IF (c(k) == zero) GO TO 130
      wk(last) = bk / c(k)
    ELSE
      
!     FIND THE PIVOT
      
      DO  i = kp1, n
        IF (ABS(c(m)) < ABS(c(i))) m = i
      END DO
      
      iwk(k) = m
      ck = c(m)
      c(m) = c(k)
      c(k) = ck
      IF (c(k) == zero) GO TO 130
      
!     FIND THE K-TH ELEMENTARY MATRIX
      
      ik = last
      DO  i = kp1, n
        wk(ik) = -c(i) / c(k)
        ik = ik + 1
      END DO
      wk(ik) = bk / c(k)
    END IF
    
!     FORM THE PRODUCT OF THE ELEMENTARY MATRICES
    
    DO  j = 1, 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)
      DO  i = kp1, np1
        ij = ij + 1
        ijold = ijold + 1
        ik = ik + 1
        wk(ij) = wk(ijold) + wk(ik) * dkj
      END DO
    END DO
  END DO
  
  last = maxwk
  IF (iflag < 0) last = maxwk - 2
  wk(n) = wk(last)
  
!     INSERT THE SOLUTION IN C
  
  c(1:n) = wk(1:n)
  
  nm1 = n - 1
  DO  i = 1, nm1
    k = n - i
    m = iwk(k)
    ck = c(k)
    c(k) = c(m)
    c(m) = ck
  END DO
  RETURN
END IF

!     THE SYSTEM IS SINGULAR

130 ierr = k
RETURN
END SUBROUTINE dple



SUBROUTINE rowk(n, k, r)

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
INTEGER, INTENT(IN)     :: n, k
REAL (dp), INTENT(OUT)  :: r(:)

INTEGER    :: i
REAL (dp)  :: temp

r(k) = 1.0_dp
IF (k > 1) THEN
  temp = 1.0_dp
  DO i = 1, k-1
    temp = 0.99_dp * temp
    r(k-i) = temp
  END DO
END IF

IF (k < n) THEN
  temp = 1.0_dp
  DO i = 1, n-k
    temp = 0.99_dp * temp
    r(k+i) = temp
  END DO
END IF

RETURN
END SUBROUTINE rowk



PROGRAM Test_DPLE
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)

REAL (dp)  :: a(100), x(100), b(100), error, temp, soln(100)
INTEGER    :: i, ierr, row

INTERFACE
  SUBROUTINE dple(rowk, n, b, c, ierr)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
    INTEGER, INTENT(IN)     :: n
    REAL (dp), INTENT(IN)   :: b(n)
    REAL (dp), INTENT(OUT)  :: c(n)
    INTEGER, INTENT(OUT)    :: ierr
    INTERFACE
      SUBROUTINE rowk(n, k, r)
        IMPLICIT NONE
        INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
        INTEGER, INTENT(IN)     :: n, k
        REAL (dp), INTENT(OUT)  :: r(:)
      END SUBROUTINE rowk
    END INTERFACE
  END SUBROUTINE dple

  SUBROUTINE rowk(n, k, r)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)
    INTEGER, INTENT(IN)     :: n, k
    REAL (dp), INTENT(OUT)  :: r(:)
  END SUBROUTINE rowk
END INTERFACE

! Generate the solution using the random number generator
CALL RANDOM_NUMBER(x)
! Calculate the right-hand side
DO row = 1, 100
  CALL rowk(100, row, a)
  b(row) = DOT_PRODUCT(a, x)
END DO

CALL dple(rowk, 100, b, soln, ierr)
IF (ierr == 0) THEN
  temp = 0.0
  DO i = 1, 100
    error = ABS(soln(i) - x(i))
    temp = MAX(temp, error)
  END DO
  WRITE(*, '(a, g13.5)') ' Max. error = ', temp
ELSE
  WRITE(*, *) ' Error = ', ierr
END IF

STOP
END PROGRAM Test_DPLE