! ==================================================================== !
! -------------------------------------------------------------------- !
!                                                                      !
         SUBROUTINE GECP_S(A,NX,N,B,EPS,ISW,P,Q,VP,VS,VW,ICON)         !
!                                                                      !
! -------------------------------------------------------------------- !
!  <<DATE>>                                                            ! 
!     MAY 31ST, 2002                                                   !
!  <<CODER>>                                                           !
!     WATANABE, YOSHITAKA                                              !
!     COMPUTING AND COMMUNICATIONS CENTER, KYUSHU UNIERSITY            !
!  <<ABSTRACT>>                                           -1           !
!     SUBROUTINE `GECP_S' SOLVES THE REAL LINEAR SYSTEM X=A B          !
!     BY GAUSSIAN ELIMINATION WITH COMPLETE PIVOTING AND SCALING.      !
!  <<PRECISION>>                                                       !
!     SINGLE                                                           !
!  <<PARAMETER>>                                                       !
!      A: (INPUT)  REAL COEFFICIENT MARIX.                             ! 
!         (OUTPUT) LU DECOMPOSED MATRIX (SEE REMARK).                  !
!     NX: (INPUT)  AJUSTED SIZE OF MATRIX A.                           !
!      N: (INPUT)  DIMENSION OF THE LINEAR SYSTEM (NX>=N>1).           !
!      B: (INPUT)  RIGHT-HAND-SIDE VECTOR                              ! 
!         (OUTPUT) NUMERICAL APPLOXIMATE SOLUTION.                     !
!    EPS: (INPUT)  PARAMETER FOR CHECKING THE SINGULALITY (>=0).       !
!                  IF EPS<=0 THEN EPS IS SET BY MACHINE EPSILON.       !
!    ISW: (INPUT)  ISW=0 ...EXECUTE ONLY LU DECOMPOSITION.             !
!                  ISW=1 ...SOLVE LINEAR SYSTEM A*X=B USUALLY.         !
!                  ISW=2 ...SKIP LU DECOMPOSITION AND SOLVE A*X=B.     !
!                  ISW=3 ...SOLVE LINEAR SYSTEM (A^T)*X=B.             !  
!                  ISW=4 ...SKIP LU DECOMPOSITION AND SOLVE A(^T)*X=B. !
!                  ISW=OTHERS .... EXECUTE THE SAME AS ISW=1.          !
!      P: (OUTPUT) PIVOTING INFORMATION VECTOR FOR COLUMN.             !
!      Q: (OUTPUT) PIVOTING INFORNATION VECTOR FOR ROW.                !
!     VP: (OUTPUT) SCALING INFORMATION VECTOR FOR NORMALIZATION.       !
!     VS: (OUTPUT) SCALING INFORMATION VECTOR OF EQUATIONS.            !
!     VW: (OUT)    WORKING VECTOR.                                     !
!   ICON: (OUTPUT) CONDITION CODE.                                     !
!                      0 ... PROGRAM ENDED NORMALY.                    !
!                   2000 ... ALL ELEMENTS OF A COLUMN OR ROW ARE ZERO  !
!                            OR PIVOT LESS THAN EPS.                   !
!                   3000 ... PARAMETER ERROR.                          !
!  <<REMARK>>                                                          !
!    A IS DECOMPOSED LU WHERE                                          !
!      L ... THE LOWER TRIANGULAR MATRIX WHOSE ALL DIAGONAL ELEMENT    !
!            ARE 1.                                                    !
!      U ... THE UPPER TRIANGULAR MATRIX WHOSE DIAGONAL ELEMENTS ARE   !
!            STORED AS 1/U(I,I).                                       !
! -------------------------------------------------------------------- !
!                                                                      !
! ==================================================================== !
!            DECRALATION                                               !
! -------------------------------------------------------------------- !
 IMPLICIT NONE
!                                                                      !
 INTEGER,INTENT(IN)                               :: N,NX,ISW
 REAL(KIND(1.0E0)),DIMENSION(NX,NX),INTENT(INOUT) :: A
 REAL(KIND(1.0E0)),DIMENSION(NX),INTENT(INOUT)    :: B,VS,VP
 REAL(KIND(1.0E0)),DIMENSION(NX),INTENT(OUT)      :: VW
 REAL(KIND(1.0E0)),INTENT(IN)                     :: EPS
 INTEGER,DIMENSION(NX),INTENT(INOUT)              :: P,Q
 INTEGER,INTENT(OUT)                              :: ICON
!                                                                      !
 REAL(KIND(1.0E0))                                :: D,EPSZ
 INTEGER                                          :: I,J,K,L,C
!                                                                      !
 INTRINSIC EPSILON,MAX,ABS
!                                                                      !
! ==================================================================== !
!            CHECKING PARAMETERS                                       !
! -------------------------------------------------------------------- !
 ICON=0      
 IF(EPS<=0.0E0)  THEN 
   EPSZ=EPSILON(1.0E0)
 ELSE
   EPSZ=EPS
 END IF
 IF(N<1.OR.N>NX) THEN
   ICON=3000
   RETURN
 END IF
 IF(ICON/=0.AND.ICON/=2.AND.ICON/=3.AND.ICON/=4) ICON=1
!                                                                      !
! ==================================================================== !
!            NORMALIZATION                                             !
! -------------------------------------------------------------------- !
 IF(ISW/=2.AND.ISW/=4) THEN
!
! ---------- FOR EACH EQUATION --------------------------------------- !
   DO I=1,N
     D=0.0E0
     DO J=1,N
       D=MAX(ABS(A(I,J)),D)
     END DO
     IF(D==0.0E0) THEN
       ICON=2000
       RETURN
     END IF
     DO J=1,N
       A(I,J)=A(I,J)/D
     END DO
     VP(I)=D
   END DO
! ---------- FOR EACH UNKNOWN VALUE----------------------------------- !
   DO J=1,N
     D=0.0E0
     DO I=1,N
       D=MAX(ABS(A(I,J)),D)
     END DO
     IF(D==0.0E0) THEN
       ICON=2000
       RETURN
     END IF
     DO I=1,N
       A(I,J)=A(I,J)/D
     END DO
     VS(J)=D
   END DO
! ==================================================================== !
!            COMPLETE PIVOTING                                         !
! -------------------------------------------------------------------- !
   DO K=1,N
     P(K)=K               ! INITIAL VALUE
     Q(K)=K
   END DO
! ---------- SERCH PIVOT --------------------------------------------- !
   DO K=1,N
     D=0.0E0
     DO J=K,N
     DO I=K,N
       IF(ABS(A(I,J))>D) THEN
         D=ABS(A(I,J))
         L=I              ! COLUMN NUMBER
         C=J              ! ROW NUMBER
       END IF
     END DO
     END DO
     IF(D==0.0E0) THEN
       ICON=2000
       RETURN
     END IF
! ---------- CHANGE INDEX -------------------------------------------- !
     IF(K/=L) THEN
        DO J=1,N
          D=A(K,J)
          A(K,J)=A(L,J)  
          A(L,J)=D
        END DO
        J=P(K)
        P(K)=P(L)
        P(L)=J
     END IF
!                                                                      !
     IF(K/=C) THEN
        DO I=1,N
          D=A(I,K)
          A(I,K)=A(I,C)    
          A(I,C)=D
        END DO
        J=Q(K)
        Q(K)=Q(C)
        Q(C)=J
     END IF
!                                                                      !
! ==================================================================== !
!            LU DECOMPOSITION                                          !
! -------------------------------------------------------------------- !
     IF(ABS(A(K,K))<EPSZ) THEN
       ICON=2000
       RETURN
     END IF
     A(K,K)=1.0E0/A(K,K)
     DO I=K+1,N
       A(I,K)=A(I,K)*A(K,K)
     END DO
     DO J=K+1,N
     DO I=K+1,N
       A(I,J)=A(I,J)-A(I,K)*A(K,J)
     END DO 
     END DO
   END DO
 END IF
!                                                                      !
! ==================================================================== !
!            SOLVE LINEAR SYSTEM L*X=B                                 !
! -------------------------------------------------------------------- !
 IF(ISW==1.OR.ISW==2) THEN 
!                                                                      !
! ---------- FORWARD SUBSTITUTION -------------------------------------!
   DO I=1,N
     VW(I)=B(P(I))/VP(P(I))
     DO J=1,I-1
       VW(I)=VW(I)-A(I,J)*VW(J)
     END DO 
   END DO 
!                                                                      !
! ---------- BACKWARD SUBSTITUTION ----------------------------------- ! 
   DO I=N,1,-1
     B(I)=VW(I)
     DO J=I+1,N
       B(I)=B(I)-A(I,J)*B(J)
     END DO
     B(I)=B(I)*A(I,I)
   END DO
!                                                                      !
! ---------- REPLACING THE SOLUTION ---------------------------------- !
   DO I=1,N
     VW(I)=B(I) 
   END DO 
!                                                                      !
   DO I=1,N
     B(Q(I))=VW(I)/VS(Q(I))
   END DO 
 END IF
!                                                                      !
! ==================================================================== !
!            SOLVE LINEAR SYSTEM (L^T)*X=B                             !
! -------------------------------------------------------------------- !
 IF(ISW==3.OR.ISW==4) THEN 
!                                                                      !
! ---------- FORWARD SUBSTITUTION -------------------------------------!
   DO I=1,N
     VW(I)=B(Q(I))/VS(Q(I))
     DO J=1,I-1
       VW(I)=VW(I)-A(J,I)*VW(J)
     END DO
     VW(I)=VW(I)*A(I,I)
   END DO
!                                                                      !
! ---------- BACKWARD SUBSTITUTION ----------------------------------- ! 
   DO I=N,1,-1
     B(I)=VW(I)
     DO J=I+1,N
       B(I)=B(I)-A(J,I)*B(J)
     END DO
   END DO
!                                                                      !
! ---------- REPLACING THE SOLUTION ---------------------------------- !
   DO I=1,N
     VW(I)=B(I)
   END DO
!                                                                      !
   DO I=1,N
     B(P(I))=VW(I)/VP(P(I))
   END DO
 END IF
!                                                                      !
END SUBROUTINE GECP_S