SUBROUTINE gaussj(a,b)
	USE nrtype; USE nrutil, ONLY : assert_eq,nrerror,outerand,outerprod,swap
	IMPLICIT NONE
	REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a,b
	INTEGER(I4B), DIMENSION(size(a,1)) :: ipiv,indxr,indxc
	LOGICAL(LGT), DIMENSION(size(a,1)) :: lpiv
	REAL(SP) :: pivinv
	REAL(SP), DIMENSION(size(a,1)) :: dumc
	INTEGER(I4B), TARGET :: irc(2)
	INTEGER(I4B) :: i,l,n
	INTEGER(I4B), POINTER :: irow,icol
	n=assert_eq(size(a,1),size(a,2),size(b,1),'gaussj')
	irow => irc(1)
	icol => irc(2)
	ipiv=0
	do i=1,n
		lpiv = (ipiv == 0)
		irc=maxloc(abs(a),outerand(lpiv,lpiv))
		ipiv(icol)=ipiv(icol)+1
		if (ipiv(icol) > 1) call nrerror('gaussj: singular matrix (1)')
		if (irow /= icol) then
			call swap(a(irow,:),a(icol,:))
			call swap(b(irow,:),b(icol,:))
		end if
		indxr(i)=irow
		indxc(i)=icol
		if (a(icol,icol) == 0.0) &
			call nrerror('gaussj: singular matrix (2)')
		pivinv=1.0_sp/a(icol,icol)
		a(icol,icol)=1.0
		a(icol,:)=a(icol,:)*pivinv
		b(icol,:)=b(icol,:)*pivinv
		dumc=a(:,icol)
		a(:,icol)=0.0
		a(icol,icol)=pivinv
		a(1:icol-1,:)=a(1:icol-1,:)-outerprod(dumc(1:icol-1),a(icol,:))
		b(1:icol-1,:)=b(1:icol-1,:)-outerprod(dumc(1:icol-1),b(icol,:))
		a(icol+1:,:)=a(icol+1:,:)-outerprod(dumc(icol+1:),a(icol,:))
		b(icol+1:,:)=b(icol+1:,:)-outerprod(dumc(icol+1:),b(icol,:))
	end do
	do l=n,1,-1
		call swap(a(:,indxr(l)),a(:,indxc(l)))
	end do
	END SUBROUTINE gaussj