Gauss-Jordan matrix inversion: Difference between revisions
Content added Content deleted
m (→version 2: changed whitespace, added a comment.) |
(Add Fortran example. Original code.) |
||
| Line 720: | Line 720: | ||
-16/23 25/69 11/138 9/46 |
-16/23 25/69 11/138 9/46 |
||
-13/23 16/69 -2/69 13/23 |
-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> |
</pre> |
||