program test
	implicit double precision (a-h, o-z)
	parameter (n=1000)
	parameter (phi=0.5d0)
	dimension x(n,n), b(n), xinv(n,n), ipiv(n)
	real cpuStart, delta
	cpuStart = secnds(0.0)
	t = 0.0d0
	lwork = n
	do i = 1, n
	    ipiv(i) = i
		do j = 1, n
			x(i,j) = phi**iabs(i-j)
		enddo
	enddo
	call decomp(n, n, x, ipiv)
	do i = 1, n
		do j = 1, n
			b(j) = 0.0d0
			if (i .eq. j) b(j) = 1.0d0
		enddo
		call solve(n, n, x, b, ipiv)
		do j = 1, n
			xinv(i, j) = b(j)
		enddo
	enddo
	delta = secnds(cpuStart)
	tr = 0.0d0
	do i = 1, n
		tr = tr + xinv(i,i)
	enddo
	write(06, 1)tr
 1	format(" trace of inverse = ",1pe20.10)
	write(06, 2)delta
 2	format("  cpu time for matrix inverse = ", f10.3)
	end
c version: Feb 23 1992
C SOLVE.
C SOLUTION OF LINEAR SYSTEM  A*X = B .
C USING MATRIX TRIANGULARIZATION.
C
C ALGORITHM 423, LINEAR EQUATION SOLVER, BY CLEVE B. MOLER,
C                COLLECTED ALGORITHMS FROM THE ASSOCIATION FOR
C                COMPUTING MACHINERY
C
C THE NUMERICAL METHODS USED ARE DESCRIBED IN THE BOOK
C COMPUTER SOLUTION OF LINEAR ALGEBRAIC SYSTEMS, FORSYTHE, G.E. AND
C MOLER, C.B.  PRENTICE-HALL, ENGLEWOOD CLIFFS, N.J., 1967.
C
C ADDITIONAL REFERENCE - MATRIX COMPUTATIONS WITH FORTRAN AND PAGING,
C     COMM. ACM 15 (1972), 268-70
C
C**********************************************************************
      SUBROUTINE SOLVE(N, NDIM, A, B, IP)
C**********************************************************************
C
C INPUT PARAMETERS
C   N - ORDER OF MATRIX
C   NDIM - DECLARED DIMENSION OF ARRAY A
C   A - TRIANGULARIZED MATRIX OBTAINED FROM SUBROUTINE DECOMP
C   B - RIGHT-HAND SIDE VECTOR
C   IP - PIVOT VECTOR OBTAINED FROM SUBROUTINE DECOMP
C        (DO NOT USE IF DECOMP HAS SET IP(N) = 0)
C
C OUTPUT PARAMETERS
C   B - THE SOLUTION VECTOR X
C
      DOUBLE PRECISION A(NDIM, NDIM), B(NDIM), T
      INTEGER IP(NDIM)
C
      IF( N .NE. 1 ) THEN
        NM1 = N - 1
        DO K = 1, NM1
          KP1 = K + 1
          M = IP(K)
          T = B(M)
          B(M) = B(K)
          B(K) = T
          DO I = KP1, N
            B(I) = B(I) + A(I, K)*T
          ENDDO
        ENDDO
        DO KB = 1, NM1
          KM1 = N - KB
          K = KM1 + 1
          B(K) = B(K)/A(K, K)
          T =  - B(K)
          DO I = 1, KM1
            B(I) = B(I) + A(I, K)*T
          ENDDO
        ENDDO
      ENDIF
      B(1) = B(1)/A(1, 1)
      RETURN
      END
C DECOMP.
C MATRIX TRIAANGULARIZATION BY GAUSSIAN ELIMINATION.
C
C**********************************************************************
      SUBROUTINE DECOMP(N, NDIM, A, IP)
C**********************************************************************
C INPUT PARAMETERS
C   N - ORDER OF MATRIX
C   NDIM - DECLARED DIMENSION OF THE ARRAY A
C   A - MATRIX TO BE TRIANGULARIZED
C
C OUTPUT PARAMETERS
C   A(I,J),  I.LE.J,  UPPER TRIANGULAR FACTOR, U
C   A(I,J),  I.GT.J,  LOWER TRIANGULAR FACTOR, 1-L (MULTIPIERS)
C   IP(K),  K.LT.N,  INDEX OF THE K-TH PIVOT ROW.
C   IP(N) = (-1)**(NUMBER OF INTERCHANGES)  OR 0
C           IF IP(N)=0,  MATRIX A IS SINGULAR AND SUBROUTINE SOLVE
C           WILL DIVIDE BY ZERO
C
C INTERCHANGES FINISHED IN U, ONLY PARTLY IN L
C
      DOUBLE PRECISION A(NDIM, NDIM), T, ZERO, DABS
      INTEGER IP(NDIM)
C INITIALIZE NUMERICAL CONSTANT
      DATA ZERO/0.0D0/
C
      IP(N) = 1
      DO K = 1, N
        IF( K .NE. N ) THEN
          KP1 = K + 1
          M = K
          DO I = KP1, N
            IF( DABS(A(I, K)) .GT. DABS(A(M, K)) ) M = I
          ENDDO
          IP(K) = M
          IF( M .NE. K ) IP(N) =  - IP(N)
          T = A(M, K)
          A(M, K) = A(K, K)
          A(K, K) = T
          IF( T .NE. ZERO ) THEN
            DO I = KP1, N
              A(I, K) =  - A(I, K)/T
            ENDDO
            DO J = KP1, N
              T = A(M, J)
              A(M, J) = A(K, J)
              A(K, J) = T
              IF( T .NE. ZERO ) THEN
                DO I = KP1, N
                  A(I, J) = A(I, J) + A(I, K)*T
                ENDDO
              ENDIF
            ENDDO
          ENDIF
        ENDIF
        IF( A(K, K) .EQ. ZERO ) IP(N) = 0
      ENDDO
      RETURN
      END