lubksb.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE lubksb
!                   *****************
!
     &(a,n,np,indx,b)
!
!***********************************************************************
! BIEF   V6P1                                  21/08/2010
!***********************************************************************
!
!brief    SOLVES THE SET OF N LINEAR EQUATIONS A \ DELTA X = B
!+                HERE A IS INPUT, NOT AS THE MATRIX A BUT AS ITS LU
!+                FACTORISATION, GIVEN BY THE SUBROUTINE LUDCMP. INDX
!+                IS INPUT AS THE PERMUTATION VECTOR RETURNED BY LUDCMP.
!+                B(1:N) IS INPUT AS THE RIGHT-HAND SIDE VECTOR B, AND
!+                RETURNS WITH THE SOLUTION VECTOR X. A, N, NP, AND
!+                INDX ARE NOT MODIFIED BY THIS SUBROUTINE AND CAN BE LEFT
!+                IN PLACE FOR SUCCESSIVE CALLS WITH DIFFERENT RIGHT-HAND
!+                SIDES B. THIS ROUTINE TAKES INTO ACCOUNT THE POSSIBILITY
!+                THAT B WILL BEGIN WITH A LOT OF 0 ELEMENTS, SO IT IS
!+                EFFICIENT FOR USE IN MATRIX INVERSION.
!
!history  CHUN WANG
!+        28/07/2006
!+        V5P7
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| MATRIX OF THE SYSTEM
!| B              |<->| RIGHT-HAND SIDE, THEN SOLUTION
!| INDX           |-->| ADDRESS IN RIGHT-HAND SIDE
!| N              |-->| SIZE OF B
!| NP             |-->| RANK OF MATRIX A
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: N,NP
      INTEGER, INTENT(IN) :: INDX(n)
      DOUBLE PRECISION, INTENT(INOUT) :: B(n)
      DOUBLE PRECISION, INTENT(IN)    :: A(np,np)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,II,J,LL
      DOUBLE PRECISION XSOM
!
!-----------------------------------------------------------------------
!
      ii=0
      ! WHEN II HAS A POSITIVE VALUE, IT WILL BECOME THE INDEX
      ! OF THE FIRST NONVANISHING ELEMENT OF B.
      ! DOES THE FORWARD SUBSTITUTION, EQUATION (2.3.6). THE ONLY
      ! NEW WRINKLE IS TO UNSCRAMBLE THE PERMUTATION AS WE GO.
!
      DO i=1,n
        ll=indx(i)
        xsom=b(ll)
        b(ll)=b(i)
        IF(ii.NE.0) THEN
          DO j=ii,i-1
            xsom=xsom-a(i,j)*b(j)
          ENDDO
        ELSEIF(xsom.NE.0.d0) THEN
          ii=i
          ! A NONZERO ELEMENT WAS ENCOUNTERED, SO FROM NOW ON
          ! WILL HAVE TO DO THE SUMS IN THE ABOVE LOOP
        ENDIF
        b(i)=xsom
      ENDDO
      DO i=n,1,-1 ! DOES THE BACKSUBSTITUTION, EQUATION (2.3.7)
        xsom=b(i)
        DO j=i+1,n
          xsom=xsom-a(i,j)*b(j)
        ENDDO
        b(i)=xsom/a(i,i) ! STORES A COMPONENT OF THE SOLUTION VECTOR X
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END