(git:ec11232)
Loading...
Searching...
No Matches
gemm_square_unittest.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7PROGRAM gemm_square_unittest
8 USE kinds, ONLY: dp
9 USE mathlib, ONLY: gemm_square
10#include "../base/base_uses.f90"
11
12 IMPLICIT NONE
13
14 COMPLEX(kind=dp), DIMENSION(3, 3) :: a_in, b_in, c_in, res_c, res_c_ref
15 REAL(kind=dp), DIMENSION(3, 3) :: x_in, y_in, z_in, res_r, res_r_ref
16 REAL(kind=dp) :: tolerance = 1.0e-6_dp
17
18 ! Prepare inputs
19 a_in(1, 1) = cmplx(0.8815928086074307, 0.6726432190297216, kind=dp)
20 a_in(1, 2) = cmplx(0.7660079579530265, 0.6663301208376479, kind=dp)
21 a_in(1, 3) = cmplx(0.8910730680466552, 0.6447684662974965, kind=dp)
22 a_in(2, 1) = cmplx(0.270178070784315, 0.9380895276020503, kind=dp)
23 a_in(2, 2) = cmplx(0.4365740872106577, 0.5843460996868933, kind=dp)
24 a_in(2, 3) = cmplx(0.07466461985206008, 0.6899750234684598, kind=dp)
25 a_in(3, 1) = cmplx(0.840974290337725, 0.8395064317543346, kind=dp)
26 a_in(3, 2) = cmplx(0.5872667752635958, 0.6233467352665024, kind=dp)
27 a_in(3, 3) = cmplx(0.5024930933188588, 0.7727803824712417, kind=dp)
28
29 b_in(1, 1) = cmplx(0.08269815253296364, 0.34184561260312574, kind=dp)
30 b_in(1, 2) = cmplx(0.9876346392802493, 0.26436123295003866, kind=dp)
31 b_in(1, 3) = cmplx(0.780810836185207, 0.376036133357872, kind=dp)
32 b_in(2, 1) = cmplx(0.4787818411690774, 0.7596241356044092, kind=dp)
33 b_in(2, 2) = cmplx(0.4298758196722595, 0.4813479548810141, kind=dp)
34 b_in(2, 3) = cmplx(0.2086685419449945, 0.3860478932514133, kind=dp)
35 b_in(3, 1) = cmplx(0.34008386216308817, 0.8353095227337101, kind=dp)
36 b_in(3, 2) = cmplx(0.7379600798045334, 0.7634442366598211, kind=dp)
37 b_in(3, 3) = cmplx(0.9840849653895581, 0.9273454280026875, kind=dp)
38
39 c_in(1, 1) = cmplx(0.731192921191078, 0.9732725403607281, kind=dp)
40 c_in(1, 2) = cmplx(0.07386957805916261, 0.14228952898305391, kind=dp)
41 c_in(1, 3) = cmplx(0.12229506342104235, 0.6298697123856768, kind=dp)
42 c_in(2, 1) = cmplx(0.007352653494114958, 0.29359318766569575, kind=dp)
43 c_in(2, 2) = cmplx(0.29087841717040863, 0.48194825561460775, kind=dp)
44 c_in(2, 3) = cmplx(0.22558916232632764, 0.9229223568661166, kind=dp)
45 c_in(3, 1) = cmplx(0.5728946948517463, 0.9149335302204014, kind=dp)
46 c_in(3, 2) = cmplx(0.20475976494474424, 0.6082208447082643, kind=dp)
47 c_in(3, 3) = cmplx(0.9060121198373113, 0.008565705864987172, kind=dp)
48
49 x_in(1, 1) = 0.42929014430726375_dp
50 x_in(1, 2) = 0.21820709659663573_dp
51 x_in(1, 3) = 0.5394292090282415_dp
52 x_in(2, 1) = 0.7828031363115503_dp
53 x_in(2, 2) = 0.1422677264194132_dp
54 x_in(2, 3) = 0.25344520034350637_dp
55 x_in(3, 1) = 0.5044049742159297_dp
56 x_in(3, 2) = 0.6969177100349894_dp
57 x_in(3, 3) = 0.6999162742203425_dp
58
59 y_in(1, 1) = 0.5331333823513378_dp
60 y_in(1, 2) = 0.8001773249628732_dp
61 y_in(1, 3) = 0.2850504760853374_dp
62 y_in(2, 1) = 0.23062673571851455_dp
63 y_in(2, 2) = 0.5013417881822918_dp
64 y_in(2, 3) = 0.07530315834987644_dp
65 y_in(3, 1) = 0.2267846125008932_dp
66 y_in(3, 2) = 0.19831160340777076_dp
67 y_in(3, 3) = 0.3050258528838238_dp
68
69 z_in(1, 1) = 0.5400800562659297_dp
70 z_in(1, 2) = 0.506259700373107_dp
71 z_in(1, 3) = 0.24342576996957088_dp
72 z_in(2, 1) = 0.3517364012861689_dp
73 z_in(2, 2) = 0.04901381134580918_dp
74 z_in(2, 3) = 0.31263102401008236_dp
75 z_in(3, 1) = 0.20684120795408456_dp
76 z_in(3, 2) = 0.8051322416754273_dp
77 z_in(3, 3) = 0.5860282518273413_dp
78
79 ! Test X * Y
80
81 CALL gemm_square(x_in, 'N', y_in, 'N', res_r)
82
83 res_r_ref(1, 1) = 0.4015275411844552_dp
84 res_r_ref(1, 2) = 0.5598796466739117_dp
85 res_r_ref(1, 3) = 0.3033408981158978_dp
86 res_r_ref(2, 1) = 0.5076266966693295_dp
87 res_r_ref(2, 2) = 0.7479672000061859_dp
88 res_r_ref(2, 3) = 0.31115895421143025_dp
89 res_r_ref(3, 1) = 0.588373227540499_dp
90 res_r_ref(3, 2) = 0.8918089125227482_dp
91 res_r_ref(3, 3) = 0.4097535412069894_dp
92
93 CALL check_ref_r(res_r, res_r_ref, tolerance)
94
95 ! Test A * B
96
97 CALL gemm_square(a_in, 'N', b_in, 'N', res_c)
98
99 res_c_ref(1, 1) = cmplx(-0.5319854477298264, 2.221497049670068, kind=dp)
100 res_c_ref(1, 2) = cmplx(0.8667540454951561, 2.708638262772094, kind=dp)
101 res_c_ref(1, 3) = cmplx(0.6169940071822635, 2.7523152479106017, kind=dp)
102 res_c_ref(2, 1) = cmplx(-1.0841486927833452, 1.0783614128260843, kind=dp)
103 res_c_ref(2, 2) = cmplx(-0.5464163853751492, 2.025430919812893, kind=dp)
104 res_c_ref(2, 3) = cmplx(-0.8426527494261556, 1.872774281691093, kind=dp)
105 res_c_ref(3, 1) = cmplx(-0.884392147923063, 1.78400551956788, kind=dp)
106 res_c_ref(3, 2) = cmplx(0.3418926087394045, 2.5559945109530244, kind=dp)
107 res_c_ref(3, 3) = cmplx(0.000721037947416292, 2.5549846237243488, kind=dp)
108
109 CALL check_ref_c(res_c, res_c_ref, tolerance)
110
111 ! Test X * Y * Z
112
113 CALL gemm_square(x_in, 'N', y_in, 'N', z_in, 'N', res_r)
114
115 res_r_ref(1, 1) = 0.4765304668978436_dp
116 res_r_ref(1, 2) = 0.47494858536191625_dp
117 res_r_ref(1, 3) = 0.45054423436947794_dp
118 res_r_ref(2, 1) = 0.6016068400643494_dp
119 res_r_ref(2, 2) = 0.5441757689127916_dp
120 res_r_ref(2, 3) = 0.5397551091346775_dp
121 res_r_ref(3, 1) = 0.7162042207878401_dp
122 res_r_ref(3, 2) = 0.6714863948435401_dp
123 res_r_ref(3, 3) = 0.6621594909204267_dp
124
125 CALL check_ref_r(res_r, res_r_ref, tolerance)
126
127 ! Test A * B * C
128
129 CALL gemm_square(a_in, 'N', b_in, 'N', c_in, 'N', res_c)
130
131 res_c_ref(1, 1) = cmplx(-5.504683782712595, 3.5222601599484484, kind=dp)
132 res_c_ref(1, 2) = cmplx(-2.9563767075011937, 2.232852141215477, kind=dp)
133 res_c_ref(1, 3) = cmplx(-3.233216866606864, 3.8464986862655572, kind=dp)
134 res_c_ref(2, 1) = cmplx(-4.63714700646176, -0.11028256160215588, kind=dp)
135 res_c_ref(2, 2) = cmplx(-2.680220510420048, 0.12215466665130881, kind=dp)
136 res_c_ref(2, 3) = cmplx(-3.583889556370547, 1.091159499729193, kind=dp)
137 res_c_ref(3, 1) = cmplx(-5.468121643036268, 2.0272651357548415, kind=dp)
138 res_c_ref(3, 2) = cmplx(-3.0054301614279586, 1.4377987781538424, kind=dp)
139 res_c_ref(3, 3) = cmplx(-3.534937025955706, 2.8681214444912606, kind=dp)
140
141 CALL check_ref_c(res_c, res_c_ref, tolerance)
142
143 ! Test X^T * Y * Z
144
145 CALL gemm_square(x_in, 'T', y_in, 'N', z_in, 'N', res_r)
146
147 res_r_ref(1, 1) = 0.6462671475738445_dp
148 res_r_ref(1, 2) = 0.5760105624808499_dp
149 res_r_ref(1, 3) = 0.5852827099256332_dp
150 res_r_ref(2, 1) = 0.36007554144181786_dp
151 res_r_ref(2, 2) = 0.40420627668516457_dp
152 res_r_ref(2, 3) = 0.36217776140102986_dp
153 res_r_ref(3, 1) = 0.5978645617240842_dp
154 res_r_ref(3, 2) = 0.600788264751924_dp
155 res_r_ref(3, 3) = 0.5673424971463421_dp
156
157 CALL check_ref_r(res_r, res_r_ref, tolerance)
158
159 ! Test A^H * B * C
160
161 CALL gemm_square(a_in, 'C', b_in, 'N', c_in, 'N', res_c)
162
163 res_c_ref(1, 1) = cmplx(3.375089298965469, 5.744913993063936, kind=dp)
164 res_c_ref(1, 2) = cmplx(2.0725172551868294, 3.258926327791143, kind=dp)
165 res_c_ref(1, 3) = cmplx(3.965529787950442, 3.621340775428089, kind=dp)
166 res_c_ref(2, 1) = cmplx(2.4231309591599897, 4.665551869666368, kind=dp)
167 res_c_ref(2, 2) = cmplx(1.5937647760286848, 2.6021783330446246, kind=dp)
168 res_c_ref(2, 3) = cmplx(2.9609793918714686, 2.92153954960111, kind=dp)
169 res_c_ref(3, 1) = cmplx(3.278689562249669, 4.308656958132163, kind=dp)
170 res_c_ref(3, 2) = cmplx(2.05357432643753, 2.5060755291807237, kind=dp)
171 res_c_ref(3, 3) = cmplx(3.646272530313196, 2.5667051324585874, kind=dp)
172
173 CALL check_ref_c(res_c, res_c_ref, tolerance)
174
175CONTAINS
176! **************************************************************************************************
177!> \brief ...
178!> \param mat ...
179!> \param ref ...
180!> \param tolerance ...
181! **************************************************************************************************
182 SUBROUTINE check_ref_r(mat, ref, tolerance)
183 REAL(kind=dp), DIMENSION(3, 3) :: mat, ref
184 REAL(kind=dp) :: tolerance
185
186 INTEGER :: i, j
187
188 DO i = 1, 3
189 DO j = 1, 3
190 cpassert(abs(mat(i, j) - ref(i, j)) <= tolerance)
191 END DO
192 END DO
193 END SUBROUTINE check_ref_r
194! **************************************************************************************************
195!> \brief ...
196!> \param mat ...
197!> \param ref ...
198!> \param tolerance ...
199! **************************************************************************************************
200 SUBROUTINE check_ref_c(mat, ref, tolerance)
201 COMPLEX(kind=dp), DIMENSION(3, 3) :: mat, ref
202 REAL(kind=dp) :: tolerance
203
204 INTEGER :: i, j
205
206 DO i = 1, 3
207 DO j = 1, 3
208 cpassert(abs(mat(i, j) - ref(i, j)) <= tolerance)
209 END DO
210 END DO
211 END SUBROUTINE check_ref_c
212END PROGRAM gemm_square_unittest
gemm_square_unittest
program gemm_square_unittest
Definition gemm_square_unittest.F:7
check_ref_r
subroutine check_ref_r(mat, ref, tolerance)
...
Definition gemm_square_unittest.F:183
mathlib::gemm_square
Definition mathlib.F:89
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15