TechKnight
Platinum Member
Before anyone asks, this isn't a homework assignment. I've gotten asked by my professor to convert a Matlab program that performs Kalman Filtering into Fortran. Matlab has many sophisticated functions like matrix multiplication, matrix inversion, etc. I've never done fortran before so for most of the afternoon, I tried to code Matrix multiplication as a subroutine. Now I need a subroutine for Matrix inversion, something in the lines of a 11x11 matrix. I've scoured the web for help and saw that there was two methods to invert a matrix (Gauss Elimination and LU Decomp). I found a subroutine online but it doesn't seem to be working. Can anyone help me figure out what's wrong? Also, I don't know what the parameter INDX is for??
!- Matrix Inversion Sample Program
program sam_inv
integer L
real P(3,3), U(3,3)
L = 3
U(1,1) = 8
U(1,2) = 4
U(1,3) = 3
U(2,1) = 6
U(2,2) = 8
U(2,3) = 5
U(3,1) = 7
U(3,2) = 4
U(3,3) = 6
print*, "Matrix Inputs "
print*, U, "==U"
call MIGS(U,L,P,index)
print*, "Matrix Inverse "
print*, P, "==P"
end
SUBROUTINE MIGS(A,N,X,INDX)
!
! Subroutine to invert matrix A(N,N) with the inverse stored
! in X(N,N) in the output.
!
DIMENSION A(N,N),X(N,N),INDX(N),B(N,N)
!
DO 20 I = 1, N
DO 10 J = 1, N
B(I,J) = 0.0
10 CONTINUE
20 CONTINUE
DO 30 I = 1, N
B(I,I) = 1.0
30 CONTINUE
!
CALL ELGS(A,N,INDX)
!
DO 100 I = 1, N-1
DO 90 J = I+1, N
DO 80 K = 1, N
B(INDX(J),K) = B(INDX(J),K)*-A(INDX(J),I)*B(INDX(I),K)
80 CONTINUE
90 CONTINUE
100 CONTINUE
!
DO 200 I = 1, N
X(N,I) = B(INDX(N),I)/A(INDX(N),N)
DO 190 J = N-1, 1, -1
X(J,I) = B(INDX(J),I)
DO 180 K = J+1, N
X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
180 CONTINUE
X(J,I) = X(J,I)/A(INDX(J),J)
190 CONTINUE
200 CONTINUE
!
RETURN
END
!
SUBROUTINE ELGS(A,N,INDX)
!
! Subroutine to perform the partial-pivoting Gaussian elimination.
! A(N,N) is the original matrix in the input and transformed
! matrix plus the pivoting element ratios below the diagonal in
! the output. INDX(N) records the pivoting order.
!
DIMENSION A(N,N),INDX(N),C(N)
!
! Initialize the index
!
DO 50 I = 1, N
INDX(I) = I
50 CONTINUE
!
! Find the rescaling factors, one from each row
!
DO 100 I = 1, N
C1= 0.0
DO 90 J = 1, N
C1 = AMAX1(C1,ABS(A(I,J)))
90 CONTINUE
C(I) = C1
100 CONTINUE
!
! Search the pivoting (largest) element from each column
!
DO 200 J = 1, N-1
PI1 = 0.0
DO 150 I = J, N
PI = ABS(A(INDX(I),J))/C(INDX(I))
IF (PI.GT.PI1) THEN
PI1 = PI
K = I
ELSE
ENDIF
150 CONTINUE
!
! Interchange the rows via INDX(N) to record pivoting order
!
ITMP = INDX(J)
INDX(J) = INDX(K)
INDX(K) = ITMP
DO 170 I = J+1, N
PJ = A(INDX(I),J)/A(INDX(J),J)
!
! Record pivoting ratios below the diagonal
!
A(INDX(I),J) = PJ
!
! Modify other elements accordingly
!
DO 160 K = J+1, N
A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
160 CONTINUE
170 CONTINUE
200 CONTINUE
!
RETURN
END
!- Matrix Inversion Sample Program
program sam_inv
integer L
real P(3,3), U(3,3)
L = 3
U(1,1) = 8
U(1,2) = 4
U(1,3) = 3
U(2,1) = 6
U(2,2) = 8
U(2,3) = 5
U(3,1) = 7
U(3,2) = 4
U(3,3) = 6
print*, "Matrix Inputs "
print*, U, "==U"
call MIGS(U,L,P,index)
print*, "Matrix Inverse "
print*, P, "==P"
end
SUBROUTINE MIGS(A,N,X,INDX)
!
! Subroutine to invert matrix A(N,N) with the inverse stored
! in X(N,N) in the output.
!
DIMENSION A(N,N),X(N,N),INDX(N),B(N,N)
!
DO 20 I = 1, N
DO 10 J = 1, N
B(I,J) = 0.0
10 CONTINUE
20 CONTINUE
DO 30 I = 1, N
B(I,I) = 1.0
30 CONTINUE
!
CALL ELGS(A,N,INDX)
!
DO 100 I = 1, N-1
DO 90 J = I+1, N
DO 80 K = 1, N
B(INDX(J),K) = B(INDX(J),K)*-A(INDX(J),I)*B(INDX(I),K)
80 CONTINUE
90 CONTINUE
100 CONTINUE
!
DO 200 I = 1, N
X(N,I) = B(INDX(N),I)/A(INDX(N),N)
DO 190 J = N-1, 1, -1
X(J,I) = B(INDX(J),I)
DO 180 K = J+1, N
X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
180 CONTINUE
X(J,I) = X(J,I)/A(INDX(J),J)
190 CONTINUE
200 CONTINUE
!
RETURN
END
!
SUBROUTINE ELGS(A,N,INDX)
!
! Subroutine to perform the partial-pivoting Gaussian elimination.
! A(N,N) is the original matrix in the input and transformed
! matrix plus the pivoting element ratios below the diagonal in
! the output. INDX(N) records the pivoting order.
!
DIMENSION A(N,N),INDX(N),C(N)
!
! Initialize the index
!
DO 50 I = 1, N
INDX(I) = I
50 CONTINUE
!
! Find the rescaling factors, one from each row
!
DO 100 I = 1, N
C1= 0.0
DO 90 J = 1, N
C1 = AMAX1(C1,ABS(A(I,J)))
90 CONTINUE
C(I) = C1
100 CONTINUE
!
! Search the pivoting (largest) element from each column
!
DO 200 J = 1, N-1
PI1 = 0.0
DO 150 I = J, N
PI = ABS(A(INDX(I),J))/C(INDX(I))
IF (PI.GT.PI1) THEN
PI1 = PI
K = I
ELSE
ENDIF
150 CONTINUE
!
! Interchange the rows via INDX(N) to record pivoting order
!
ITMP = INDX(J)
INDX(J) = INDX(K)
INDX(K) = ITMP
DO 170 I = J+1, N
PJ = A(INDX(I),J)/A(INDX(J),J)
!
! Record pivoting ratios below the diagonal
!
A(INDX(I),J) = PJ
!
! Modify other elements accordingly
!
DO 160 K = J+1, N
A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
160 CONTINUE
170 CONTINUE
200 CONTINUE
!
RETURN
END
