# FORTRAN 90: Matrix Inversion

Page 1 of 1

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

### #1 Louisda16th

• dream.in.assembly.code

Reputation: 15
• 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

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

Reputation: 0
• 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

### #3 masaif2009

Reputation: 0
• 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!!!

### #4 Louisda16th

• dream.in.assembly.code

Reputation: 15
• 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?

### #5 rhamdan

Reputation: 0
• 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

### #6 Louisda16th

• dream.in.assembly.code

Reputation: 15
• 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

### #7 Chris.Webster

Reputation: 0
• 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

### #8 vb10000

Reputation: 0
• 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.

### #9 vb10000

Reputation: 0
• 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.

### #10 Louisda16th

• dream.in.assembly.code

Reputation: 15
• 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.

### #11 Meister_Chicho

Reputation: 0
• 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

### #12 bfgreen

Reputation: 0
• 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.

### #13 mitek

Reputation: 0
• 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

### #14 mitek

Reputation: 0
• 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

### #15 mitek

Reputation: 0
• 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

Page 1 of 1

 .related ul { list-style-type: circle; font-size: 12px; font-weight: bold; } .related li { margin-bottom: 5px; background-position: left 7px !important; margin-left: -35px; } .related h2 { font-size: 18px; font-weight: bold; } .related a { color: blue; }