-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmatrix-inverse.f90
83 lines (70 loc) · 2.14 KB
/
matrix-inverse.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
!About:
!This program introduces a 3x3 matrix inverse function that uses the
!array expression for cross-product. See matrix_inv below
!
!Author: Somajit Dey <[email protected]> January 2021
!Copyright (C) 2021 Somajit Dey
!License: GNU LGPL v2.1-or-later
program matrix_inversion
implicit none
real(8),dimension(3,3)::matrix,inverted
real(8)::determinant
logical::is_invertible
write(*,*)'Input matrix:'
call input(matrix)
write(*,*)
inverted=matrix_inv(matrix,determinant,is_invertible)
write(*,*)'Determinant:'
write(*,*)determinant
write(*,*)
if(.NOT.is_invertible)STOP 'Not invertible'
write(*,*)'Inverse:'
call output(inverted)
write(*,*)
write(*,*)'Product:'
call output(matmul(matrix,inverted),tolerance=1E-8)
contains
function matrix_inv(matrix,determinant,is_invertible)
real(8),dimension(3,3),intent(in)::matrix
real(8),dimension(3,3)::matrix_inv
real(8),intent(out),optional::determinant
logical,intent(out),optional::is_invertible
real(8),dimension(3,3)::tmp1,tmp2,tmp
real(8)::dotproduct
real(8),parameter::tol=2*tiny(tol) !Tolerance
tmp1=cshift(matrix,1,2)
tmp2=cshift(matrix,-1,2)
tmp=cshift((cshift(tmp1,-1,1)*tmp2-tmp1*cshift(tmp2,-1,1)),-1,1)
dotproduct=dot_product(tmp(:,1),matrix(:,1)) !Could be 2 or 3 instead of 1
if(abs(dotproduct)>tol)then
matrix_inv=transpose(tmp)/dotproduct
else
matrix_inv=0.0 !To return a null-matrix in case of non-invertibility
endif
!Optional extras follow:
if(present(determinant))determinant=dotproduct
if(present(is_invertible))then
if(abs(dotproduct)>tol)then
is_invertible=.true.
else
is_invertible=.false.
endif
endif
end function matrix_inv
subroutine input(matrix)
real(8),dimension(3,3),intent(out)::matrix
integer::row,col
do row=1,3
read(*,*)(matrix(row,col),col=1,3)
enddo
end subroutine input
subroutine output(matrix,tolerance)
real(8),dimension(3,3)::matrix
integer::row,col
real,intent(in),optional::tolerance
if(present(tolerance))where(abs(matrix)<tolerance)matrix=0.0
do row=1,3
write(*,*)(real(matrix(row,col)),col=1,3)
enddo
end subroutine output
end program matrix_inversion