SUBROUTINE dple(rowk, n, b, c, ierr) ! Code from the Naval Surface Warfare Center Mathematics Library. ! Solves A.c = b where the n x n matrix is not stored but read in, ! one row at a time. ! The amount of storage required is slightly over (n/2)^2 for large n, ! or just over a quarter of that which would be required to store A. ! The user must provide the routine ROWK to read in row K of the matrix. ! The rows are read sequentially; each row is read once only. ! N.B. Arguments D & IP have been removed. ! Code converted using TO_F90 by Alan Miller ! Date: 2000-10-20 Time: 22:49:06 ! ****************************************************************** ! SOLUTION OF LINEAR EQUATIONS WITH REDUCED STORAGE ! ****************************************************************** ! The partial pivot Henderson-Wassyng algorithm is used. ! Reference: ! Wassyng, A. "Solving Ax = b: A method with reduced storage requirements", ! SIAM J. Numerical Analysis, vol.19 (1982), pp.197-204. IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: b(:) REAL (dp), INTENT(OUT) :: c(:) INTEGER, INTENT(OUT) :: ierr INTERFACE SUBROUTINE rowk(n, k, row) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) INTEGER, INTENT(IN) :: n, k REAL (dp), INTENT(OUT) :: row(:) END SUBROUTINE rowk END INTERFACE ! Local variables REAL (dp) :: bk, cj, ck, c1, dkj INTEGER :: i, iflag, ij, ijold, ik, ip(n-1), j, k, kjold, km1, kp1, & last, lastm1, lcol, lcolp1, m, MAX, mjold, nm1, np1 REAL (dp), PARAMETER :: zero = 0.0_dp REAL (dp), ALLOCATABLE :: d(:) ! SET THE NECESSARY CONSTANTS ierr = 0 np1 = n + 1 MAX = n * n / 4 + n + 3 ALLOCATE( d(MAX) ) 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 ip(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 d(i-1) = -c(i) / c(1) END DO d(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 = ip(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 = MAX - n + k IF (k /= 2) THEN lastm1 = MAX - 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) + d(ij) * cj END DO bk = bk - d(ij+1) * cj END DO ! K=N CASE m = k IF (k >= n) THEN IF (c(k) == zero) GO TO 130 d(last) = bk / c(k) ELSE ! FIND THE PIVOT DO i = kp1, n IF (ABS(c(m)) < ABS(c(i))) m = i END DO ip(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 d(ik) = -c(i) / c(k) ik = ik + 1 END DO d(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 = d(mjold) d(mjold) = d(kjold) DO i = kp1, np1 ij = ij + 1 ijold = ijold + 1 ik = ik + 1 d(ij) = d(ijold) + d(ik) * dkj END DO END DO END DO last = MAX IF (iflag < 0) last = MAX - 2 d(n) = d(last) ! INSERT THE SOLUTION IN C c(1:n) = d(1:n) nm1 = n - 1 DO i = 1, nm1 k = n - i m = ip(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, row) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) INTEGER, INTENT(IN) :: n, k REAL (dp), INTENT(OUT) :: row(:) READ(UNIT=9) row(1:n) RETURN END SUBROUTINE rowk PROGRAM test_big_solve IMPLICIT NONE INTERFACE SUBROUTINE dple(rowk, n, b, c, ierr) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: b(:) REAL (dp), INTENT(OUT) :: c(:) INTEGER, INTENT(OUT) :: ierr INTERFACE SUBROUTINE rowk(n, k, row) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) INTEGER, INTENT(IN) :: n, k REAL (dp), INTENT(OUT) :: row(:) END SUBROUTINE rowk END INTERFACE END SUBROUTINE dple SUBROUTINE rowk(n, k, row) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) INTEGER, INTENT(IN) :: n, k REAL (dp), INTENT(OUT) :: row(:) END SUBROUTINE rowk END INTERFACE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) INTEGER :: i, ierr, n, row REAL (dp), ALLOCATABLE :: b(:), c(:) REAL (dp) :: errmax, total REAL :: finish, start WRITE(*, '(a)', ADVANCE='NO') ' Enter no. of rows (n): ' READ(*, *) n ALLOCATE( b(n), c(n) ) ! The solution will be: ! c = 1, 2, 3, ..., n ! Generate a matrix of random numbers which is written to a scratch file OPEN(UNIT=9, STATUS='SCRATCH', FORM='UNFORMATTED') DO row = 1, n CALL RANDOM_NUMBER( c ) total = 0 DO i = 1, n total = total + i * c(i) END DO b(row) = total WRITE(9) c END DO REWIND (9) CALL CPU_TIME(start) CALL dple(rowk, n, b, c, ierr) CALL CPU_TIME(finish) WRITE(*, '(a, F9.2, a)') ' Time taken = ', finish-start, 'sec.' errmax = 0.0_dp DO i = 1, n errmax = MAX( errmax, ABS(c(i)-i) ) END DO WRITE(*, '(a, g12.3)') ' Max. error = ', errmax STOP END PROGRAM test_big_solve