!**********************************************************************!
!
! File: lineqp.f90 (01-Jun-2001)
!                  (02-Jun-2005)
!
! 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 lineq_m

  implicit none

  public::gaussp

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

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(1)::mvect
    integer,dimension(size(b))::iperm
    real(kind=double)::pivot,pmax
    real(kind=double),dimension(size(b),size(b)+1)::a

    n=size(b)   ! 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

      mvect=maxloc(abs(a(iperm(i:n),i)))
      m=mvect(1)+i-1
      pmax=abs(a(iperm(m),i))
      if(pmax<=eps) then                       ! Matrix singulaer?
        write(unit=*,fmt="("" 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(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 lineq_m

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

program lineqp

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

  use lineq_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="("" n="")",advance="no")
  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="("" x(:)="",5(tr1,es14.7))") x

end program lineqp

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