      subroutine mat_gau(amat,n,bmat,m) ! bernd berg aug 3 2004.
c invert square matrix amat by gauss-jordan elimination.
c subsequently calculation of solution vector for m linear eqns.
c for m=0 only inverse matrix calculated.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(nmax=100)
      dimension amat(n,n),bmat(n,m),bin(nmax),ik(nmax),jk(nmax)
      if(n.gt.nmax) stop "mat_gau: enlarge nmax (n.gt.nmax)!"
      ltest=.true.
      ltest=.false.
c - find largest element.
      do k=1,n
        amax=zero
1       do i=k,n
          do j=k,n
            if(abs(amat(i,j)).gt.abs(amax)) then
              amax=amat(i,j)
              ik(k)=i
              jk(k)=j
            end if
          end do ! do j
        end do ! do i
        if(amax.eq.zero) stop "mat_gau: singular matrix."
        det=one
c - interchange rows and columns to put amax in amat(k,k).
        i=ik(k)
        if(i.lt.k) goto 1
        if(i.gt.k) then
          do j=1,n
            aa=amat(k,j)
            amat(k,j)=amat(i,j)
            amat(i,j)=-aa
          end do ! do j
        endif ! if i.gt.k.
        j=jk(k)
        if(j.lt.k) goto 2
        if(j.gt.k) then
          do i=1,n
            aa=amat(i,k)
            amat(i,k)=amat(i,j)
            amat(i,j)=-aa
          end do ! do i.
        endif ! if j.gt.k.
c - accumulate elements of the inverse matrix.
        do i=1,n
          if(i.ne.k) amat(i,k)=-amat(i,k)/amax
        end do ! do i.
        do i=1,n
          do j=1,n
          if((i.ne.k).and.(j.ne.k)) 
     1        amat(i,j)=amat(i,j)+amat(i,k)*amat(k,j)
          end do ! do j.
        end do ! do i.
        do j=1,n
          if(j.ne.k) amat(k,j)=amat(k,j)/amax
        end do ! do j.
        amat(k,k)=one/amax
        det=det*amax ! determinant (presently not returned on output).
      end do ! do k.
2     continue  
      if(ltest) stop "mat_gau - ltest."
c - restore ordering of the matrix.
      do il=1,n
        k=n+1-il
        j=ik(k)
        if(j.gt.k) then
          do i=1,n
            aa=amat(i,k)
            amat(i,k)=-amat(i,j)
            amat(i,j)=aa
          end do ! do i.
        endif ! if j.gt.k.
        i=jk(k)
        if(i.gt.k) then
          do j=1,n
            aa=amat(k,j)
            amat(k,j)=-amat(i,j)
            amat(i,j)=aa
          end do ! do j.
        endif ! if i.gt.k.
      end do ! do il.
c - calculate the solution vector.
      if(m.eq.0) return
      do k=1,m
        do i=1,n
          bin(i)=bmat(i,k)
        end do ! do i.
        do i=1,n
          bmat(i,k)=zero
          do j=1,n
            bmat(i,k)=bmat(i,k)+amat(i,j)*bin(j)
          end do ! do j.
        end do ! do i.
      end do ! do k.

      return
      end
