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