Gauss-Jordan matrix inversion: Difference between revisions
Content deleted Content added
m →version 2: changed whitespace, added a comment. |
Add Fortran example. Original code. |
||
Line 720:
-16/23 25/69 11/138 9/46
-13/23 16/69 -2/69 13/23
</pre>
=={{header|Fortran}}==
Note that the Crout algorithm is a more efficient way to invert a matrix.
<lang Fortran>!-----------------------------------------------------------------------
! gjinv - Invert a matrix, Gauss-Jordan algorithm
! A is destroyed.
!
!___Name_______Type_______________In/Out____Description_________________
! A(LDA,N) Real In An N by N matrix
! LDA Integer In Row bound of A
! N Integer In Order of matrix
! B(LDB,N) Real Out Inverse of A
! LDB Integer In Row bound of B
! IERR Integer Out 0 = no errors
! 1 = singular matrix
!-----------------------------------------------------------------------
SUBROUTINE GJINV (A, LDA, N, B, LDB, IERR)
IMPLICIT NONE
INTEGER LDA, N, LDB, IERR
REAL A(LDA,N), B(LDB,N)
REAL EPS ! machine constant
PARAMETER (EPS = 1.1920929E-07)
INTEGER I, J, K, P ! local variables
REAL F, TOL
!-----------------------------------------------------------------------
! Begin.
!-----------------------------------------------------------------------
IF (N < 1) THEN ! Validate.
IERR = -1
RETURN
ELSE IF (N > LDA .OR. N > LDB) THEN
IERR = -2
RETURN
END IF
IERR = 0
F = 0. ! Frobenius norm of A
DO J = 1, N
DO I = 1, N
F = F + A(I,J)**2
END DO
END DO
F = SQRT(F)
TOL = F * EPS
DO J = 1, N ! Set B to identity matrix.
DO I = 1, N
IF (I .EQ. J) THEN
B(I,J) = 1.
ELSE
B(I,J) = 0.
END IF
END DO
END DO
! Main loop
DO K = 1, N
F = ABS(A(K,K)) ! Find pivot.
P = K
DO I = K+1, N
IF (ABS(A(I,K)) > F) THEN
F = ABS(A(I,K))
P = I
END IF
END DO
IF (F < TOL) THEN ! Matrix is singular.
IERR = 1
RETURN
END IF
IF (P .NE. K) THEN ! Swap rows.
DO J = K, N
F = A(K,J)
A(K,J) = A(P,J)
A(P,J) = F
END DO
DO J = 1, N
F = B(K,J)
B(K,J) = B(P,J)
B(P,J) = F
END DO
END IF
F = 1. / A(K,K) ! Scale row so pivot is 1.
DO J = K, N
A(K,J) = A(K,J) * F
END DO
DO J = 1, N
B(K,J) = B(K,J) * F
END DO
DO 10 I = 1, N ! Subtract to get zeros.
IF (I .EQ. K) GO TO 10
F = A(I,K)
DO J = K, N
A(I,J) = A(I,J) - A(K,J) * F
END DO
DO J = 1, N
B(I,J) = B(I,J) - B(K,J) * F
END DO
10 CONTINUE
END DO
RETURN
END ! of gjinv</lang>
Test program:
{{works with|gfortran|8.3.0}}
<lang Fortran>! Test matrix inversion
PROGRAM TINV
IMPLICIT NONE
INTEGER I, J, K, N, IERR
PARAMETER (N = 4)
REAL AORIG(N,N), A(N,N), B(N,N), C(N,N)
DATA AORIG / -1., -4., 7., 1.,
$ -2., -1., -8., -2.,
$ 3., 6., 9., 1.,
$ 2., 2., 1., 3. /
20 FORMAT (F8.3, $)
30 FORMAT ()
DO J = 1, N
DO I = 1, N
A(I,J) = AORIG(I,J)
END DO
END DO
CALL GJINV (A, N, N, B, N, IERR)
PRINT *, 'gjinv returns #', IERR
PRINT 30
PRINT *, 'matrix:'
DO I = 1, N
DO J = 1, N
PRINT 20, AORIG(I,J)
END DO
PRINT 30
END DO
PRINT 30
PRINT *, 'inverse:'
DO I = 1, N
DO J = 1, N
PRINT 20, B(I,J)
END DO
PRINT 30
END DO
PRINT 30
DO K = 1, N
DO J = 1, N
C(J,K) = 0.
DO I = 1, N
C(J,K) = C(J,K) + AORIG(J,I) * B(I,K)
END DO
END DO
END DO
PRINT *, 'matrix @ inverse:'
DO I = 1, N
DO J = 1, N
PRINT 20, C(I,J)
END DO
PRINT 30
END DO
PRINT 30
CALL GJINV (B, N, N, A, N, IERR)
PRINT *, 'gjinv returns #', IERR
PRINT 30
PRINT *, 'inverse of inverse:'
DO I = 1, N
DO J = 1, N
PRINT 20, A(I,J)
END DO
PRINT 30
END DO
END ! of test program</lang>
{{out}}
<pre>
gjinv returns # 0
matrix:
-1.000 -2.000 3.000 2.000
-4.000 -1.000 6.000 2.000
7.000 -8.000 9.000 1.000
1.000 -2.000 1.000 3.000
inverse:
-0.913 0.246 0.094 0.413
-1.652 0.652 0.043 0.652
-0.696 0.362 0.080 0.196
-0.565 0.232 -0.029 0.565
matrix @ inverse:
1.000 0.000 -0.000 0.000
0.000 1.000 -0.000 0.000
-0.000 -0.000 1.000 0.000
-0.000 -0.000 0.000 1.000
gjinv returns # 0
inverse of inverse:
-1.000 -2.000 3.000 2.000
-4.000 -1.000 6.000 2.000
7.000 -8.000 9.000 1.000
1.000 -2.000 1.000 3.000
</pre>
|