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