14 Replies - 4188 Views - Last Post: 10 February 2011 - 08:05 AM

#1 Louisda16th  Icon User is offline

  • dream.in.assembly.code
  • member icon

Reputation: 15
  • View blog
  • Posts: 1,967
  • Joined: 03-August 06

FORTRAN 90: Matrix Inversion

Posted 17 September 2007 - 08:09 AM

Description: Use the subroutine like this: !Example program PROGRAM INVERSE IMPLICIT NONE REAL, ALLOCATABLE, DIMENSION(:,:) :: Matrix, invMatrix INTEGER :: i, j, n, ErrorFlag PRINT*,"Enter size of square matrix" READ*,n ALLOCATE(Matrix(n,n)) ALLOCATE(invMatrix(n,n)) PRINT*,"Enter elements of matrix Row wise" DO i = 1, n READ*,(Matrix(i,j), j = 1, n) END DO CALL FindInv(Matrix, invMatrix, n, ErrorFlag) Print*,"Inverse Matrix:" DO i = 1, n PRINT*,(invMatrix(i,j), j = 1, n) END DO END PROGRAM INVERSE !Subroutine comes hereThis subroutine finds the inverse of a matrix
!Subroutine to find the inverse of a square matrix
!Author : Louisda16th a.k.a Ashwith J. Rego
!Reference : Algorithm has been well explained in:
!http://math.uww.edu/~mcfarlat/inverse.htm           
!http://www.tutor.ms.unimelb.edu.au/matrix/matrix_inverse.html
SUBROUTINE FINDInv(matrix, inverse, n, errorflag)
	IMPLICIT NONE
	!Declarations
	INTEGER, INTENT(IN) :: n
	INTEGER, INTENT(OUT) :: errorflag  !Return error status. -1 for error, 0 for normal
	REAL, INTENT(IN), DIMENSION(n,n) :: matrix  !Input matrix
	REAL, INTENT(OUT), DIMENSION(n,n) :: inverse !Inverted matrix
	
	LOGICAL :: FLAG = .TRUE.
	INTEGER :: i, j, k, l
	REAL :: m
	REAL, DIMENSION(n,2*n) :: augmatrix !augmented matrix
	
	!Augment input matrix with an identity matrix
	DO i = 1, n
		DO j = 1, 2*n
			IF (j <= n ) THEN
				augmatrix(i,j) = matrix(i,j)
			ELSE IF ((i+n) == j) THEN
				augmatrix(i,j) = 1
			Else
				augmatrix(i,j) = 0
			ENDIF
		END DO
	END DO
	
	!Reduce augmented matrix to upper traingular form
	DO k =1, n-1
		IF (augmatrix(k,k) == 0) THEN
			FLAG = .FALSE.
			DO i = k+1, n
				IF (augmatrix(i,k) /= 0) THEN
					DO j = 1,2*n
						augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j)
					END DO
					FLAG = .TRUE.
					EXIT
				ENDIF
				IF (FLAG .EQV. .FALSE.) THEN
					PRINT*, "Matrix is non - invertible"
					inverse = 0
					errorflag = -1
					return
				ENDIF
			END DO
		ENDIF
		DO j = k+1, n			
			m = augmatrix(j,k)/augmatrix(k,k)
			DO i = k, 2*n
				augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i)
			END DO
		END DO
	END DO
	
	!Test for invertibility
	DO i = 1, n
		IF (augmatrix(i,i) == 0) THEN
			PRINT*, "Matrix is non - invertible"
			inverse = 0
			errorflag = -1
			return
		ENDIF
	END DO
	
	!Make diagonal elements as 1
	DO i = 1 , n
		m = augmatrix(i,i)
		DO j = i , (2 * n)				
			   augmatrix(i,j) = (augmatrix(i,j) / m)
		END DO
	END DO
	
	!Reduced right side half of augmented matrix to identity matrix
	DO k = n-1, 1, -1
		DO i =1, k
		m = augmatrix(i,k+1)
			DO j = k, (2*n)
				augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m
			END DO
		END DO
	END DO				
	
	!store answer
	DO i =1, n
		DO j = 1, n
			inverse(i,j) = augmatrix(i,j+n)
		END DO
	END DO
	errorflag = 0
END SUBROUTINE FINDinv


Is This A Good Question/Topic? 0
  • +

Replies To: FORTRAN 90: Matrix Inversion

#2 VAR1318  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 24-February 09

Re: FORTRAN 90: Matrix Inversion

Posted 24 February 2009 - 01:55 PM

!ADDED IDENTITY MATRIX CHECK & !FORMATED FOR MY SIMPLE FORTRAN PROGRAM INVERSE IMPLICIT NONE REAL, ALLOCATABLE, DIMENSION(:,:) :: Matrix, invMatrix REAL, ALLOCATABLE, DIMENSION(:,:) :: IDENTITY INTEGER :: i, j, n, ErrorFlag ,K REAL SUM PRINT*,"Enter size of square matrix" READ*,n ALLOCATE(Matrix(n,n)) ALLOCATE(invMatrix(n,n)) ALLOCATE(IDENTITY(n,n)) PRINT*,"Enter elements of matrix Row wise" DO i = 1, n READ*,(Matrix(i,j), j = 1, n) END DO CALL FindInv(Matrix, invMatrix, n, ErrorFlag) Print*,"Inverse Matrix:" DO i = 1, n PRINT*,(invMatrix(i,j), j = 1, n) END DO C OPTIONAL CHECK, COMMENT OUT IF NOT WANTED DO I=1,N DO K=1,N SUM=0 DO J=1,N SUM=SUM+MATRIX(I,J)*INVMATRIX(J,K) ENDDO ! LOOP J IDENTITY(I,K)=SUM ENDDO ! LOOP K ENDDO ! LOOP I PRINT* Print*,"IDENTITY Matrix:" DO i = 1, n PRINT*,(IDENTITY(i,j), j = 1, n) END DO STOP END PROGRAM INVERSE !Subroutine to find the inverse of a square matrix !Author : Louisda16th a.k.a Ashwith J. Rego !Reference : Algorithm has been well explained in: !http://math.uww.edu/~mcfarlat/inverse.htm !http://www.tutor.ms.unimelb.edu.au/matrix/matrix_inverse.html SUBROUTINE FINDInv(matrix, inverse, n, errorflag) IMPLICIT NONE !Declarations INTEGER, INTENT(IN) :: n INTEGER, INTENT(OUT) :: errorflag !Return error status. -1 for error, 0 for normal REAL, INTENT(IN), DIMENSION(n,n) :: matrix !Input matrix REAL, INTENT(OUT), DIMENSION(n,n) :: inverse !Inverted matrix LOGICAL :: FLAG = .TRUE. INTEGER :: i, j, k, l REAL :: m REAL, DIMENSION(n,2*n) :: augmatrix !augmented matrix !Augment input matrix with an identity matrix DO i = 1, n DO j = 1, 2*n IF (j <= n ) THEN augmatrix(i,j) = matrix(i,j) ELSE IF ((i+n) == j) THEN augmatrix(i,j) = 1 Else augmatrix(i,j) = 0 ENDIF END DO END DO !Reduce augmented matrix to upper traingular form DO 101 k =1, n-1 IF (augmatrix(k,k) == 0) THEN ! LABEL=201 FLAG = .FALSE. DO 102 i = k+1, n IF (augmatrix(i,k) /= 0) THEN DO j = 1,2*n augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j) ENDDO FLAG = .TRUE. EXIT ENDIF IF (FLAG .EQV. .FALSE.) THEN PRINT*, "Matrix is non - invertible" inverse = 0 errorflag = -1 return ENDIF 102 CONTINUE ENDIF ! LABEL=201 DO j = k+1, n m = augmatrix(j,k)/augmatrix(k,k) DO i = k, 2*n augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i) END DO END DO 101 CONTINUE !END REDUCE TO UPPER TRIANGULAR FORM !Test for invertibility DO i = 1, n IF (augmatrix(i,i) == 0) THEN PRINT*, "Matrix is non - invertible" inverse = 0 errorflag = -1 return ENDIF END DO !Make diagonal elements as 1 DO i = 1 , n m = augmatrix(i,i) DO j = i , (2 * n) augmatrix(i,j) = (augmatrix(i,j) / m) END DO END DO !Reduced right side half of augmented matrix to identity matrix DO k = n-1, 1, -1 DO i =1, k m = augmatrix(i,k+1) DO j = k, (2*n) augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m END DO END DO END DO !store answer DO i =1, n DO j = 1, n inverse(i,j) = augmatrix(i,j+n) END DO END DO errorflag = 0 RETURN END C END SUBROUTINE FINDinv
Was This Post Helpful? 0
  • +
  • -

#3 masaif2009  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 23-June 09

Re: FORTRAN 90: Matrix Inversion

Posted 03 August 2009 - 09:27 PM

I try to use this program to find the inverse matrix of matrix A =10 31 5 12 35 5 5 15 3 The inverse of A = 0.2562 0.1563 -0,6875 0.6875 -0.3125 -0.6250 -0.3125 -0.3125 1.3750 To test the result A*A-1(inverse gives identity) 1.0 3.0 0.0 -0.0 1.0 0.0 0.0 0.0 1.0 In the identity matrix 3.0 must be 0.0 This fortran to find an inverse of a matrix gives wrong result. Please advise what to do!!!
Was This Post Helpful? 0
  • +
  • -

#4 Louisda16th  Icon User is offline

  • dream.in.assembly.code
  • member icon

Reputation: 15
  • View blog
  • Posts: 1,967
  • Joined: 03-August 06

Re: FORTRAN 90: Matrix Inversion

Posted 03 August 2009 - 10:18 PM

I got the correct result when I ran it with your example. Which compiler are you using?
Was This Post Helpful? 0
  • +
  • -

#5 rhamdan  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 17-September 09

Re: FORTRAN 90: Matrix Inversion

Posted 17 September 2009 - 10:30 AM

There might be a bug in the code not sure please check the condition IF (FLAG .EQV. .FALSE.) should be outside the loop Do i =k+1,n
Was This Post Helpful? 0
  • +
  • -

#6 Louisda16th  Icon User is offline

  • dream.in.assembly.code
  • member icon

Reputation: 15
  • View blog
  • Posts: 1,967
  • Joined: 03-August 06

Re: FORTRAN 90: Matrix Inversion

Posted 17 September 2009 - 10:40 AM

Not really. That block of code converts the matrix to upper-triangular form. FLAG is initially set to False if an element in the diagonal is zero. After this row operations are performed. If the diagonal element is still a zero then Flag remains false and the matrix is non-invertible. The condition cannot be used outside the loop. The variable Flag is used only in this loop and nowhere else in the procedure
Was This Post Helpful? 0
  • +
  • -

#7 Chris.Webster  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 14-October 09

Re: FORTRAN 90: Matrix Inversion

Posted 14 October 2009 - 04:29 AM

Brilliant. There are some very clever people about. I wish that I were one of them. Chris Webster
Was This Post Helpful? 0
  • +
  • -

#8 vb10000  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 09-January 10

Re: FORTRAN 90: Matrix Inversion

Posted 09 January 2010 - 10:07 AM

I tried to invert the following matrix: 0 1 0 0 1 1 1 1 1 and got an error. I had to change the matrix to 1 1 1 0 1 1 0 1 0 to get the routine to work.
Was This Post Helpful? 0
  • +
  • -

#9 vb10000  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 09-January 10

Re: FORTRAN 90: Matrix Inversion

Posted 09 January 2010 - 10:09 AM

Sorry, the matrix did not show well. First 3 #s are 1st row, 2nd set is 2nd row and 3rd is the third row. I am using Salford fortran compiler. Thanks.
Was This Post Helpful? 0
  • +
  • -

#10 Louisda16th  Icon User is offline

  • dream.in.assembly.code
  • member icon

Reputation: 15
  • View blog
  • Posts: 1,967
  • Joined: 03-August 06

Re: FORTRAN 90: Matrix Inversion

Posted 13 January 2010 - 10:45 AM

Thanks for pointing that out. I'll look into it.
Was This Post Helpful? 0
  • +
  • -

#11 Meister_Chicho  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 13
  • Joined: 30-March 10

Re: FORTRAN 90: Matrix Inversion

Posted 30 March 2010 - 02:19 AM

Thanks a lot, i'm going to give it a try inmediately
Was This Post Helpful? 0
  • +
  • -

#12 bfgreen  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 27-October 10

Re: FORTRAN 90: Matrix Inversion

Posted 27 October 2010 - 07:12 AM

I compared the results of the code with some other subroutines for inversion in Compaq VF. Definitely , it does not work well. Please verify it.
Was This Post Helpful? 0
  • +
  • -

#13 mitek  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 10-February 11

Re: FORTRAN 90: Matrix Inversion

Posted 10 February 2011 - 08:02 AM

There are errors concerning comparison of float with zero, and position of one of the conditions. The corrected code is: !Subroutine to find the inverse of a square matrix !Author : Louisda16th a.k.a Ashwith J. Rego !Reference : Algorithm has been well explained in: !http://math.uww.edu/~mcfarlat/inverse.htm !http://www.tutor.ms.unimelb.edu.au/matrix/matrix_inverse.html !Code was corrected by Dmitry Svinkin at Ioffe Inst SUBROUTINE FINDInv(matrix, inverse, n, errorflag) IMPLICIT NONE !Declarations INTEGER, INTENT(IN) :: n INTEGER, INTENT(OUT) :: errorflag !Return error status. -1 for error, 0 for normal REAL(8), INTENT(IN), DIMENSION(n,n) :: matrix !Input matrix REAL(8), INTENT(OUT), DIMENSION(n,n) :: inverse !Inverted matrix LOGICAL :: FLAG = .TRUE. INTEGER :: i, j, k REAL(8) :: m REAL(8), DIMENSION(n,2*n) :: augmatrix !augmented matrix REAL(8), PARAMETER :: epsilon = 1.0e-8 !Augment input matrix with an identity matrix DO i = 1, n DO j = 1, 2*n IF (j <= n ) THEN augmatrix(i,j) = matrix(i,j) ELSE IF ((i+n) == j) THEN augmatrix(i,j) = 1 Else augmatrix(i,j) = 0 ENDIF END DO END DO !Reduce augmented matrix to upper traingular form DO k =1, n-1 IF (dabs(augmatrix(k,k)) <= epsilon) THEN FLAG = .FALSE. DO i = k+1, n IF (dabs(augmatrix(i,k)) > epsilon) THEN DO j = 1,2*n augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j) END DO FLAG = .TRUE. EXIT ENDIF END DO IF (FLAG .EQV. .FALSE.) THEN PRINT*, "Matrix is non - invertible" inverse = 0 errorflag = -1 return ENDIF ENDIF DO j = k+1, n m = augmatrix(j,k)/augmatrix(k,k) DO i = k, 2*n augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i) END DO END DO END DO !Test for invertibility DO i = 1, n IF ( dabs(augmatrix(i,i)) <= epsilon ) THEN PRINT*, "Matrix is non - invertible" inverse = 0 errorflag = -1 return ENDIF END DO !Make diagonal elements as 1 DO i = 1 , n m = augmatrix(i,i) DO j = i , (2 * n) augmatrix(i,j) = (augmatrix(i,j) / m) END DO END DO !Reduced right side half of augmented matrix to identity matrix DO k = n-1, 1, -1 DO i =1, k m = augmatrix(i,k+1) DO j = k, (2*n) augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m END DO END DO END DO !store answer DO i =1, n DO j = 1, n inverse(i,j) = augmatrix(i,j+n) END DO END DO errorflag = 0 END SUBROUTINE FINDinv
Was This Post Helpful? 0
  • +
  • -

#14 mitek  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 10-February 11

Re: FORTRAN 90: Matrix Inversion

Posted 10 February 2011 - 08:03 AM

Sorry. comparison with zero should be like this dabs(augmatrix(k,k)) <= epsilon
Was This Post Helpful? 0
  • +
  • -

#15 mitek  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 0
  • Joined: 10-February 11

Re: FORTRAN 90: Matrix Inversion

Posted 10 February 2011 - 08:05 AM

And IF (FLAG .EQV. .FALSE.) THEN must be out of condition IF (dabs(augmatrix(k,k)) .gt. epsilon) THEN... ENDIF
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1