C+
C NAME:
C	MATRIX_INV_R8
C PURPOSE:
C	Matrix inversion
C CATEGORY:
C	Math: matrix algebra
C CALLING SEQUENCE:
	function MATRIX_INV_R8(N,A,IK,JK)
C INPUTS:
C	N	integer			dimension of the NxN matrix A
C	A	double precision	matrix
C	IK,JK	integer			scratch arrays
C OUTPUTS:
C	MATRIX_INV_R8  integer		0 if inversion not succesfull; 1 if succesfull
C	A	double precision	inverse matrix
C PROCEDURE:
C	Matrix inversion by Gauss-Jordan elimination with full pivoting
C MODIFICATION HISTORY:
C	Joop Gronenschild (SRON)
C-
	    integer		N
	    double precision	A(N,N)
	    integer		IK(N)
	    integer		JK(N)

	integer			EXPO
	double precision	EMAX
	double precision	FACTOR
	double precision	AMAX
	double precision	ADUM

	EMAX = 0.
	do I=1,N
	    do J=1,N
		if (abs(A(I,J)) .gt. abs(EMAX)) EMAX  = A(I,J)
	    end do
	end do
	EMAX = abs(EMAX)	!Largest element in absolute value
	EXPO = 0
	do while (EMAX .lt. 1.) 
	    EXPO = EXPO-1
	    EMAX = EMAX*10
	end do
	do while (EMAX .ge. 10.)
	    EXPO = EXPO+1
	    EMAX = EMAX*.1
	end do
	RMACH = 1E-16
	if (EXPO .ne. 0) then
	    FACTOR = 1.
		if (EXPO .gt. 0) then
		do I=1,EXPO
		    FACTOR = FACTOR*10.
		end do
	    else
		do I=1,-EXPO
		    FACTOR = FACTOR/10.
		end do
	    end if
	    RMACH = RMACH*FACTOR
	end if

	MATRIX_INV_R8 = 0		!Suppose there is no inverse
	do K=1,N
	    AMAX = 0.	!Find the largest element A(I,J) in rest of matrix
	    do I=K,N
		do J=K,N
		    if (abs(A(I,J)) .gt. abs(AMAX)) then
			AMAX  = A(I,J)
			IK(K) = I	!Row index largest element
			JK(K) = J	!Column index largest element
		    end if
		end do
	    end do
	    if (abs(AMAX) .lt. RMACH) return	!Singular matrix; no inverse
	    I = IK(K)
	    if (I .ne. K) then		!If largest element not in row K ...
		do J=1,N		!... interchange rows I and K
		    ADUM   = A(K,J)
		    A(K,J) = A(I,J)
		    A(I,J) = -ADUM
		end do
	    end if
	    J = JK(K)
	    if (J .ne. K) then		!If largest element not in column K ...
		do I=1,N		!... interchange columns J and K
		    ADUM   = A(I,K)
		    A(I,K) = A(I,J)
		    A(I,J) = -ADUM
		end do
	    end if			!Now AMAX is in position (K,K)
	    do I=1,N	!Accumulate elements of inverse matrix in column K
		A(I,K) = -A(I,K)/AMAX	!Mult. factors for cleaning column K
	    end do
	    do I=1,N
		if (I .ne. K) then	!Subtract multiple of ...
		    do J=1,N		!. row K from all other rows
			if (J .ne. K) A(I,J)= A(I,J)+A(I,K)*A(K,J)
		    end do
		end if
	    end do
	    do J=1,N
		A(K,J) = A(K,J)/AMAX !Divide row K by AMAX to get 1 on diagonal
	    end do
	    A(K,K) = 1./AMAX
	end do

	do K=N,1,-1		!Restore ordering of matrix
	    J = IK(K)		!Run through sequence of row and column ...
	    if (J .gt. K) then	!. permutations in reverse order.
		do I=1,N	!Previous row changes IK become column changes
		    ADUM   = A(I,K)
		    A(I,K) = -A(I,J)
		    A(I,J) = ADUM
		end do
	    end if
	    I = JK(K)		!Column changes JK become row changes
	    if (I .gt. K) then
	        do J=1,N
		    ADUM   = A(K,J)
		    A(K,J) = -A(I,J)
		    A(I,J) = ADUM
		end do
	    end if
	end do
	MATRIX_INV_R8 = 1		!Inverse matrix exists
	return
	end