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