      SUBROUTINE PrintVector(v,n)
      INCLUDE 'stdhdr.for'
      REAL n, v(0:maxv)
      INTEGER i
      DO i=0,n-1
        PRINT *,v(i)
      END DO
      END


      SUBROUTINE GaussSwap (Swapx ,Swapy )
      REAL Swapx, Swapy,temp

	   temp = Swapx
	   Swapx = Swapy
	   Swapy = temp
      RETURN
      END




      SUBROUTINE GaussJordan( coefary,  constary,numcol,
     +		solcoef, invary, det)

       INCLUDE 'stdhdr.for'
       REAL coefary(0:maxc, 0:maxc), constary(0:maxc )
       REAL det, solcoef( 0:maxc), invary(0:maxc,0:maxc)
       REAL leval, piv, t
       INTEGER  numcol, i, j, lerow, lecol, k, l, iErr
       INTEGER pivlst[ALLOCATABLE](:,:)
       LOGICAL pivchk[ALLOCATABLE](:)

       ALLOCATE(pivlst(0:maxc,0:1),STAT=iErr)
       CALL CheckMem(iErr)
       ALLOCATE(pivchk(0:maxc),STAT=iErr)
       CALL CheckMem(iErr)

       det = 1.0
       DO i  = 0, numcol  - 1
	pivchk(i) = .FALSE.
	 DO j  = 0, numcol  - 1
	   invary(i, j ) = coefary(i, j )
	 END DO
       END DO

       DO i  = 0, numcol  - 1
	  leval = 0.0
	  DO j	= 0, numcol  - 1
	    IF ( .not.(pivchk(j) )) THEN
	     DO k  = 0, numcol	- 1
	       IF (.not.(pivchk(k))) THEN
		   IF (ABS(invary(j, k )) .GT. leval) THEN
		      lerow  = j
		      lecol  = k
		      leval = ABS(invary(j, k ))
		    END IF
	       END IF
	     END DO
	    END IF
	END DO

	 pivchk(lecol ) = .TRUE.
	 pivlst(i, 0) = lerow
	 pivlst(i, 1) = lecol

	 IF (lerow  .NE. lecol ) THEN
	   det = -det
	   DO l  = 0, numcol  - 1
	     call GaussSwap( invary(lerow, l ), invary(lecol, l ))
	   END DO
	   call GaussSwap (constary(lerow ), constary(lecol ))
	 END IF

	 piv = invary(lecol, lecol )
	 det = det * piv
	 IF (det .GT. 1.0E30) det = 1.0E30 !prevents numeric overflow
	 invary(lecol, lecol ) = 1.0
	 DO l  = 0, numcol  - 1
	   invary(lecol, l ) = invary(lecol, l ) / piv
	 END DO
	 constary(lecol ) = constary(lecol ) / piv

	 DO l1	= 0, numcol  - 1
	   IF (l1  .NE. lecol ) THEN
	     t = invary(l1, lecol )
	     invary(l1, lecol ) = 0.0
	     DO l  = 0, numcol	- 1
	       invary(l1, l ) = invary(l1, l ) - invary(lecol, l ) * t
	     END DO
	     constary(l1 ) = constary(l1 ) - constary(lecol ) * t
	   END IF
	 END DO
       END DO

      DO i  = 0, numcol  - 1
	 l  = numcol  - i  - 1
	 IF (pivlst(l, 0) .NE. pivlst(l, 1)) THEN
	   lerow  = pivlst(l, 0)
	   lecol  = pivlst(l, 1)
	   DO k  = 0, numcol  - 1
	     CALL GaussSwap (invary(k, lerow ), invary(k, lecol ))
	   END DO
	  END IF
      END DO

      DO i  = 0, numcol  - 1
	  solcoef(i ) = constary(i )
      END DO
      DEALLOCATE(pivlst,STAT=iErr)
      CALL CheckMem(iErr)
      DEALLOCATE(pivchk,STAT=iErr)
      CALL CheckMem(iErr)

      END  !GaussJordan


