Programmi

Algoritmo di Henderson-Wassyng

Soluzione di un sistema lineare di equazioni attraverso l'Algoritmo di Henderson-Wassyng

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(:)

• » Invia tramite EMAIL
• » Versione per la STAMPA
• » Le vostre opinioni

I link correlati all'argomento

• Software from Alan J. Miller