!**********************************************************************!
!
! File: lineqp.f90 (01-Jun-2001)
!                  (08-Jun-2006)
!                  (29-Mar-2008)
!
! Loesung eines linearen Gleichungssystems
!
!   A*x = b
!
! mit dem Gauss'schen Eliminationsverfahren und Spaltenpivotsuche
!
! n.......Anzahl der Gleichungen
! A,b,x...Koeffizientenmatrix, rechte Seite und Loesungsvektor
!
!**********************************************************************!

module lineqp_m

!**********************************************************************!

  implicit none

  private

  integer,parameter,public::double=selected_real_kind(15)
  real(kind=double),parameter,public::eps=1.0e-08

  public::gaussp

contains

!**********************************************************************!

  subroutine gaussp(a_in,b,x)

!**********************************************************************!

    real(kind=double),dimension(:,:),intent(in)::a_in
    real(kind=double),dimension(:),intent(in)::b
    real(kind=double),dimension(:),intent(out)::x

    integer::i,itemp,j,k,m,n
    integer,dimension(size(b))::iperm
    real(kind=double)::pivot,pmax
    real(kind=double),dimension(size(b),size(b)+1)::a

    n=size(a_in,dim=1)   ! Dimension des Gleichungssystems

! Kopie der Koeffizientenmatrix a_in und der rechten Seite b
! in EINE n*(n+1) Matrix a (d.h. a_in wird nicht ueberschrieben)

    a(:,1:n)=a_in
    a(:,n+1)=b

! Permutationsvektor initialisieren

    iperm=(/(j,j=1,n)/)

    do i=1,n-1

! Maximales Pivotelement in Spalte i suchen

      m=maxloc(abs(a(iperm(i:n),i)),dim=1)+i-1
      pmax=abs(a(iperm(m),i))
      if(pmax<eps) then                        ! Matrix singulaer?
        write(unit=*,fmt="(a)") " pmax<eps"
        stop
      end if

! Zeilen (nur Indizes!) vertauschen

      if(m/=i) then
        itemp=iperm(i)
        iperm(i)=iperm(m)
        iperm(m)=itemp
      end if

! Elimination

      pivot=a(iperm(i),i)
      do k=i+1,n
       a(iperm(k),i+1:n+1)=a(iperm(k),i+1:n+1)- &
            a(iperm(i),i+1:n+1)*a(iperm(k),i)/pivot
      end do

    end do

! Ruecksubstitution

    if(abs(a(iperm(n),n))<eps) then
      write(unit=*,fmt="(a)") " |a(n,n)|<eps"
      stop
    end if

    x(n)=a(iperm(n),n+1)/a(iperm(n),n)
    do i=n-1,1,-1
      x(i)=(a(iperm(i),n+1)-sum(a(iperm(i),i+1:n)*x(i+1:n)))/a(iperm(i),i)
    end do

    return

  end subroutine gaussp

end module lineqp_m

!**********************************************************************!

program lineqp

!**********************************************************************!

  use lineqp_m

  implicit none

  integer::i,n
  real(kind=double),dimension(:),allocatable::b,x
  real(kind=double),dimension(:,:),allocatable::a

! Dimension des Gleichungssystems?

  write(unit=*,fmt="(a)",advance="no") " n="
  read(unit=*,fmt=*) n

  allocate(a(n,n),b(n),x(n))

! Koeffizientenmatrix und rechte Seite zeilenweise einlesen

  do i=1,n
    read(unit=*,fmt=*) a(i,:),b(i)
  end do

  call gaussp(a,b,x)

! Loesungsvektor ausgeben

  write(unit=*,fmt="(a,5(tr1,es14.7))") " x(:)=",x

end program lineqp

!**********************************************************************!