module GECP
  interface linear_solve
        function linear_solve_s(A,b,eps,sv,delta) result(x)
          real(kind(1E0)),dimension(:,:),intent(in) :: A
          real(kind(1E0)),dimension(:),intent(in)   :: b
          real(kind(1E0)),intent(in),optional       :: eps
          real(kind(1E0)),intent(out),optional      :: sv
          real(kind(1E0)),intent(in),optional       :: delta
          real(kind(1E0)),dimension(size(b))        :: x
        end function linear_solve_s

        function linear_solve_d(A,b,eps,sv,delta) result(x)
          real(kind(1D0)),dimension(:,:),intent(in) :: A
          real(kind(1D0)),dimension(:),intent(in)   :: b
          real(kind(1D0)),intent(in),optional       :: eps
          real(kind(1D0)),intent(out),optional      :: sv
          real(kind(1D0)),intent(in),optional       :: delta
          real(kind(1D0)),dimension(size(b))        :: x
        end function linear_solve_d

        function linear_solve_q(A,b,eps,sv,delta) result(x)
          real(kind(1Q0)),dimension(:,:),intent(in) :: A
          real(kind(1Q0)),dimension(:),intent(in)   :: b
          real(kind(1Q0)),intent(in),optional       :: eps
          real(kind(1Q0)),intent(out),optional      :: sv
          real(kind(1Q0)),intent(in),optional       :: delta
          real(kind(1Q0)),dimension(size(b))        :: x
        end function linear_solve_q
  end interface
end module GECP
! =========================================================================== !
!                                                                             !
!                SOLVE THE LINEAR EQUATIONS                                   !
!                                                                             !
! =========================================================================== !

! --------------------- SINGLE PRECISION REAL ------------------------------- !

 function linear_solve_s(A,b,eps,sv,delta) result(x)
   implicit none
   real(kind(1E0)),dimension(:,:),intent(in) :: A
   real(kind(1E0)),dimension(:),intent(in)   :: b
   real(kind(1E0)),intent(in),optional       :: eps
   real(kind(1E0)),intent(out),optional      :: sv
   real(kind(1E0)),intent(in),optional       :: delta
   real(kind(1E0)),dimension(size(b))        :: x

   integer,parameter                         :: IOUT=6
   integer                                   :: n,isw=1,icon,k1,k2,KMAX=100,k
   real(kind(1E0))                           :: epsz,deltaz,s
   real(kind(1E0)),dimension(:,:),allocatable:: Z
   real(kind(1E0)),dimension(:),allocatable  :: vp,vs,vw,y,u
   integer,dimension(:),allocatable          :: ip,iq

   intrinsic size,present,sum,sqrt,abs,maxval
   external GECP_S

   n=size(b); k1=size(A,1); k2=size(A,2); epsz=0.0E0; deltaz=1.0E-5
   if(present(eps)) epsz=eps
   if(present(delta)) deltaz=delta

   if( k1<n .or. k2<n .or. n<1 .or. epsz<0.0E0 .or. deltaz<0.0E0 ) then
     write(IOUT,*) '*** PARAMETER ERROR IN LINEAR_SOLVE'
     return
   end if

   allocate(Z(k1,k2),vp(n),vs(n),vw(n),ip(n),iq(n))
   x=b
   Z=A

   call GECP_S(Z,k1,n,x,epsz,isw,ip,iq,vp,vs,vw,icon)
     if( icon /= 0 ) then
       write(IOUT,*) '*** LINEAR_SOLVE FAILED;'
       write(IOUT,*) '*** IT IS HIGHLY PROBABLE THAT THE MATRIX IS SINGULAR.'
       deallocate(Z,vp,vs,vw,ip,iq)
       return
     end if

! -------------- inverse power method ------------------- !
   if(present(sv)) then 
     allocate(y(n),u(n))
     y=1; y=y/sqrt(sum(y**2))
     if(present(delta)) deltaz=delta
     do k=1,KMAX
       u=y
       isw=2; call GECP_S(Z,k1,n,y,epsz,isw,ip,iq,vp,vs,vw,icon)
       isw=4; call GECP_S(Z,k1,n,y,epsz,isw,ip,iq,vp,vs,vw,icon)
       s=sqrt(sum(y**2))
       y=y/s
       if(maxval(abs(y-u))/maxval(abs(u))<deltaz) exit
       if(k==KMAX) then
        write(IOUT,*) '*** INVERSE ITERATION DID NOT CONVERGE UNTIL',k,'-TIMES'
        write(IOUT,*) '*** IN LINEAR_SOLVE'
       end if
     end do
     sv=1.0E0/sqrt(s)
     deallocate(y,u)
   end if

   deallocate(Z,vp,vs,vw,ip,iq)
 end function linear_solve_s


! --------------------- DOUBLE PRECISION REAL ------------------------------- !

 function linear_solve_d(A,b,eps,sv,delta) result(x)
   implicit none
   real(kind(1D0)),dimension(:,:),intent(in) :: A
   real(kind(1D0)),dimension(:),intent(in)   :: b
   real(kind(1D0)),intent(in),optional       :: eps
   real(kind(1D0)),intent(out),optional      :: sv
   real(kind(1D0)),intent(in),optional       :: delta
   real(kind(1D0)),dimension(size(b))        :: x

   integer,parameter                         :: IOUT=6
   integer                                   :: n,isw=1,icon,k1,k2,KMAX=100,k
   real(kind(1D0))                           :: epsz,deltaz,s
   real(kind(1D0)),dimension(:,:),allocatable:: Z
   real(kind(1D0)),dimension(:),allocatable  :: vp,vs,vw,y,u
   integer,dimension(:),allocatable          :: ip,iq

   intrinsic size,present,sum,sqrt,abs,maxval
   external GECP_D

   n=size(b); k1=size(A,1); k2=size(A,2); epsz=0.0D0; deltaz=1.0D-10
   if(present(eps)) epsz=eps
   if(present(delta)) deltaz=delta

   if( k1<n .or. k2<n .or. n<1 .or. epsz<0.0D0 .or. deltaz<0.0D0 ) then
     write(IOUT,*) '*** PARAMETER ERROR IN LINEAR_SOLVE'
     return
   end if

   allocate(Z(k1,k2),vp(n),vs(n),vw(n),ip(n),iq(n))
   x=b
   Z=A

   call GECP_D(Z,k1,n,x,epsz,isw,ip,iq,vp,vs,vw,icon)
     if( icon /= 0 ) then
       write(IOUT,*) '*** LINEAR_SOLVE FAILED;'
       write(IOUT,*) '*** IT IS HIGHLY PROBABLE THAT THE MATRIX IS SINGULAR.'
       deallocate(Z,vp,vs,vw,ip,iq)
       return
     end if

! -------------- inverse power method ------------------- !
   if(present(sv)) then 
     allocate(y(n),u(n))
     y=1; y=y/sqrt(sum(y**2))
     if(present(delta)) deltaz=delta
     do k=1,KMAX
       u=y
       isw=2; call GECP_D(Z,k1,n,y,epsz,isw,ip,iq,vp,vs,vw,icon)
       isw=4; call GECP_D(Z,k1,n,y,epsz,isw,ip,iq,vp,vs,vw,icon)
       s=sqrt(sum(y**2))
       y=y/s
       if(maxval(abs(y-u))/maxval(abs(u))<deltaz) exit
       if(k==KMAX) then
        write(IOUT,*) '*** INVERSE ITERATION DID NOT CONVERGE UNTIL',k,'-TIMES'
        write(IOUT,*) '*** IN LINEAR_SOLVE'
       end if
     end do
     sv=1.0D0/sqrt(s)
     deallocate(y,u)
   end if

   deallocate(Z,vp,vs,vw,ip,iq)
 end function linear_solve_d


! ---------------------  QUADRUPLE PRECISION REAL --------------------------- !

 function linear_solve_q(A,b,eps,sv,delta) result(x)
   implicit none
   real(kind(1Q0)),dimension(:,:),intent(in) :: A
   real(kind(1Q0)),dimension(:),intent(in)   :: b
   real(kind(1Q0)),intent(in),optional       :: eps
   real(kind(1Q0)),intent(out),optional      :: sv
   real(kind(1Q0)),intent(in),optional       :: delta
   real(kind(1Q0)),dimension(size(b))        :: x

   integer,parameter                         :: IOUT=6
   integer                                   :: n,isw=1,icon,k1,k2,KMAX=100,k
   real(kind(1Q0))                           :: epsz,deltaz,s
   real(kind(1Q0)),dimension(:,:),allocatable:: Z
   real(kind(1Q0)),dimension(:),allocatable  :: vp,vs,vw,y,u
   integer,dimension(:),allocatable          :: ip,iq

   intrinsic size,present,sum,sqrt,abs,maxval
   external GECP_Q

   n=size(b); k1=size(A,1); k2=size(A,2); epsz=0.0Q0; deltaz=1.0Q-20
   if(present(eps)) epsz=eps
   if(present(delta)) deltaz=delta

   if( k1<n .or. k2<n .or. n<1 .or. epsz<0.0Q0 .or. deltaz<0.0Q0 ) then
     write(IOUT,*) '*** PARAMETER ERROR IN LINEAR_SOLVE'
     return
   end if

   allocate(Z(k1,k2),vp(n),vs(n),vw(n),ip(n),iq(n))
   x=b
   Z=A

   call GECP_Q(Z,k1,n,x,epsz,isw,ip,iq,vp,vs,vw,icon)
     if( icon /= 0 ) then
       write(IOUT,*) '*** LINEAR_SOLVE FAILED;'
       write(IOUT,*) '*** IT IS HIGHLY PROBABLE THAT THE MATRIX IS SINGULAR.'
       deallocate(Z,vp,vs,vw,ip,iq)
       return
     end if

! -------------- inverse power method ------------------- !
   if(present(sv)) then 
     allocate(y(n),u(n))
     y=1; y=y/sqrt(sum(y**2))
     if(present(delta)) deltaz=delta
     do k=1,KMAX
       u=y
       isw=2; call GECP_Q(Z,k1,n,y,epsz,isw,ip,iq,vp,vs,vw,icon)
       isw=4; call GECP_Q(Z,k1,n,y,epsz,isw,ip,iq,vp,vs,vw,icon)
       s=sqrt(sum(y**2))
       y=y/s
       if(maxval(abs(y-u))/maxval(abs(u))<deltaz) exit
       if(k==KMAX) then
        write(IOUT,*) '*** INVERSE ITERATION DID NOT CONVERGE UNTIL',k,'-TIMES'
        write(IOUT,*) '*** IN LINEAR_SOLVE'
       end if
     end do
     sv=1.0Q0/sqrt(s)
     deallocate(y,u)
   end if

   deallocate(Z,vp,vs,vw,ip,iq)
 end function linear_solve_q
! ==================================================================== !
! -------------------------------------------------------------------- !
!                                                                      !
         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
! ==================================================================== !
! -------------------------------------------------------------------- !
!                                                                      !
         SUBROUTINE GECP_D(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_D' SOLVES THE REAL LINEAR SYSTEM X=A B          !
!     BY GAUSSIAN ELIMINATION WITH COMPLETE PIVOTING AND SCALING.      !
!  <<PRECISION>>                                                       !
!     DOUBLE                                                           !
!  <<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.0D0)),DIMENSION(NX,NX),INTENT(INOUT) :: A
 REAL(KIND(1.0D0)),DIMENSION(NX),INTENT(INOUT)    :: B,VS,VP
 REAL(KIND(1.0D0)),DIMENSION(NX),INTENT(OUT)      :: VW
 REAL(KIND(1.0D0)),INTENT(IN)                     :: EPS
 INTEGER,DIMENSION(NX),INTENT(INOUT)              :: P,Q
 INTEGER,INTENT(OUT)                              :: ICON
!                                                                      !
 REAL(KIND(1.0D0))                                :: D,EPSZ
 INTEGER                                          :: I,J,K,L,C
!                                                                      !
 INTRINSIC EPSILON,MAX,ABS
!                                                                      !
! ==================================================================== !
!            CHECKING PARAMETERS                                       !
! -------------------------------------------------------------------- !
 ICON=0      
 IF(EPS<=0.0D0)  THEN 
   EPSZ=EPSILON(1.0D0)
 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.0D0
     DO J=1,N
       D=MAX(ABS(A(I,J)),D)
     END DO
     IF(D==0.0D0) 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.0D0
     DO I=1,N
       D=MAX(ABS(A(I,J)),D)
     END DO
     IF(D==0.0D0) 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.0D0
     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.0D0) 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.0D0/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_D
! ==================================================================== !
! -------------------------------------------------------------------- !
!                                                                      !
         SUBROUTINE GECP_Q(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_Q' SOLVES THE REAL LINEAR SYSTEM X=A B          !
!     BY GAUSSIAN ELIMINATION WITH COMPLETE PIVOTING AND SCALING.      !
!  <<PRECISION>>                                                       !
!     QUADRUPLE                                                        !
!  <<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.0Q0)),DIMENSION(NX,NX),INTENT(INOUT) :: A
 REAL(KIND(1.0Q0)),DIMENSION(NX),INTENT(INOUT)    :: B,VS,VP
 REAL(KIND(1.0Q0)),DIMENSION(NX),INTENT(OUT)      :: VW
 REAL(KIND(1.0Q0)),INTENT(IN)                     :: EPS
 INTEGER,DIMENSION(NX),INTENT(INOUT)              :: P,Q
 INTEGER,INTENT(OUT)                              :: ICON
!                                                                      !
 REAL(KIND(1.0Q0))                                :: D,EPSZ
 INTEGER                                          :: I,J,K,L,C
!                                                                      !
 INTRINSIC EPSILON,MAX,ABS
!                                                                      !
! ==================================================================== !
!            CHECKING PARAMETERS                                       !
! -------------------------------------------------------------------- !
 ICON=0      
 IF(EPS<=0.0Q0)  THEN 
   EPSZ=EPSILON(1.0Q0)
 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.0Q0
     DO J=1,N
       D=MAX(ABS(A(I,J)),D)
     END DO
     IF(D==0.0Q0) 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.0Q0
     DO I=1,N
       D=MAX(ABS(A(I,J)),D)
     END DO
     IF(D==0.0Q0) 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.0Q0
     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.0Q0) 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.0Q0/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_Q