Dinv_gj Subroutine

subroutine Dinv_gj(a)
  real(8), dimension(:,:), intent(inout) :: a
  integer, dimension(size(a,1))      :: ipiv,indxr,indxc
  !these arrays are used for bookkeeping on the pivoting.
  !integer                            :: nn
  logical, dimension(size(a,1))      :: lpiv
  real(8)                            :: pivinv
  real(8), dimension(size(a,1))      :: dumc
  integer, target                    :: irc(2)
  integer                            :: i,l,n
  integer, pointer                   :: irow,icol
  real(8)                            :: zero=0.d0, one=1.d0
  n=size(a,1)
  irow => irc(1)
  icol => irc(2)
  ipiv=0
  do i=1,n
     !main loop over columns to be reduced.
     lpiv = (ipiv == 0)
     !begin search for a pivot element.
     irc=maxloc(abs(a),outerand(lpiv,lpiv))
     ipiv(icol)=ipiv(icol)+1
     if (ipiv(icol) > 1) stop 'gaussj:singular matrix (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 inverse matrix will be
     !scrambled by
     !columns.
     if (irow /= icol) call swap(a(irow,:),a(icol,:))
     indxr(i)=irow !we are now ready to divide the pivot row by the pivot element,
     !located at irow and icol.
     indxc(i)=icol
     if (a(icol,icol) == zero) stop 'gaussj:singular matrix (2)'
     pivinv=one/a(icol,icol)
     a(icol,icol)= one !cmplx(one,zero)
     a(icol,:)=a(icol,:)*pivinv
     dumc=a(:,icol)
     !next, we reduce the rows, except for the pivot one, of course.
     a(:,icol)     = zero !cmplx
     a(icol,icol)  = pivinv
     a(1:icol-1,:) = a(1:icol-1,:) - outerprod(dumc(1:icol-1),a(icol,:))
     a(icol+1:,:)  = a(icol+1:,:)  - outerprod(dumc(icol+1:),a(icol,:))
  end do
  !it only remains to unscramble the solution in view of the column
  !interchanges.
  !we do this by interchanging pairs of columns in the reverse order that the
  !permutation
  !was built up.
  do l=n,1,-1
     call swap(a(:,indxr(l)),a(:,indxc(l)))
  end do
contains
  function outerand(a,b)
    implicit none
    logical, dimension(:), intent(in)   :: a,b
    logical, dimension(size(a),size(b)) :: outerand
    outerand = spread(a,dim=2,ncopies=size(b)).and.spread(b,dim=1,ncopies=size(a))
  end function outerand
  function outerprod(a,b)
    real(8), dimension(:), intent(in) :: a,b
    real(8), dimension(size(a),size(b)) :: outerprod
    outerprod = spread(a,dim=2,ncopies=size(b)) * &
         spread(b,dim=1,ncopies=size(a))
  end function outerprod
end subroutine Dinv_gj