gaussji Function

  function gaussji(c,d,n) result(a)
    !Linear equation solution by Gauss-Jordan elimination, used to find inverse of matrix c.
    !possibly slow for large "c" (LAPACK better?)
    !c needs to be square and have dimension n
    !c(1:n,1:n) is an input matrix stored in an array of physical dimensions n by n.
    !d(1:n,1:n) is an input matrix containing the n by n identity matrix.
    !On  output, a(1:n,1:n) (and b(1:n,1:n)) are the inverse of c
    !Parameter: NMAX is  the  largest  anticipated  value  of n.

    integer :: n
    real, intent(in) :: c(n,n) !WILL BE OVERWRITTEN!!
    real, intent(in) :: d(n, n)
    real :: a(n,n),b(n,n)
    integer, parameter :: NMAX = 50
    integer :: m, i, icol, irow, j, k, l, ll, indxc(NMAX), indxr(NMAX), ipiv(NMAX)
    !The integer arrays ipiv, indxr, and indxc are  used for bookkeeping  on the pivoting.
    REAL :: big, dum, pivinv
    a=c
    b=d
    m=n
    do j = 1, n
      ipiv(j) = 0
    end do
    do i = 1, n !This  is  the  main  loop  over  the  columns  to  be  reduced.
      big = 0.
      do j = 1, n !This  is  the  outer  loop  of  the  search  for  a  pivot  element.
        if (ipiv(j) .ne. 1) then
          do k = 1, n
            if (ipiv(k) .eq. 0) then
              if (abs(a(j, k)) .ge. big) then
                big = abs(a(j, k))
                irow = j
                icol = k
              endif
              !else if (ipiv(k).gt.1) then
              !pause 'singular matrix in gaussj'
            end if
          end do
        end if
      end do
      ipiv(icol) = ipiv(icol) + 1
      !We  now  have  the  pivot  element,  so  we  interchange  rows,  if  needed,  to  put  the  pivot
      !element  on  the  diagonal.  The  columns  are  not  physically  interchanged,  only  relabeled:
      !indxc(i), the column of the ith pivot element, is the ith column that is reduced, while
      !indxr(i) is  the  row in  which  that  pivot  element  was  originally  located.  If
      !indxr(i) /= indxc(i) there  is  an  implied  column  interchange.  With  this  form  of  bookkeeping,  the
      !solution b's  will  end  up  in  the  correct  order,  and  the  inverse  matrix  will  be  scrambled by  columns
      if (irow .ne. icol) then
        do l = 1, n
          dum = a(irow, l)
          a(irow, l) = a(icol, l)
          a(icol, l) = dum
        end do
        do l = 1, m
          dum = b(irow, l)
          b(irow, l) = b(icol, l)
          b(icol, l) = dum
        enddo
      endif
      !We are now ready to divide the pivot row by the pivot element, located at irow and icol.
      indxr(i) = irow
      indxc(i) = icol
      !if (a(icol,icol).eq.0.) pause 'singular matrix in gaussj'
      pivinv = 1./a(icol, icol)
      a(icol, icol) = 1.
      do l = 1, n
        a(icol, l) = a(icol, l)*pivinv
      end do
      do l = 1, m
        b(icol, l) = b(icol, l)*pivinv
      end do
      do ll = 1, n
        !Next,  we  reduce  the  rows, except for the  pivot  one, of course.
        if (ll .ne. icol) then
          dum = a(ll, icol)
          a(ll, icol) = 0.
          do l = 1, n
            a(ll, l) = a(ll, l) - a(icol, l)*dum
          end do
          do l = 1, m
            b(ll, l) = b(ll, l) - b(icol, l)*dum
          end do
        end if
      end do
    end do
    !This is the end of the main loop over columns of the reduction.
    do l = n, 1, -1
      !It  only  remains  to  unscramble  the  solution  in  view
      !of  the  column  interchanges.  We  do  this  by  in-
      !terchanging pairs of columns in the reverse order
      !that the permutation was built  up.
      if (indxr(l) .ne. indxc(l)) then
        do k = 1, n
          dum = a(k, indxr(l))
          a(k, indxr(l)) = a(k, indxc(l))
          a(k, indxc(l)) = dum
        end do
      end if
    end do
    return
    !And  we  are  done.
  end function gaussji