File: N:\mfix\model\qmomk\qmomk_tools_mod.f
1
2
3
4
5
6
7
8
9
10
11 #include "version.inc"
12
13 MODULE qmomk_tools
14
15 IMPLICIT NONE
16
17 PRIVATE
18
19 PUBLIC :: det3
20 PUBLIC :: inv3
21 PUBLIC :: diag3
22 PUBLIC :: cholesky3
23 PUBLIC :: transpose3
24 PUBLIC :: multiplyMatrix3
25
26 CONTAINS
27
28 SUBROUTINE DET3 (A, det)
29
30 IMPLICIT NONE
31
32 DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
33 DOUBLE PRECISION, INTENT(OUT) :: det
34
35 det = a(1,1)*a(2,2)*a(3,3) + a(1,2)*a(2,3)*a(3,1) + a(1,3)*a(2,1)*a(2,3) &
36 - a(3,1)*a(2,2)*a(1,3) - a(3,2)*a(2,3)*a(1,1) - a(3,3)*a(2,1)*a(1,2)
37
38 END SUBROUTINE DET3
39
40 SUBROUTINE INV3 (A, B)
41
42 IMPLICIT none
43
44 DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
45 DOUBLE PRECISION, iNTENT(OUT), DIMENSION(3,3) :: B
46
47 DOUBLE PRECISION :: detA;
48
49 CALL det3(A, detA)
50
51 IF (detA == 0.D0) THEN
52 PRINT *,'QMOMK: Null determinant in matrix'
53 ERROR_STOP
54 END IF
55
56 B(1,1) = (A(3,3)*A(2,2) - A(3,2)*A(2,3))/detA
57 B(1,2) = -(A(3,3)*A(1,2) - A(3,2)*A(1,3))/detA
58 B(1,3) = (A(2,3)*A(1,2) - A(2,2)*A(1,3))/detA
59
60 B(2,1) = -(A(3,3)*A(2,1) - A(3,1)*A(2,3))/detA
61 B(2,2) = (A(3,3)*A(1,1) - A(3,1)*A(1,3))/detA
62 B(2,3) = -(A(2,3)*A(1,1) - A(2,1)*A(1,3))/detA
63
64 B(3,1) = (A(3,2)*A(2,1) - A(3,1)*A(2,2))/detA
65 B(3,2) = -(A(3,2)*A(1,1) - A(3,1)*A(1,2))/detA
66 B(3,3) = (A(2,2)*A(1,1) - A(2,1)*A(1,2))/detA
67 END SUBROUTINE INV3
68
69 SUBROUTINE DIAG3 (a, b, c, mat)
70
71 IMPLICIT NONE
72
73 INTEGER :: i, j
74 DOUBLE PRECISION, INTENT(IN) :: a, b, c
75 DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: mat
76
77 DO i = 1, 3
78 DO j = 1, 3
79 mat(i,j) = 0.D0
80 END DO
81 END DO
82
83 mat(1,1) = a
84 mat(2,2) = b
85 mat(3,3) = c
86 END SUBROUTINE DIAG3
87
88 SUBROUTINE CHOLESKY3 (A, L)
89
90 IMPLICIT NONE
91
92 DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
93 DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: L
94
95 INTEGER :: i, j
96
97 DO i = 1, 3
98 DO j = 1, 3
99 L(i,j) = 0.D0
100 END DO
101 END DO
102
103 L(1,1) = SQRT(A(1,1))
104 IF (L(1,1) == 0.D0) THEN
105 PRINT *,'Impossible to find Cholesky decomposition.'
106 ERROR_STOP
107 ELSE
108 L(2,1) = A(2,1)/L(1,1)
109 L(2,2) = SQRT(A(2,2) - (L(2,1))**2)
110 L(3,1) = A(3,1)/L(1,1)
111
112 IF (L(2,2) == 0.D0) THEN
113 PRINT *,'Impossible to find Cholesky decomposition.'
114 ERROR_STOP
115 ELSE
116 L(3,2) = (A(3,2) - L(3,1)*L(2,1))/L(2,2)
117 L(3,3) = SQRT(A(3,3) - (L(3,1))**2 - (L(3,2))**2)
118 END IF
119 END IF
120 END SUBROUTINE CHOLESKY3
121
122 SUBROUTINE TRANSPOSE3 (A, T)
123
124 IMPLICIT NONE
125
126 DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
127 DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: T
128
129 INTEGER :: i, j
130
131 DO i = 1, 3
132 DO j = 1, 3
133 T(i,j) = A(j,i)
134 END DO
135 END DO
136
137 END SUBROUTINE TRANSPOSE3
138
139 SUBROUTINE MULTIPLYMATRIX3 (A, B, P)
140
141 IMPLICIT NONE
142
143 DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
144 DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: B
145 DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: P
146
147 INTEGER :: i, j, k
148 DOUBLE PRECISION ::sum
149
150 sum = 0.D0
151
152 DO i = 1, 3
153 DO j = 1, 3
154 DO k = 1, 3
155 sum = sum + A(i,k)*B(k,j)
156 END DO
157 P(i,j) = sum
158 sum = 0.D0
159 END DO
160 END DO
161 END SUBROUTINE MULTIPLYMATRIX3
162
163 END MODULE qmomk_tools
164