Line data Source code
1 : module m_wrapper
2 : interface packmat
3 : module procedure packmat_d, packmat_z
4 : end interface
5 :
6 : interface packmatcoul
7 : module procedure packmatcoul_d, packmatcoul_z
8 : end interface
9 :
10 : interface unpackmat
11 : module procedure unpackmat_d, unpackmat_z
12 : end interface
13 :
14 : interface matvec
15 : module procedure matvec_dpd, matvec_dpz, matvec_zpd, matvec_zpz
16 : end interface
17 :
18 : interface matmat
19 : module procedure &
20 : matmat_dpdp, matmat_dpzp, matmat_zpdp, matmat_zpzp, &
21 : matmat_dmzp, matmat_dpzm, matmat_dmzm, &
22 : matmat_zmdp, matmat_zpdm, matmat_zmdm, &
23 : matmat_zmzp, matmat_zpzm, matmat_zmzm
24 : end interface
25 : contains
26 :
27 : ! --------
28 :
29 0 : function identity(n)
30 : implicit none
31 : integer, intent(in) :: n
32 : integer :: identity(n, n)
33 : integer :: i
34 0 : identity = 0
35 0 : do i = 1, n
36 0 : identity(i, i) = 1
37 : enddo
38 : end function identity
39 :
40 : ! --------
41 :
42 0 : function packmat_d(mat)
43 : use m_juDFT
44 : implicit none
45 : real, intent(in) :: mat(:, :)
46 : real :: packmat_d(size(mat, 1)*(size(mat, 1) + 1)/2)
47 : integer :: n, nn, i, j, k
48 0 : n = size(mat, 1); nn = n*(n + 1)/2
49 0 : if (size(mat, 2) /= n) call juDFT_error( 'packmat_d: array dimensions differ.')
50 0 : k = 0
51 0 : do j = 1, n
52 0 : do i = 1, j
53 0 : k = k + 1
54 0 : packmat_d(k) = mat(i, j)
55 : enddo
56 : enddo
57 : end function packmat_d
58 :
59 0 : function packmatcoul_d(mat)
60 : use m_juDFT
61 : implicit none
62 : real, intent(in) :: mat(:, :)
63 : real :: packmatcoul_d( &
64 : size(mat, 1)*(size(mat, 1) + 1)/2)
65 : integer :: n, nn, i, j, k
66 0 : n = size(mat, 1); nn = n*(n + 1)/2
67 0 : if (size(mat, 2) /= n) call juDFT_error( 'packmat_d: array dimensions differ.')
68 0 : k = 0
69 0 : do j = 1, n
70 0 : do i = 1, j
71 0 : k = k + 1
72 :
73 0 : packmatcoul_d(k) = (mat(i, j) + mat(j, i))/2.
74 :
75 : ! if(abs(mat(j,i)-mat(i,j)).gt.1e-6) then
76 : ! write(*,*) 'packmatcoul_d: input matrix not symmetric; deviation .gt. 1E-06'
77 : ! endif
78 : enddo
79 : enddo
80 : end function packmatcoul_d
81 :
82 0 : function unpackmat_d(mat)
83 : use m_juDFT
84 : implicit none
85 : real, intent(in) :: mat(:)
86 : real :: unpackmat_d( &
87 : nint(sqrt(0.25 + 2*size(mat)) - 0.5), &
88 : nint(sqrt(0.25 + 2*size(mat)) - 0.5))
89 : integer :: n, nn, i, j, k
90 0 : nn = size(mat); n = nint(sqrt(0.25 + 2*nn) - 0.5)
91 0 : k = 0
92 0 : do j = 1, n
93 0 : do i = 1, j
94 0 : k = k + 1
95 0 : unpackmat_d(i, j) = mat(k)
96 0 : unpackmat_d(j, i) = mat(k)
97 : enddo
98 : enddo
99 : end function unpackmat_d
100 :
101 0 : function packmat_z(mat)
102 : use m_juDFT
103 : implicit none
104 : complex, intent(in) :: mat(:, :)
105 : complex :: packmat_z(size(mat, 1)*(size(mat, 1) + 1)/2)
106 : integer :: n, nn, i, j, k
107 0 : n = size(mat, 1); nn = n*(n + 1)/2
108 0 : if (size(mat, 2) /= n) call juDFT_error( 'packmat_z: array dimensions differ.')
109 0 : k = 0
110 0 : do j = 1, n
111 0 : do i = 1, j
112 0 : k = k + 1
113 0 : packmat_z(k) = mat(i, j)
114 : enddo
115 : enddo
116 : end function packmat_z
117 :
118 0 : function packmatcoul_z(mat)
119 : use m_juDFT
120 : implicit none
121 : complex, intent(in) :: mat(:, :)
122 : complex :: packmatcoul_z( &
123 : size(mat, 1)*(size(mat, 1) + 1)/2)
124 : integer :: n, nn, i, j, k
125 0 : n = size(mat, 1); nn = n*(n + 1)/2
126 0 : if (size(mat, 2) /= n) call juDFT_error( 'packmat_z: array dimensions differ.')
127 0 : k = 0
128 0 : do j = 1, n
129 0 : do i = 1, j
130 0 : k = k + 1
131 0 : packmatcoul_z(k) = (mat(i, j) + conjg(mat(j, i)))/2.
132 :
133 0 : if (abs(conjg(mat(j, i)) - mat(i, j)) > 1e-4) then
134 : call juDFT_error( 'packmatcoul_z: input matrix not Hermitian; '&
135 0 : // 'deviation > 1E-04.')
136 : endif
137 : enddo
138 : enddo
139 : end function packmatcoul_z
140 :
141 0 : function unpackmat_z(mat)
142 : use m_juDFT
143 : implicit none
144 : complex, intent(in) :: mat(:)
145 : complex :: unpackmat_z( &
146 : nint(sqrt(0.25 + 2*size(mat)) - 0.5), &
147 : nint(sqrt(0.25 + 2*size(mat)) - 0.5))
148 : integer :: n, nn, i, j, k
149 0 : nn = size(mat); n = nint(sqrt(0.25 + 2*nn) - 0.5)
150 0 : k = 0
151 0 : do j = 1, n
152 0 : do i = 1, j
153 0 : k = k + 1
154 0 : unpackmat_z(i, j) = mat(k)
155 0 : unpackmat_z(j, i) = conjg(mat(k))
156 : enddo
157 : enddo
158 : end function unpackmat_z
159 : ! --------
160 :
161 0 : function matvec_dpd(mat, vec)
162 : use m_juDFT
163 : implicit none
164 : real, intent(in) :: mat(:), vec(:)
165 : real :: matvec_dpd(size(vec))
166 : integer :: nn, n
167 0 : n = size(vec)
168 0 : nn = n*(n + 1)/2
169 0 : if (size(mat) /= nn) call juDFT_error( 'matvec_dpd: input array has wrong size.')
170 0 : call dspmv('U', n, 1.0, mat, vec, 1, 0.0, matvec_dpd, 1)
171 : end function matvec_dpd
172 :
173 0 : function matvec_dpz(mat, vec)
174 : use m_juDFT
175 : implicit none
176 : real, intent(in) :: mat(:)
177 : complex, intent(in) :: vec(:)
178 : complex :: matvec_dpz(size(vec))
179 0 : real, allocatable :: vecr(:), veci(:)
180 : integer :: nn, n
181 0 : n = size(vec); allocate(vecr(n), veci(n))
182 0 : nn = n*(n + 1)/2
183 0 : if (size(mat) /= nn) call juDFT_error( 'matvec_dpz: input array has wrong size.')
184 0 : call dspmv('U', n, 1.0, mat, real(vec), 1, 0.0, vecr, 1)
185 0 : call dspmv('U', n, 1.0, mat, aimag(vec), 1, 0.0, veci, 1)
186 0 : matvec_dpz = vecr + (0.0, 1.0)*veci
187 0 : deallocate(vecr, veci)
188 : end function matvec_dpz
189 :
190 0 : function matvec_zpd(mat, vec)
191 : use m_juDFT
192 : implicit none
193 : complex, intent(in) :: mat(:)
194 : real, intent(in) :: vec(:)
195 : complex :: matvec_zpd(size(vec))
196 0 : real, allocatable :: vecr(:), veci(:)
197 : integer :: nn, n
198 0 : n = size(vec); allocate(vecr(n), veci(n))
199 0 : nn = n*(n + 1)/2
200 0 : if (size(mat) /= nn) call juDFT_error( 'matvec_zpd: input array has wrong size.')
201 0 : call dspmv('U', n, 1.0, real(mat), vec, 1, 0.0, vecr, 1)
202 0 : call dspmv('U', n, 1.0, aimag(mat), vec, 1, 0.0, veci, 1)
203 0 : matvec_zpd = vecr + (0.0, 1.0)*veci
204 0 : deallocate(vecr, veci)
205 : end function matvec_zpd
206 :
207 0 : function matvec_zpz(mat, vec)
208 : use m_juDFT
209 : implicit none
210 : complex, intent(in) :: mat(:), vec(:)
211 : complex :: matvec_zpz(size(vec))
212 : integer :: nn, n
213 0 : n = size(vec)
214 0 : nn = n*(n + 1)/2
215 0 : if (size(mat) /= nn) call juDFT_error( 'matvec_zpz: input array has wrong size.')
216 0 : call zhpmv('U', n, (1.0, 0.0), mat, vec, 1, (0.0, 0.0), matvec_zpz, 1)
217 : end function matvec_zpz
218 :
219 : ! --------
220 :
221 0 : function matmat_dpdp(mat1, mat2)
222 : use m_juDFT
223 : implicit none
224 : real, intent(in) :: mat1(:), mat2(:)
225 : real :: matmat_dpdp( &
226 : nint(sqrt(0.25 + 2*size(mat1)) - 0.5), &
227 : nint(sqrt(0.25 + 2*size(mat1)) - 0.5))
228 0 : real, allocatable :: vec(:), vec2(:)
229 : integer :: nn, n, k1, i, j, k
230 0 : nn = size(mat1)
231 0 : n = nint(sqrt(0.25 + 2*nn) - 0.5); allocate(vec(n), vec2(n))
232 0 : if (size(mat2) /= nn) &
233 0 : call juDFT_error( 'matmat_dpdp: second input array has wrong size.')
234 0 : k = 0
235 0 : do i = 1, n
236 0 : vec2(:i) = mat2(k + 1:k + i)
237 0 : k1 = k + 2*i
238 0 : do j = i + 1, n
239 0 : vec2(j) = mat2(k1)
240 0 : k1 = k1 + j
241 : enddo
242 0 : call dspmv('U', n, 1.0, mat1, vec2, 1, 0.0, vec, 1)
243 0 : matmat_dpdp(:, i) = vec
244 0 : k = k + i
245 : enddo
246 0 : deallocate(vec, vec2)
247 : end function matmat_dpdp
248 :
249 0 : function matmat_dpzp(mat1, mat2)
250 : use m_juDFT
251 : implicit none
252 : real, intent(in) :: mat1(:)
253 : complex, intent(in) :: mat2(:)
254 : complex :: matmat_dpzp( &
255 : nint(sqrt(0.25 + 2*size(mat1)) - 0.5), &
256 : nint(sqrt(0.25 + 2*size(mat1)) - 0.5))
257 0 : real, allocatable :: vecr(:), veci(:)
258 0 : complex, allocatable :: vec2(:)
259 : integer :: nn, n, k1, i, j, k
260 0 : nn = size(mat1)
261 0 : n = nint(sqrt(0.25 + 2*nn) - 0.5)
262 0 : allocate(vecr(n), veci(n), vec2(n))
263 0 : if (size(mat2) /= nn) &
264 0 : call juDFT_error( 'matmat_dpzp: second input array has wrong size.')
265 0 : k = 0
266 0 : do i = 1, n
267 0 : vec2(:i) = mat2(k + 1:k + i)
268 0 : k1 = k + 2*i
269 0 : do j = i + 1, n
270 0 : vec2(j) = conjg(mat2(k1))
271 0 : k1 = k1 + j
272 : enddo
273 0 : call dspmv('U', n, 1.0, mat1, real(vec2), 1, 0.0, vecr, 1)
274 0 : call dspmv('U', n, 1.0, mat1, aimag(vec2), 1, 0.0, veci, 1)
275 0 : matmat_dpzp(:, i) = vecr + (0.0, 1.0)*veci
276 0 : k = k + i
277 : enddo
278 0 : deallocate(vecr, veci, vec2)
279 : end function matmat_dpzp
280 :
281 0 : function matmat_zpdp(mat1, mat2)
282 : use m_juDFT
283 : implicit none
284 : complex, intent(in) :: mat1(:)
285 : real, intent(in) :: mat2(:)
286 : complex :: matmat_zpdp( &
287 : nint(sqrt(0.25 + 2*size(mat1)) - 0.5), &
288 : nint(sqrt(0.25 + 2*size(mat1)) - 0.5))
289 0 : real, allocatable :: vecr(:), veci(:)
290 0 : complex, allocatable :: vec1(:)
291 : integer :: nn, n, k1, i, j, k
292 0 : nn = size(mat1)
293 0 : n = nint(sqrt(0.25 + 2*nn) - 0.5)
294 0 : allocate(vecr(n), veci(n), vec1(n))
295 0 : if (size(mat2) /= nn) &
296 0 : call juDFT_error( 'matmat_zpdp: second input array has wrong size.')
297 0 : k = 0
298 0 : do i = 1, n
299 0 : vec1(:i) = conjg(mat1(k + 1:k + i))
300 0 : k1 = k + 2*i
301 0 : do j = i + 1, n
302 0 : vec1(j) = mat1(k1)
303 0 : k1 = k1 + j
304 : enddo
305 0 : call dspmv('U', n, 1.0, mat2, real(vec1), 1, 0.0, vecr, 1)
306 0 : call dspmv('U', n, 1.0, mat2, aimag(vec1), 1, 0.0, veci, 1)
307 0 : matmat_zpdp(i, :) = vecr + (0.0, 1.0)*veci
308 0 : k = k + i
309 : enddo
310 0 : deallocate(vecr, veci, vec1)
311 : end function matmat_zpdp
312 :
313 0 : function matmat_zpzp(mat1, mat2)
314 : use m_juDFT
315 : implicit none
316 : complex, intent(in) :: mat1(:), mat2(:)
317 : complex :: matmat_zpzp( &
318 : nint(sqrt(0.25 + 2*size(mat1)) - 0.5), &
319 : nint(sqrt(0.25 + 2*size(mat1)) - 0.5))
320 0 : complex, allocatable :: vec(:), vec2(:)
321 : integer :: nn, n, k1, i, j, k
322 0 : nn = size(mat1)
323 0 : n = nint(sqrt(0.25 + 2*nn) - 0.5); allocate(vec(n), vec2(n))
324 0 : if (size(mat2) /= nn) &
325 0 : call juDFT_error( 'matmat_zpzp: second input array has wrong size.')
326 0 : k = 0
327 0 : do i = 1, n
328 0 : vec2(:i) = mat2(k + 1:k + i)
329 0 : k1 = k + 2*i
330 0 : do j = i + 1, n
331 0 : vec2(j) = conjg(mat2(k1))
332 0 : k1 = k1 + j
333 : enddo
334 0 : call zhpmv('U', n, (1.0, 0.0), mat1, vec2, 1, (0.0, 0.0), vec, 1)
335 0 : matmat_zpzp(:, i) = vec
336 0 : k = k + i
337 : enddo
338 0 : deallocate(vec, vec2)
339 : end function matmat_zpzp
340 :
341 0 : function matmat_dpdm(mat1, mat2)
342 : use m_juDFT
343 : implicit none
344 : real, intent(in) :: mat1(:), mat2(:, :)
345 : real :: matmat_dpdm(size(mat2, 1), size(mat2, 1))
346 0 : real, allocatable :: vec(:), vec2(:)
347 : integer :: nn, n, i
348 0 : n = size(mat2, 1); nn = n*(n + 1)/2; allocate(vec(n), vec2(n))
349 0 : if (size(mat2, 2) /= n) &
350 0 : call juDFT_error( 'matmat_dpdm: dimensions of second input array differ.')
351 0 : if (size(mat1) /= nn) &
352 0 : call juDFT_error( 'matmat_dpdm: first input array has wrong size.')
353 0 : do i = 1, n
354 0 : vec2 = mat2(:, i)
355 0 : call dspmv('U', n, 1.0, mat1, vec2, 1, 0.0, vec, 1)
356 0 : matmat_dpdm(:, i) = vec
357 : enddo
358 0 : deallocate(vec, vec2)
359 : end function matmat_dpdm
360 :
361 0 : function matmat_dmdp(mat1, mat2)
362 : use m_juDFT
363 : implicit none
364 : real, intent(in) :: mat1(:, :), mat2(:)
365 : real :: matmat_dmdp(size(mat1, 1), size(mat1, 1))
366 0 : real, allocatable :: vec(:), vec2(:)
367 : integer :: nn, n, i
368 0 : n = size(mat1, 1); nn = n*(n + 1)/2; allocate(vec(n), vec2(n))
369 0 : if (size(mat1, 2) /= n) &
370 0 : call juDFT_error( 'matmat_dmdp: dimensions of first input array differ.')
371 0 : if (size(mat2) /= nn) &
372 0 : call juDFT_error( 'matmat_dmdp: second input array has wrong size.')
373 0 : do i = 1, n
374 0 : vec2 = mat1(i, :)
375 0 : call dspmv('U', n, 1.0, mat2, vec2, 1, 0.0, vec, 1)
376 0 : matmat_dmdp(i, :) = vec
377 : enddo
378 0 : deallocate(vec, vec2)
379 : end function matmat_dmdp
380 :
381 0 : function matmat_dmdm(mat1, mat2)
382 : use m_juDFT
383 : implicit none
384 : real, intent(in) :: mat1(:, :), mat2(:, :)
385 : real :: matmat_dmdm(size(mat1, 1), size(mat1, 1))
386 : integer :: n
387 0 : n = size(mat1, 1)
388 0 : if (size(mat1, 2) /= n) &
389 0 : call juDFT_error( 'matmat_dmdm: dimensions of first input array differ.')
390 0 : if (size(mat2, 1) /= n) &
391 0 : call juDFT_error( 'matmat_dmdm: second input array has wrong dimensions.')
392 0 : if (size(mat2, 2) /= n) &
393 0 : call juDFT_error( 'matmat_dmdm: dimensions of second input array differ.')
394 0 : call dgemm('N', 'N', n, n, n, 1.0, mat1, n, mat2, n, 0.0, matmat_dmdm, n)
395 : end function matmat_dmdm
396 :
397 0 : function matmat_dpzm(mat1, mat2)
398 : use m_juDFT
399 : implicit none
400 : real, intent(in) :: mat1(:)
401 : complex, intent(in) :: mat2(:, :)
402 : complex :: matmat_dpzm(size(mat2, 1), size(mat2, 1))
403 0 : real, allocatable :: vecr(:), veci(:)
404 0 : complex, allocatable :: vec2(:)
405 : integer :: nn, n, i
406 0 : n = size(mat2, 1)
407 0 : nn = n*(n + 1)/2; allocate(vecr(n), veci(n), vec2(n))
408 0 : if (size(mat2, 2) /= n) &
409 0 : call juDFT_error( 'matmat_dpzm: dimensions of second input array differ.')
410 0 : if (size(mat1) /= nn) &
411 0 : call juDFT_error( 'matmat_dpzm: first input array has wrong size.')
412 0 : do i = 1, n
413 0 : vec2 = mat2(:, i)
414 0 : call dspmv('U', n, 1.0, mat1, real(vec2), 1, 0.0, vecr, 1)
415 0 : call dspmv('U', n, 1.0, mat1, aimag(vec2), 1, 0.0, veci, 1)
416 0 : matmat_dpzm(:, i) = vecr + (0.0, 1.0)*veci
417 : enddo
418 0 : deallocate(vecr, veci, vec2)
419 : end function matmat_dpzm
420 :
421 0 : function matmat_dmzp(mat1, mat2)
422 : use m_juDFT
423 : implicit none
424 : real, intent(in) :: mat1(:, :)
425 : complex, intent(in) :: mat2(:)
426 : complex :: matmat_dmzp(size(mat1, 1), size(mat1, 1))
427 0 : complex, allocatable :: vec1(:), vec(:)
428 : integer :: nn, n, i
429 0 : n = size(mat1, 1); nn = n*(n + 1)/2; allocate(vec(n), vec1(n))
430 0 : if (size(mat1, 2) /= n) &
431 0 : call juDFT_error( 'matmat_dmzp: dimensions of first input array differ.')
432 0 : if (size(mat2) /= nn) &
433 0 : call juDFT_error( 'matmat_dmzp: second input array has wrong size.')
434 0 : do i = 1, n
435 0 : vec1 = mat1(i, :)
436 0 : call zhpmv('U', n, (1.0, 0.0), mat2, vec1, 1, (0.0, 0.0), vec, 1)
437 0 : matmat_dmzp(i, :) = conjg(vec)
438 : enddo
439 0 : deallocate(vec, vec1)
440 : end function matmat_dmzp
441 :
442 0 : function matmat_dmzm(mat1, mat2)
443 : use m_juDFT
444 : implicit none
445 : real, intent(in) :: mat1(:, :)
446 : complex, intent(in) :: mat2(:, :)
447 : complex :: matmat_dmzm(size(mat1, 1), size(mat2, 2))
448 0 : real :: matr(size(mat1, 1), size(mat2, 2)), &
449 0 : mati(size(mat1, 1), size(mat2, 2))
450 : integer :: n, n1, n2
451 0 : n1 = size(mat1, 1)
452 0 : n = size(mat1, 2)
453 0 : n2 = size(mat2, 2)
454 0 : if (size(mat2, 1) /= n) &
455 0 : call juDFT_error( 'matmat_dmzm: dimensions of matrices are inconsistent.')
456 0 : call dgemm('N', 'N', n1, n2, n, 1.0, mat1, n1, real(mat2), n, 0.0, matr, n1)
457 0 : call dgemm('N', 'N', n1, n2, n, 1.0, mat1, n1, aimag(mat2), n, 0.0, mati, n1)
458 0 : matmat_dmzm = matr + (0.0, 1.0)*mati
459 : end function matmat_dmzm
460 :
461 0 : function matmat_zpdm(mat1, mat2)
462 : use m_juDFT
463 : implicit none
464 : complex, intent(in) :: mat1(:)
465 : real, intent(in) :: mat2(:, :)
466 : complex :: matmat_zpdm(size(mat2, 1), size(mat2, 1))
467 0 : complex, allocatable :: vec(:), vec2(:)
468 : integer :: nn, n, i
469 0 : n = size(mat2, 1); nn = n*(n + 1)/2; allocate(vec(n), vec2(n))
470 0 : if (size(mat2, 2) /= n) &
471 0 : call juDFT_error( 'matmat_zpdm: dimensions of second input array differ.')
472 0 : if (size(mat1) /= nn) &
473 0 : call juDFT_error( 'matmat_zpdm: first input array has wrong size.')
474 0 : do i = 1, n
475 0 : vec2 = mat2(:, i)
476 0 : call zhpmv('U', n, (1.0, 0.0), mat1, vec2, 1, (0.0, 0.0), vec, 1)
477 0 : matmat_zpdm(:, i) = vec
478 : enddo
479 0 : deallocate(vec, vec2)
480 : end function matmat_zpdm
481 :
482 0 : function matmat_zmdp(mat1, mat2)
483 : use m_juDFT
484 : implicit none
485 : complex, intent(in) :: mat1(:, :)
486 : real, intent(in) :: mat2(:)
487 : complex :: matmat_zmdp(size(mat1, 1), size(mat1, 1))
488 0 : complex, allocatable :: vec1(:)
489 0 : real, allocatable :: vecr(:), veci(:)
490 : integer :: nn, n, i
491 0 : n = size(mat1, 1); nn = n*(n + 1)/2
492 0 : allocate(vecr(n), veci(n), vec1(n))
493 0 : if (size(mat1, 2) /= n) &
494 0 : call juDFT_error( 'matmat_zmdp: dimensions of first input array differ.')
495 0 : if (size(mat2) /= nn) &
496 0 : call juDFT_error( 'matmat_zmdp: second input array has wrong size.')
497 0 : do i = 1, n
498 0 : vec1 = conjg(mat1(i, :))
499 0 : call dspmv('U', n, 1.0, mat2, real(vec1), 1, 0.0, vecr, 1)
500 0 : call dspmv('U', n, 1.0, mat2, aimag(vec1), 1, 0.0, veci, 1)
501 0 : matmat_zmdp(i, :) = vecr - (0.0, 1.0)*veci
502 : enddo
503 0 : deallocate(vecr, veci, vec1)
504 : end function matmat_zmdp
505 :
506 0 : function matmat_zmdm(mat1, mat2)
507 : use m_juDFT
508 : implicit none
509 : complex, intent(in) :: mat1(:, :)
510 : real, intent(in) :: mat2(:, :)
511 : complex :: matmat_zmdm(size(mat1, 1), size(mat2, 2))
512 0 : real :: matr(size(mat1, 1), size(mat2, 2)), &
513 0 : mati(size(mat1, 1), size(mat2, 2))
514 : integer :: n, n1, n2
515 0 : n1 = size(mat1, 1)
516 0 : n = size(mat1, 2)
517 0 : n2 = size(mat2, 2)
518 0 : if (size(mat2, 1) /= n) &
519 0 : call juDFT_error( 'matmat_zmdm: dimensions of matrices are inconsistent.')
520 0 : call dgemm('N', 'N', n1, n2, n, 1.0, real(mat1), n1, mat2, n, 0.0, matr, n1)
521 0 : call dgemm('N', 'N', n1, n2, n, 1.0, aimag(mat1), n1, mat2, n, 0.0, mati, n1)
522 0 : matmat_zmdm = matr + (0.0, 1.0)*mati
523 : end function matmat_zmdm
524 :
525 0 : function matmat_zpzm(mat1, mat2)
526 : use m_juDFT
527 : implicit none
528 : complex, intent(in) :: mat1(:), mat2(:, :)
529 : complex :: matmat_zpzm(size(mat2, 1), size(mat2, 2))
530 0 : complex, allocatable :: vec(:), vec2(:)
531 : integer :: nn, n, i, n2
532 0 : n = size(mat2, 1); nn = n*(n + 1)/2; allocate(vec(n), vec2(n))
533 0 : n2 = size(mat2, 2)
534 0 : if (size(mat1) /= nn) &
535 0 : call juDFT_error( 'matmat_zpzm: first input array has wrong size.')
536 0 : do i = 1, n2
537 0 : vec2 = mat2(:, i)
538 0 : call zhpmv('U', n, (1.0, 0.0), mat1, vec2, 1, (0.0, 0.0), vec, 1)
539 0 : matmat_zpzm(:, i) = vec
540 : enddo
541 0 : deallocate(vec, vec2)
542 : end function matmat_zpzm
543 :
544 0 : function matmat_zmzp(mat1, mat2)
545 : use m_juDFT
546 : implicit none
547 : complex, intent(in) :: mat1(:, :), mat2(:)
548 : complex :: matmat_zmzp(size(mat1, 1), size(mat1, 1))
549 0 : complex, allocatable :: vec(:), vec2(:)
550 : integer :: nn, n, i
551 0 : n = size(mat1, 1); nn = n*(n + 1)/2; allocate(vec(n), vec2(n))
552 0 : if (size(mat1, 2) /= n) &
553 0 : call juDFT_error( 'matmat_zmzp: dimensions of first input array differ.')
554 0 : if (size(mat2) /= nn) &
555 0 : call juDFT_error( 'matmat_zmzp: second input array has wrong size.')
556 0 : do i = 1, n
557 0 : vec2 = conjg(mat1(i, :))
558 0 : call zhpmv('U', n, (1.0, 0.0), mat2, vec2, 1, (0.0, 0.0), vec, 1)
559 0 : matmat_zmzp(i, :) = conjg(vec)
560 : enddo
561 0 : deallocate(vec, vec2)
562 : end function matmat_zmzp
563 :
564 0 : function matmat_zmzm(mat1, mat2)
565 : use m_juDFT
566 : implicit none
567 : complex, intent(in) :: mat1(:, :), mat2(:, :)
568 : complex :: matmat_zmzm(size(mat1, 1), size(mat2, 2))
569 : integer :: n1, n, n2
570 : complex, parameter :: one = (1, 0), zero = 0
571 0 : n1 = size(mat1, 1)
572 0 : n = size(mat1, 2)
573 0 : n2 = size(mat2, 2)
574 0 : if (size(mat2, 1) /= n) &
575 0 : call juDFT_error( 'matmat_zmzm: dimensions of matrices are inconsistent.')
576 0 : call zgemm('N', 'N', n1, n2, n, one, mat1, n1, mat2, n, zero, matmat_zmzm, n1)
577 : end function matmat_zmzm
578 : end module m_wrapper
|
