Line data Source code
1 : !!****m* ABINIT/m_hide_blas
2 : !! NAME
3 : !! m_hide_blas
4 : !!
5 : !! FUNCTION
6 : !! This module defines interfaces for overloading BLAS routines.
7 : !! whose goal is twofold. On one hand, using generic interfaces renders
8 : !! the code more readable, especially when the routine can be compiled with
9 : !! different precision type (single-precision or double precision as done for example in the GW code)
10 : !! On the other hand, the generic interfaces defined here introduce a programming
11 : !! layer that can be exploited for interfacing non-standard libraries such as for
12 : !! example CUBLAS routines for GPU computations.
13 : !!
14 : !! COPYRIGHT
15 : !! Copyright (C) 1992-2026 ABINIT group (MG)
16 : !! This file is distributed under the terms of the
17 : !! GNU General Public License, see ~abinit/COPYING
18 : !! or http://www.gnu.org/copyleft/gpl.txt .
19 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
20 : !!
21 : !! NOTES
22 : !!
23 : !! The convention about names of interfaced routine is: x<name>,
24 : !! where <name> is usually equal to the name of the standard routine
25 : !! without the first character specifying the type and kind.
26 : !! The full list of names is reported below.
27 : !! BLAS procedures interfaced in this module are marked with an asterisk.
28 : !! A complete list of possible overloaded interfaces is provided as guide for future additions.
29 : !!
30 : !! ================
31 : !! ==== BLAS 1 ====
32 : !! ================
33 : !! FUNCTION idamax isamax icamax izamax ---> XIAMAX(n,dx,incx)
34 : !! * FUNCTION snrm2 dnrm2 scnrm2 dznmr2 ---> XNRM2(n,x,incx)
35 : !! FUNCTION sasum dasum scasum dzasum ---> XASUM(n,x,incx)
36 : !! * FUNCTION cdotu zdotu ---> XDOTU(n,x,incx,y,incy)
37 : !! * FUNCTION cdotc zdotc ---> XDOTC(n,x,incx,y,incy)
38 : !! FUNCTION sdot ddot ---> XDOT(n,x,incx,y,incy)
39 : !! FUNCTION sdsdot sdot ---> XDSDOT(n,x,incx,y,incy)
40 : !! SUBROUTINE saxpy daxpy caxpy zaxpy ---> XAXPY(n,ca,cx,incx,cy,incy)
41 : !! * SUBROUTINE scopy dcopy ccopy zcopy ---> XCOPY(n,cx,incx,cy,incy)
42 : !! SUBROUTINE srotg drotg crotg zrotg ---> XROTG(a,b,c,s)
43 : !! SUBROUTINE srot drot csrot zdrot ---> XROT(n,x,incx,y,incy,c,s)
44 : !! * SUBROUTINE sscal dscal cscal zscal
45 : !! csscal zdscal ---> XSCAL(n,a,x,incx)
46 : !! SUBROUTINE sswap dswap cswap zswap ---> XSWAP(n,x,incx,y,incy)
47 : !!
48 : !! ================
49 : !! ==== BLAS 2 ====
50 : !! ================
51 : !! SUBROUTINE sgbmv dgbmv cgbmv zgbmv ---> XGBMV(trans,m,kl,ku,n,alpha,A,lda,X,incx,beta,Y,incy)
52 : !! * SUBROUTINE sgemv dgemv cgemv zgemv ---> XGEMV(trans,m,n,alpha,A,lda,X,incx,beta,Y,incy)
53 : !! * SUBROUTINE cgerc zgerc ---> XGERC(m,n,alpha,x,incx,y,incy,A,lda)
54 : !! SUBROUTINE cgeru zgeru ---> XGERU(m,n,alpha,x,incx,y,incy,A,lda)
55 : !! SUBROUTINE chbmv zhbmv ---> XHBMV(uplo,n,k,alpha,A,lda,X,incx,beta,Y,incy)
56 : !! SUBROUTINE chemv zhemv ---> XHEMV(uplo,n,alpha,A,lda,X,incx,beta,Y,incy)
57 : !! * SUBROUTINE cher zher ---> XHER(uplo,n,alpha,x,incx,A,lda)
58 : !! SUBROUTINE cher2 zher2 ---> XHER2(uplo,n,alpha,x,incx,y,incy,A,lda)
59 : !! SUBROUTINE chpr zhpr ---> XHPR(uplo,n,alpha,x,incx,AP)
60 : !! SUBROUTINE chpr2 zhpr2 ---> XHPR2(uplo,n,alpha,x,incx,y,incy,AP)
61 : !! SUBROUTINE chpmv zhpmv ---> XHPMV(uplo,n,alpha,AP,X,incx,beta,Y,incy)
62 : !! SUBROUTINE stbmv dtbmv ctbmv ztbmv ---> XTBMV(uplo,trans,diag,n,k,A,lda,X,incx)
63 : !! SUBROUTINE stpmv dtpmv ctpmv ztpmv ---> XTPMV(uplo,trans,diag,n,AP,X,incx)
64 : !! SUBROUTINE strmv dtrmv ctrmv ztrmv ---> XTRMV(uplo,trans,diag,n,A,lda,X,incx)
65 : !! SUBROUTINE ssymv dsymv ---> XSYMV(uplo,n,alpha,A,lda,X,incx,beta,Y,incy)
66 : !! SUBROUTINE ssbmv dsbmv ---> XSBMV(uplo,n,k,alpha,A,lda,X,incx,beta,Y,incy)
67 : !! SUBROUTINE sspmv dspmv ---> XSPMV(uplo,n,alpha,AP,X,incx,beta,Y,incy)
68 : !! SUBROUTINE stbsv dtbsv ctbsv ztbsv ---> XTBSV(uplo,trans,diag,n,k,A,lda,X,incx)
69 : !! SUBROUTINE stpsv dtpsv ctpsv ztpsv ---> XTPSV(uplo,trans,diag,n,AP,X,incx)
70 : !! SUBROUTINE strsv dtrsv ctrsv ztrsv ---> XTRSV(uplo,trans,diag,n,A,lda,X,incx)
71 : !! SUBROUTINE sger dger ---> XGER(m,n,alpha,x,incx,y,incy,A,lda)
72 : !! SUBROUTINE sspr dspr ---> XSPR(uplo,n,alpha,x,incx,AP)
73 : !! SUBROUTINE sspr2 dspr2 ---> XSPR2(uplo,n,alpha,x,incx,y,incy,AP)
74 : !! SUBROUTINE ssyr dsyr ---> XSYR(uplo,n,alpha,x,incx,A,lda)
75 : !! SUBROUTINE ssyr2 dsyr2 ---> XSYR2(uplo,n,alpha,x,incx,y,incy,A,lda)
76 : !!
77 : !! ================
78 : !! ==== BLAS 3 ====
79 : !! ================
80 : !! * SUBROUTINE sgemm dgemm cgemm zgemm ---> XGEMM(transA,transB,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc)
81 : !! SUBROUTINE chemm zhemm ---> XHEMM(side,uplo,m,n,alpha,A,lda,B,ldb,beta,C,ldc)
82 : !! SUBROUTINE cher2k zher2k ---> XHER2K(uplo,trans,n,k,alpha,A,lda,B,ldb,beta,C,ldc)
83 : !! * SUBROUTINE cherk zherk ---> XHERK(uplo,trans,n,k,alpha,A,lda,beta,C,ldc)
84 : !! SUBROUTINE ssymm dsymm csymm zsymm ---> XSYMM(side,uplo,m,n,alpha,A,lda,B,ldb,beta,C,ldc)
85 : !! SUBROUTINE ssyr2k dsyr2k csyr2k zsyr2k ---> XSYR2K(uplo,trans,n,k,alpha,A,lda,B,ldb,beta,C,ldc)
86 : !! SUBROUTINE ssyrk dsyrk csyrk zsyrk ---> XSYRK(uplo,trans,n,k,alpha,A,lda,beta,C,ldc)
87 : !! SUBROUTINE strmm dtrmm ctrmm ztrmm ---> XTRMM(side,uplo,transa,diag,m,n,alpha,A,lda,B,ldb)
88 : !! SUBROUTINE strsm dtrsm ctrsm ztrsm ---> XTRSM(side,uplo,transa,diag,m,n,alpha,A,lda,B,ldb)
89 : !!-------------------------------------------------------------------------------
90 : !!
91 : !! SOURCE
92 :
93 : #if defined HAVE_CONFIG_H
94 : #include "config.h"
95 : #endif
96 :
97 : #include "abi_common.h"
98 :
99 : module m_hide_blas
100 :
101 : use defs_basis
102 : use m_abicore
103 : use m_errors
104 :
105 : implicit none
106 :
107 : private
108 :
109 : !BLAS1
110 : public :: xnrm2
111 : public :: xscal
112 : public :: xdotu
113 : public :: xdotc
114 : public :: xcopy
115 :
116 : !BLAS2
117 : public :: xgemv
118 : public :: xgerc
119 : public :: xher
120 :
121 : !BLAS3
122 : public :: xgemm
123 : public :: xherk
124 :
125 : ! Helper functions
126 : public :: blas_cholesky_ortho ! Cholesky orthogonalization.
127 :
128 : public :: sqmat_itranspose ! In-place transposition of a square matrix.
129 : public :: sqmat_otranspose ! out-of-place transposition of a square matrix.
130 :
131 : public :: sqmat_iconjgtrans ! in-place conjugate transpose of a square matrix.
132 : public :: sqmat_oconjgtrans ! out-of-place conjugate transpose of a square matrix.
133 :
134 : !----------------------------------------------------------------------
135 :
136 : interface xnrm2
137 : !
138 : function snrm2 ( n, x, incx )
139 : use defs_basis
140 : real(sp) :: snrm2
141 : integer,intent(in) :: incx, n
142 : real(sp),intent(in) :: x( * )
143 : end function snrm2
144 : !
145 : function dnrm2 ( n, x, incx )
146 : use defs_basis
147 : real(dp) :: dnrm2
148 : integer,intent(in) :: incx, n
149 : real(dp),intent(in) :: x( * )
150 : end function dnrm2
151 : !
152 : function scnrm2( n, x, incx )
153 : use defs_basis
154 : real(sp) :: scnrm2
155 : integer,intent(in) :: incx, n
156 : complex(sp),intent(in) :: x( * )
157 : end function scnrm2
158 : !
159 : function dznrm2( n, x, incx )
160 : use defs_basis
161 : real(dp) :: dznrm2
162 : integer,intent(in) :: incx, n
163 : complex(dp),intent(in) :: x( * )
164 : end function dznrm2
165 : !
166 : end interface xnrm2
167 :
168 : !-------------------------------------------------------------------------------
169 :
170 : interface xscal
171 : !
172 : subroutine sscal(n,sa,sx,incx)
173 : use defs_basis
174 : integer :: incx
175 : integer :: n
176 : real(sp) :: sa
177 : real(sp) :: sx(*)
178 : end subroutine sscal
179 : !
180 : subroutine dscal(n,da,dx,incx)
181 : use defs_basis
182 : integer :: incx
183 : integer :: n
184 : real(dp):: da
185 : real(dp):: dx(*)
186 : end subroutine dscal
187 : !
188 : subroutine cscal(n,ca,cx,incx)
189 : use defs_basis
190 : integer :: incx
191 : integer :: n
192 : complex(sp) :: ca
193 : complex(sp) :: cx(*)
194 : end subroutine cscal
195 : !
196 : subroutine zscal(n,za,zx,incx)
197 : use defs_basis
198 : integer :: incx
199 : integer :: n
200 : complex(dp) :: za
201 : complex(dp) :: zx(*)
202 : end subroutine zscal
203 : !
204 : subroutine csscal(n,sa,cx,incx)
205 : use defs_basis
206 : integer :: incx
207 : integer :: n
208 : real(sp) :: sa
209 : complex(sp) :: cx(*)
210 : end subroutine csscal
211 : !
212 : subroutine zdscal(n,da,zx,incx)
213 : use defs_basis
214 : integer :: incx
215 : integer :: n
216 : real(dp) :: da
217 : complex(dp) :: zx(*)
218 : end subroutine zdscal
219 : !
220 : end interface xscal
221 :
222 : !-------------------------------------------------------------------------------
223 :
224 : interface xdotu
225 : !
226 : #ifdef HAVE_LINALG_ZDOTU_BUG
227 : module procedure cdotu
228 : module procedure zdotu
229 : #else
230 : function cdotu(n,cx,incx,cy,incy)
231 : use defs_basis
232 : complex(sp) :: cdotu
233 : complex(sp),intent(in) :: cx(*),cy(*)
234 : integer,intent(in) :: incx,incy,n
235 : end function cdotu
236 : !
237 : function zdotu(n,zx,incx,zy,incy)
238 : use defs_basis
239 : complex(dp) :: zdotu
240 : complex(dp),intent(in) :: zx(*),zy(*)
241 : integer,intent(in) :: incx,incy,n
242 : end function zdotu
243 : #endif
244 : !
245 : end interface xdotu
246 :
247 : !-------------------------------------------------------------------------------
248 :
249 :
250 : ! CDOTC, CDOTU, ZDOTC, and ZDOTU are problematic if Mac OS X's Vec lib is used.
251 : ! See http://developer.apple.com/hardwaredrivers/ve/errata.html.
252 : ! If needed, we replace them with plain Fortran code.
253 :
254 : interface xdotc
255 : !
256 : #ifdef HAVE_LINALG_ZDOTC_BUG
257 : module procedure cdotc
258 : module procedure zdotc
259 : #else
260 : function cdotc(n,cx,incx,cy,incy)
261 : use defs_basis
262 : complex(sp) :: cdotc
263 : complex(sp),intent(in) :: cx(*),cy(*)
264 : integer,intent(in) :: incx,incy,n
265 : end function cdotc
266 : !
267 : function zdotc(n,zx,incx,zy,incy)
268 : use defs_basis
269 : complex(dp) :: zdotc
270 : complex(dp),intent(in) :: zx(*),zy(*)
271 : integer,intent(in) :: incx,incy,n
272 : end function zdotc
273 : #endif
274 : !
275 : end interface xdotc
276 :
277 : !-------------------------------------------------------------------------------
278 :
279 : interface xcopy
280 : !module procedure ABI_xcopy
281 : !
282 : subroutine scopy(n,sx,incx,sy,incy)
283 : use defs_basis
284 : integer,intent(in) :: incx
285 : integer,intent(in) :: incy
286 : integer,intent(in) :: n
287 : real(sp),intent(in) :: sx(*)
288 : real(sp),intent(inout) :: sy(*)
289 : end subroutine scopy
290 : !
291 : subroutine dcopy(n,dx,incx,dy,incy)
292 : use defs_basis
293 : integer,intent(in) :: incx
294 : integer,intent(in) :: incy
295 : integer,intent(in) :: n
296 : real(dp),intent(in) :: dx(*)
297 : real(dp),intent(inout) :: dy(*)
298 : end subroutine dcopy
299 : !
300 : subroutine ccopy(n,cx,incx,cy,incy)
301 : use defs_basis
302 : integer,intent(in) :: incx
303 : integer,intent(in) :: incy
304 : integer,intent(in) :: n
305 : complex(sp),intent(in) :: cx(*)
306 : complex(sp),intent(inout) :: cy(*)
307 : end subroutine ccopy
308 : !
309 : subroutine zcopy(n,cx,incx,cy,incy)
310 : use defs_basis
311 : integer,intent(in) :: incx
312 : integer,intent(in) :: incy
313 : integer,intent(in) :: n
314 : complex(dp),intent(in) :: cx(*)
315 : complex(dp),intent(inout) :: cy(*)
316 : end subroutine zcopy
317 : !
318 : end interface xcopy
319 :
320 : !-------------------------------------------------------------------------------
321 :
322 : interface xgemv
323 : !
324 : subroutine sgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
325 : use defs_basis
326 : real(sp),intent(in) :: alpha, beta
327 : integer,intent(in) :: incx, incy, lda, m, n
328 : character(len=1),intent(in) :: trans
329 : real(sp),intent(in) :: a( lda, * ), x( * )
330 : real(sp),intent(inout) :: y( * )
331 : end subroutine sgemv
332 : !
333 : subroutine dgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
334 : use defs_basis
335 : real(dp),intent(in) :: alpha, beta
336 : integer,intent(in) :: incx, incy, lda, m, n
337 : character(len=1),intent(in) :: trans
338 : real(dp),intent(in) :: a( lda, * ), x( * )
339 : real(dp),intent(inout) :: y( * )
340 : end subroutine dgemv
341 : !
342 : subroutine cgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
343 : use defs_basis
344 : character(len=1),intent(in) :: trans
345 : integer,intent(in) :: incx, incy, lda, m, n
346 : complex(sp),intent(in) :: alpha, beta
347 : complex(sp),intent(in) :: a( lda, * ), x( * )
348 : complex(sp),intent(inout) :: y( * )
349 : end subroutine cgemv
350 : !
351 : subroutine zgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
352 : use defs_basis
353 : character(len=1),intent(in) :: trans
354 : integer,intent(in) :: incx, incy, lda, m, n
355 : complex(dp),intent(in) :: alpha, beta
356 : complex(dp),intent(in) :: a( lda, * ), x( * )
357 : complex(dp),intent(inout) :: y( * )
358 : end subroutine zgemv
359 : !
360 : end interface xgemv
361 :
362 : !-------------------------------------------------------------------------------
363 :
364 : interface xgerc
365 : !
366 : subroutine cgerc ( m, n, alpha, x, incx, y, incy, a, lda )
367 : use defs_basis
368 : complex(sp),intent(in) :: alpha
369 : integer,intent(in) :: incx, incy, lda, m, n
370 : complex(sp),intent(inout) :: a( lda, * )
371 : complex(sp),intent(in) :: x( * ), y( * )
372 : end subroutine cgerc
373 : !
374 : subroutine zgerc ( m, n, alpha, x, incx, y, incy, a, lda )
375 : use defs_basis
376 : complex(dp),intent(in) :: alpha
377 : integer,intent(in) :: incx, incy, lda, m, n
378 : complex(dp),intent(inout) :: a( lda, * )
379 : complex(dp),intent(in) :: x( * ), y( * )
380 : end subroutine zgerc
381 : !
382 : end interface xgerc
383 :
384 : !-------------------------------------------------------------------------------
385 :
386 : interface xher
387 : !
388 : subroutine cher ( uplo, n, alpha, x, incx, a, lda )
389 : use defs_basis
390 : character(len=1),intent(in) :: uplo
391 : real(sp),intent(in) :: alpha
392 : integer,intent(in) :: incx, lda, n
393 : complex(sp),intent(inout) :: a( lda, * )
394 : complex(sp),intent(in) :: x( * )
395 : end subroutine cher
396 : !
397 : subroutine zher ( uplo, n, alpha, x, incx, a, lda )
398 : use defs_basis
399 : character(len=1),intent(in) :: uplo
400 : real(dp),intent(in) :: alpha
401 : integer,intent(in) :: incx, lda, n
402 : complex(dp),intent(inout) :: a( lda, * )
403 : complex(dp),intent(in) :: x( * )
404 : end subroutine zher
405 : !
406 : end interface xher
407 :
408 : !-------------------------------------------------------------------------------
409 :
410 : interface xgemm
411 : !
412 : subroutine sgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
413 : use defs_basis
414 : character(len=1),intent(in) :: transa, transb
415 : integer,intent(in) :: m, n, k, lda, ldb, ldc
416 : real(sp),intent(in) :: alpha, beta
417 : real(sp),intent(in) :: a( lda, * ), b( ldb, * )
418 : real(sp),intent(inout) :: c( ldc, * )
419 : end subroutine sgemm
420 : !
421 : subroutine dgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
422 : use defs_basis
423 : character(len=1),intent(in) :: transa, transb
424 : integer,intent(in) :: m, n, k, lda, ldb, ldc
425 : real(dp),intent(in) :: alpha, beta
426 : real(dp),intent(in) :: a( lda, * ), b( ldb, * )
427 : real(dp),intent(inout) :: c( ldc, * )
428 : end subroutine dgemm
429 : !
430 : subroutine cgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
431 : use defs_basis
432 : character(len=1),intent(in) :: transa, transb
433 : integer,intent(in) :: m, n, k, lda, ldb, ldc
434 : complex(sp),intent(in) :: alpha, beta
435 : complex(sp),intent(in) :: a( lda, * ), b( ldb, * )
436 : complex(sp),intent(inout) :: c( ldc, * )
437 : end subroutine cgemm
438 : !
439 : subroutine zgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
440 : use defs_basis
441 : character(len=1),intent(in) :: transa, transb
442 : integer,intent(in) :: m, n, k, lda, ldb, ldc
443 : complex(dp),intent(in) :: alpha, beta
444 : complex(dp),intent(in) :: a( lda, * ), b( ldb, * )
445 : complex(dp),intent(inout) :: c( ldc, * )
446 : end subroutine zgemm
447 : !
448 : end interface xgemm
449 :
450 : interface xherk
451 : !
452 : subroutine cherk( uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
453 : use defs_basis
454 : character(len=1),intent(in) :: uplo
455 : character(len=1),intent(in) :: trans
456 : integer,intent(in) :: n,k,lda,ldc
457 : real(sp),intent(in) :: alpha
458 : complex(sp),intent(in) :: a( lda, * )
459 : real(sp),intent(in) :: beta
460 : complex(sp),intent(inout) :: c( ldc, * )
461 : end subroutine cherk
462 : !
463 : subroutine zherk( uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
464 : use defs_basis
465 : character(len=1),intent(in) :: uplo
466 : character(len=1),intent(in) :: trans
467 : integer,intent(in) :: n,k,lda,ldc
468 : real(dp),intent(in) :: alpha
469 : complex(dp),intent(in) :: a( lda, * )
470 : real(dp),intent(in) :: beta
471 : complex(dp),intent(inout) :: c( ldc, * )
472 : end subroutine zherk
473 : !
474 : end interface xherk
475 :
476 : !-------------------------------------------------------------------------------
477 :
478 : interface blas_cholesky_ortho
479 : module procedure blas_cholesky_ortho_spc
480 : module procedure blas_cholesky_ortho_dpc
481 : end interface blas_cholesky_ortho
482 :
483 : interface sqmat_itranspose
484 : module procedure sqmat_itranspose_sp
485 : module procedure sqmat_itranspose_dp
486 : module procedure sqmat_itranspose_spc
487 : module procedure sqmat_itranspose_dpc
488 : end interface sqmat_itranspose
489 :
490 : interface sqmat_otranspose
491 : module procedure sqmat_otranspose_sp
492 : module procedure sqmat_otranspose_dp
493 : module procedure sqmat_otranspose_spc
494 : module procedure sqmat_otranspose_dpc
495 : end interface sqmat_otranspose
496 :
497 : interface sqmat_iconjgtrans
498 : module procedure sqmat_iconjgtrans_spc
499 : module procedure sqmat_iconjgtrans_dpc
500 : end interface sqmat_iconjgtrans
501 :
502 : interface sqmat_oconjgtrans
503 : module procedure sqmat_oconjgtrans_spc
504 : module procedure sqmat_oconjgtrans_dpc
505 : end interface sqmat_oconjgtrans
506 :
507 : real(dp),private,parameter :: zero_dp = 0._dp
508 : real(dp),private,parameter :: one_dp = 1._dp
509 :
510 : complex(dp),private,parameter :: czero_dpc = (0._dp,0._dp)
511 : complex(dp),private,parameter :: cone_dpc = (1._dp,0._dp)
512 :
513 : CONTAINS !========================================================================================
514 :
515 : ! CDOTC, CDOTU, ZDOTC, and ZDOTU are problematic if Mac OS X's Vec lib is used.
516 : ! See http://developer.apple.com/hardwaredrivers/ve/errata.html.
517 : ! Here we replace them with plain Fortran code.
518 :
519 : #ifdef HAVE_LINALG_ZDOTC_BUG
520 : !#warning "Using internal replacement for zdotc. External library cannot be used"
521 : #include "replacements/cdotc.f"
522 : #include "replacements/zdotc.f"
523 : #endif
524 :
525 : #ifdef HAVE_LINALG_ZDOTU_BUG
526 : !#warning "Using internal replacement for zdotu. External library cannot be used"
527 : #include "replacements/cdotu.f"
528 : #include "replacements/zdotu.f"
529 : #endif
530 :
531 : !----------------------------------------------------------------------
532 :
533 : !!***
534 :
535 : !!****f* m_hide_blas/blas_cholesky_ortho_spc
536 : !! NAME
537 : !! blas_cholesky_ortho_spc
538 : !!
539 : !! FUNCTION
540 : !! Performs the Cholesky orthonormalization of the vectors stored in iomat.
541 : !!
542 : !! INPUTS
543 : !! vec_size=Size of each vector.
544 : !! nvec=Number of vectors in iomat
545 : !!
546 : !! OUTPUT
547 : !! cf_ovlp=Cholesky factorization of the overlap matrix. ovlp = U^H U with U upper triangle matrix returned in cf_ovlp
548 : !!
549 : !! SIDE EFFECTS
550 : !! iomat(vec_size,nvec)
551 : !! input: Input set of vectors.
552 : !! output: Orthonormalized set.
553 : !!
554 : !! SOURCE
555 :
556 0 : subroutine blas_cholesky_ortho_spc(vec_size,nvec,iomat,cf_ovlp,use_gemm)
557 :
558 : !Arguments ------------------------------------
559 : integer,intent(in) :: vec_size,nvec
560 : logical,optional,intent(in) :: use_gemm
561 : complex(sp),intent(inout) :: iomat(vec_size,nvec)
562 : complex(sp),intent(out) :: cf_ovlp(nvec,nvec)
563 :
564 : !Local variables ------------------------------
565 : !scalars
566 : integer :: ierr
567 : logical :: my_usegemm
568 : character(len=500) :: msg
569 : ! *************************************************************************
570 :
571 : ! 1) Calculate overlap_ij = <phi_i|phi_j>
572 : ! TODO: use dsyrk
573 0 : my_usegemm = .FALSE.; if (PRESENT(use_gemm)) my_usegemm = use_gemm
574 :
575 0 : if (my_usegemm) then
576 0 : call xgemm("Conjugate","Normal",nvec,nvec,vec_size,cone_sp,iomat,vec_size,iomat,vec_size,czero_sp,cf_ovlp,nvec)
577 : else
578 0 : call xherk("U","C", nvec, vec_size, one_sp, iomat, vec_size, zero_sp, cf_ovlp, nvec)
579 : end if
580 : !
581 : ! 2) Cholesky factorization: ovlp = U^H U with U upper triangle matrix.
582 0 : call CPOTRF('U',nvec,cf_ovlp,nvec,ierr)
583 0 : if (ierr/=0) then
584 0 : write(msg,'(a,i0)')' ZPOTRF returned info= ',ierr
585 0 : ABI_ERROR(msg)
586 : end if
587 : !
588 : ! 3) Solve X U = io_mat. On exit iomat is orthonormalized.
589 0 : call CTRSM('Right','Upper','Normal','Normal',vec_size,nvec,cone_sp,cf_ovlp,nvec,iomat,vec_size)
590 :
591 0 : end subroutine blas_cholesky_ortho_spc
592 : !!***
593 :
594 : !----------------------------------------------------------------------
595 :
596 : !!****f* m_hide_blas/blas_cholesky_ortho_dpc
597 : !! NAME
598 : !! blas_cholesky_ortho_dpc
599 : !!
600 : !! FUNCTION
601 : !! Performs the Cholesky orthonormalization of the vectors stored in iomat.
602 : !!
603 : !! INPUTS
604 : !! vec_size=Size of each vector.
605 : !! nvec=Number of vectors in iomat
606 : !!
607 : !! OUTPUT
608 : !! cf_ovlp=Cholesky factorization of the overlap matrix. ovlp = U^H U with U upper triangle matrix returned in cf_ovlp
609 : !!
610 : !! SIDE EFFECTS
611 : !! iomat(vec_size,nvec)
612 : !! input: Input set of vectors.
613 : !! output: Orthonormalized set.
614 : !!
615 : !! SOURCE
616 :
617 0 : subroutine blas_cholesky_ortho_dpc(vec_size,nvec,iomat,cf_ovlp,use_gemm)
618 :
619 : !Arguments ------------------------------------
620 : integer,intent(in) :: vec_size,nvec
621 : logical,optional,intent(in) :: use_gemm
622 : complex(dp),intent(inout) :: iomat(vec_size,nvec)
623 : complex(dp),intent(out) :: cf_ovlp(nvec,nvec)
624 :
625 : !Local variables ------------------------------
626 : !scalars
627 : integer :: ierr
628 : logical :: my_usegemm
629 : character(len=500) :: msg
630 : ! *************************************************************************
631 :
632 : ! 1) Calculate overlap_ij = <phi_i|phi_j>
633 0 : my_usegemm = .FALSE.; if (PRESENT(use_gemm)) my_usegemm = use_gemm
634 :
635 0 : if (my_usegemm) then
636 0 : call xgemm("Conjugate","Normal",nvec,nvec,vec_size,cone_dpc,iomat,vec_size,iomat,vec_size,czero_dpc,cf_ovlp,nvec)
637 : else
638 0 : call xherk("U","C", nvec, vec_size, one_dp, iomat, vec_size, zero_dp, cf_ovlp, nvec)
639 : end if
640 : !
641 : ! 2) Cholesky factorization: ovlp = U^H U with U upper triangle matrix.
642 0 : call ZPOTRF('U',nvec,cf_ovlp,nvec,ierr)
643 0 : if (ierr/=0) then
644 0 : write(msg,'(a,i0)')' ZPOTRF returned info= ',ierr
645 0 : ABI_ERROR(msg)
646 : end if
647 : !
648 : ! 3) Solve X U = io_mat. On exit io_mat is orthonormalized.
649 0 : call ZTRSM('Right','Upper','Normal','Normal',vec_size,nvec,cone_dpc,cf_ovlp,nvec,iomat,vec_size)
650 :
651 0 : end subroutine blas_cholesky_ortho_dpc
652 : !!***
653 :
654 : !----------------------------------------------------------------------
655 :
656 : !!****f* m_hide_blas/sqmat_itranspose_sp
657 : !! NAME
658 : !! sqmat_itranspose_sp
659 : !!
660 : !! FUNCTION
661 : !! Compute alpha * mat^T in place. target: single precision real matrix.
662 : !!
663 : !! INPUTS
664 : !! n=size of the matrix
665 : !! [alpha]=scalar, set to 1.0 if not present
666 : !!
667 : !! SIDE EFFECTS
668 : !! mat(n,n)=in output, it contains alpha * mat^T.
669 : !!
670 : !! SOURCE
671 :
672 0 : subroutine sqmat_itranspose_sp(n,mat,alpha)
673 :
674 : !Arguments ------------------------------------
675 : !scalars
676 : integer,intent(in) :: n
677 : real(sp),optional,intent(in) :: alpha
678 : !arrays
679 : real(sp),intent(inout) :: mat(n,n)
680 : ! *************************************************************************
681 :
682 : #ifdef HAVE_LINALG_MKL_IMATCOPY
683 0 : if (PRESENT(alpha)) then
684 0 : call mkl_simatcopy("Column", "Trans", n, n, alpha, mat, n, n)
685 : else
686 0 : call mkl_simatcopy("Column", "Trans", n, n, one_sp, mat, n, n)
687 : end if
688 : #else
689 : ! Fallback to Fortran.
690 : if (PRESENT(alpha)) then
691 : mat = alpha * TRANSPOSE(mat)
692 : else
693 : mat = TRANSPOSE(mat)
694 : end if
695 : #endif
696 :
697 0 : end subroutine sqmat_itranspose_sp
698 : !!***
699 :
700 : !----------------------------------------------------------------------
701 :
702 : !!****f* m_hide_blas/sqmat_itranspose_dp
703 : !! NAME
704 : !! sqmat_itranspose_dp
705 : !!
706 : !! FUNCTION
707 : !! Compute alpha * mat^T in place. target: double precision real matrix.
708 : !!
709 : !! INPUTS
710 : !! n=size of the matrix
711 : !! [alpha]=scalar, set to 1.0 if not present
712 : !!
713 : !! SIDE EFFECTS
714 : !! mat(n,n)=in output, it contains alpha * mat^T.
715 : !!
716 : !! SOURCE
717 :
718 0 : subroutine sqmat_itranspose_dp(n,mat,alpha)
719 :
720 : !Arguments ------------------------------------
721 : !scalars
722 : integer,intent(in) :: n
723 : real(dp),optional,intent(in) :: alpha
724 : !arrays
725 : real(dp),intent(inout) :: mat(n,n)
726 : ! *************************************************************************
727 :
728 : #ifdef HAVE_LINALG_MKL_IMATCOPY
729 0 : if (PRESENT(alpha)) then
730 0 : call mkl_dimatcopy("Column", "Trans", n, n, alpha, mat, n, n)
731 : else
732 0 : call mkl_dimatcopy("Column", "Trans", n, n, one_dp, mat, n, n)
733 : end if
734 : #else
735 : ! Fallback to Fortran.
736 : if (PRESENT(alpha)) then
737 : mat = alpha * TRANSPOSE(mat)
738 : else
739 : mat = TRANSPOSE(mat)
740 : end if
741 : #endif
742 :
743 0 : end subroutine sqmat_itranspose_dp
744 : !!***
745 :
746 : !----------------------------------------------------------------------
747 :
748 : !!****f* m_hide_blas/sqmat_itranspose_spc
749 : !! NAME
750 : !! sqmat_itranspose_spc
751 : !!
752 : !! FUNCTION
753 : !! Compute alpha * mat^T in place. target: single precision complex matrix.
754 : !!
755 : !! INPUTS
756 : !! n=size of the matrix
757 : !! [alpha]=scalar, set to 1.0 if not present
758 : !!
759 : !! SIDE EFFECTS
760 : !! mat(n,n)=in output, it contains alpha * mat^T.
761 : !!
762 : !! SOURCE
763 :
764 0 : subroutine sqmat_itranspose_spc(n,mat,alpha)
765 :
766 : !Arguments ------------------------------------
767 : !scalars
768 : integer,intent(in) :: n
769 : complex(sp),optional,intent(in) :: alpha
770 : !arrays
771 : complex(sp),intent(inout) :: mat(n,n)
772 : ! *************************************************************************
773 :
774 : #ifdef HAVE_LINALG_MKL_IMATCOPY
775 0 : if (PRESENT(alpha)) then
776 0 : call mkl_cimatcopy("Column", "Trans", n, n, alpha, mat, n, n)
777 : else
778 0 : call mkl_cimatcopy("Column", "Trans", n, n, cone_sp, mat, n, n)
779 : end if
780 : #else
781 : ! Fallback to Fortran.
782 : if (PRESENT(alpha)) then
783 : mat = alpha * TRANSPOSE(mat)
784 : else
785 : mat = TRANSPOSE(mat)
786 : end if
787 : #endif
788 :
789 0 : end subroutine sqmat_itranspose_spc
790 : !!***
791 :
792 : !----------------------------------------------------------------------
793 :
794 : !!****f* m_hide_blas/sqmat_itranspose_dpc
795 : !! NAME
796 : !! sqmat_itranspose_dpc
797 : !!
798 : !! FUNCTION
799 : !! Compute alpha * mat^T in place. target: double precision complex matrix.
800 : !!
801 : !! INPUTS
802 : !! n=size of the matrix
803 : !! [alpha]=scalar, set to 1.0 if not present
804 : !!
805 : !! SIDE EFFECTS
806 : !! mat(n,n)=in output, it contains alpha * mat^T.
807 : !!
808 : !! SOURCE
809 :
810 2719 : subroutine sqmat_itranspose_dpc(n,mat,alpha)
811 :
812 : !Arguments ------------------------------------
813 : !scalars
814 : integer,intent(in) :: n
815 : complex(dp),optional,intent(in) :: alpha
816 : !arrays
817 : complex(dp),intent(inout) :: mat(n,n)
818 : ! *************************************************************************
819 :
820 : #ifdef HAVE_LINALG_MKL_IMATCOPY
821 2719 : if (PRESENT(alpha)) then
822 0 : call mkl_zimatcopy("Column", "Trans", n, n, alpha, mat, n, n)
823 : else
824 2719 : call mkl_zimatcopy("Column", "Trans", n, n, cone_dpc, mat, n, n)
825 : end if
826 : #else
827 : ! Fallback to Fortran.
828 : if (PRESENT(alpha)) then
829 : mat = alpha * TRANSPOSE(mat)
830 : else
831 : mat = TRANSPOSE(mat)
832 : end if
833 : #endif
834 :
835 2719 : end subroutine sqmat_itranspose_dpc
836 : !!***
837 :
838 : !----------------------------------------------------------------------
839 :
840 : !!****f* m_hide_blas/sqmat_otranspose_sp
841 : !! NAME
842 : !! sqmat_otranspose_sp
843 : !!
844 : !! FUNCTION
845 : !! Compute alpha * mat^T out-of-place. target: single precision real matrix.
846 : !!
847 : !! INPUTS
848 : !! n=size of the matrix
849 : !! [alpha]=scalar, set to 1.0 if not present
850 : !! imat(n,n)=Input matrix.
851 : !!
852 : !! OUTPUT
853 : !! omat(n,n)=contains alpha * imat^T.
854 : !!
855 : !! SOURCE
856 :
857 0 : subroutine sqmat_otranspose_sp(n,imat,omat,alpha)
858 :
859 : !Arguments ------------------------------------
860 : !scalars
861 : integer,intent(in) :: n
862 : real(sp),optional,intent(in) :: alpha
863 : !arrays
864 : real(sp),intent(in) :: imat(n,n)
865 : real(sp),intent(out) :: omat(n,n)
866 : ! *************************************************************************
867 :
868 : #ifdef HAVE_LINALG_MKL_OMATCOPY
869 0 : if (PRESENT(alpha)) then
870 0 : call mkl_somatcopy("Column", "Transpose", n, n, alpha, imat, n, omat, n)
871 : else
872 0 : call mkl_somatcopy("Column", "Transpose", n, n, one_sp, imat, n, omat, n)
873 : end if
874 : #else
875 : ! Fallback to Fortran.
876 : if (PRESENT(alpha)) then
877 : omat = alpha * TRANSPOSE(imat)
878 : else
879 : omat = TRANSPOSE(imat)
880 : end if
881 : #endif
882 :
883 0 : end subroutine sqmat_otranspose_sp
884 : !!***
885 :
886 : !----------------------------------------------------------------------
887 :
888 : !!****f* m_hide_blas/sqmat_otranspose_dp
889 : !! NAME
890 : !! sqmat_otranspose_dp
891 : !!
892 : !! FUNCTION
893 : !! Compute alpha * mat^T out-of-place. target: double precision real matrix.
894 : !!
895 : !! INPUTS
896 : !! n=size of the matrix
897 : !! [alpha]=scalar, set to 1.0 if not present
898 : !! imat(n,n)=Input matrix.
899 : !!
900 : !! OUTPUT
901 : !! omat(n,n)=contains alpha * imat^T.
902 : !!
903 : !! SOURCE
904 :
905 0 : subroutine sqmat_otranspose_dp(n,imat,omat,alpha)
906 :
907 : !Arguments ------------------------------------
908 : !scalars
909 : integer,intent(in) :: n
910 : real(dp),optional,intent(in) :: alpha
911 : !arrays
912 : real(dp),intent(in) :: imat(n,n)
913 : real(dp),intent(out) :: omat(n,n)
914 : ! *************************************************************************
915 :
916 : #ifdef HAVE_LINALG_MKL_OMATCOPY
917 0 : if (PRESENT(alpha)) then
918 0 : call mkl_domatcopy("Column", "Transpose", n, n, alpha, imat, n, omat, n)
919 : else
920 0 : call mkl_domatcopy("Column", "Transpose", n, n, one_dp, imat, n, omat, n)
921 : end if
922 : #else
923 : ! Fallback to Fortran.
924 : if (PRESENT(alpha)) then
925 : omat = alpha * TRANSPOSE(imat)
926 : else
927 : omat = TRANSPOSE(imat)
928 : end if
929 : #endif
930 :
931 0 : end subroutine sqmat_otranspose_dp
932 : !!***
933 :
934 : !----------------------------------------------------------------------
935 :
936 : !!****f* m_hide_blas/sqmat_otranspose_spc
937 : !! NAME
938 : !! sqmat_otranspose_spc
939 : !!
940 : !! FUNCTION
941 : !! Compute alpha * mat^T out-of-place. target: single precision complex matrix.
942 : !!
943 : !! INPUTS
944 : !! n=size of the matrix
945 : !! [alpha]=scalar, set to 1.0 if not present
946 : !! imat(n,n)=Input matrix.
947 : !!
948 : !! OUTPUT
949 : !! omat(n,n)=contains alpha * imat^T.
950 : !!
951 : !! SOURCE
952 :
953 0 : subroutine sqmat_otranspose_spc(n,imat,omat,alpha)
954 :
955 : !Arguments ------------------------------------
956 : !scalars
957 : integer,intent(in) :: n
958 : complex(sp),optional,intent(in) :: alpha
959 : !arrays
960 : complex(sp),intent(in) :: imat(n,n)
961 : complex(sp),intent(out) :: omat(n,n)
962 : ! *************************************************************************
963 :
964 : #ifdef HAVE_LINALG_MKL_OMATCOPY
965 0 : if (PRESENT(alpha)) then
966 0 : call mkl_comatcopy("Column", "Transpose", n, n, alpha, imat, n, omat, n)
967 : else
968 0 : call mkl_comatcopy("Column", "Transpose", n, n, cone_sp, imat, n, omat, n)
969 : end if
970 : #else
971 : ! Fallback to Fortran.
972 : if (PRESENT(alpha)) then
973 : omat = alpha * TRANSPOSE(imat)
974 : else
975 : omat = TRANSPOSE(imat)
976 : end if
977 : #endif
978 :
979 0 : end subroutine sqmat_otranspose_spc
980 : !!***
981 :
982 : !----------------------------------------------------------------------
983 :
984 : !!****f* m_hide_blas/sqmat_otranspose_dpc
985 : !! NAME
986 : !! sqmat_otranspose_dpc
987 : !!
988 : !! FUNCTION
989 : !! Compute alpha * mat^T out-of-place. target: double precision complex matrix.
990 : !!
991 : !! INPUTS
992 : !! n=size of the matrix
993 : !! [alpha]=scalar, set to 1.0 if not present
994 : !! imat(n,n)=Input matrix.
995 : !!
996 : !! OUTPUT
997 : !! omat(n,n)=contains alpha * imat^T.
998 : !!
999 : !! SOURCE
1000 :
1001 0 : subroutine sqmat_otranspose_dpc(n,imat,omat,alpha)
1002 :
1003 : !Arguments ------------------------------------
1004 : !scalars
1005 : integer,intent(in) :: n
1006 : complex(dp),optional,intent(in) :: alpha
1007 : !arrays
1008 : complex(dp),intent(in) :: imat(n,n)
1009 : complex(dp),intent(out) :: omat(n,n)
1010 : ! *************************************************************************
1011 :
1012 : #ifdef HAVE_LINALG_MKL_OMATCOPY
1013 0 : if (PRESENT(alpha)) then
1014 0 : call mkl_zomatcopy("Column", "Transpose", n, n, alpha, imat, n, omat, n)
1015 : else
1016 0 : call mkl_zomatcopy("Column", "Transpose", n, n, cone_dpc, imat, n, omat, n)
1017 : end if
1018 : #else
1019 : ! Fallback to Fortran.
1020 : if (PRESENT(alpha)) then
1021 : omat = alpha * TRANSPOSE(imat)
1022 : else
1023 : omat = TRANSPOSE(imat)
1024 : end if
1025 : #endif
1026 :
1027 0 : end subroutine sqmat_otranspose_dpc
1028 : !!***
1029 :
1030 : !----------------------------------------------------------------------
1031 :
1032 : !!****f* m_hide_blas/sqmat_iconjgtrans_spc
1033 : !! NAME
1034 : !! sqmat_iconjgtrans_spc
1035 : !!
1036 : !! FUNCTION
1037 : !! Compute alpha * CONJG(mat^T) in place. target: single precision complex matrix.
1038 : !!
1039 : !! INPUTS
1040 : !! n=size of the matrix
1041 : !! [alpha]=scalar, set to 1.0 if not present
1042 : !!
1043 : !! SIDE EFFECTS
1044 : !! mat(n,n)=in output, it contains alpha * CONJG(mat^T).
1045 : !!
1046 : !! SOURCE
1047 :
1048 0 : subroutine sqmat_iconjgtrans_spc(n,mat,alpha)
1049 :
1050 : !Arguments ------------------------------------
1051 : !scalars
1052 : integer,intent(in) :: n
1053 : complex(sp),optional,intent(in) :: alpha
1054 : !arrays
1055 : complex(sp),intent(inout) :: mat(n,n)
1056 : ! *************************************************************************
1057 :
1058 : #ifdef HAVE_LINALG_MKL_IMATCOPY
1059 0 : if (PRESENT(alpha)) then
1060 0 : call mkl_cimatcopy("Column", "C", n, n, alpha, mat, n, n)
1061 : else
1062 0 : call mkl_cimatcopy("Column", "C", n, n, cone_sp, mat, n, n)
1063 : end if
1064 : #else
1065 : ! Fallback to Fortran.
1066 : if (PRESENT(alpha)) then
1067 : mat = alpha * TRANSPOSE(CONJG(mat))
1068 : else
1069 : mat = TRANSPOSE(CONJG(mat))
1070 : end if
1071 : #endif
1072 :
1073 0 : end subroutine sqmat_iconjgtrans_spc
1074 : !!***
1075 :
1076 : !----------------------------------------------------------------------
1077 :
1078 : !!****f* m_hide_blas/sqmat_iconjgtrans_dpc
1079 : !! NAME
1080 : !! sqmat_iconjgtrans_dpc
1081 : !!
1082 : !! FUNCTION
1083 : !! Compute alpha * CONJG(mat^T) in place. target: double precision complex matrix.
1084 : !!
1085 : !! INPUTS
1086 : !! n=size of the matrix
1087 : !! [alpha]=scalar, set to 1.0 if not present
1088 : !!
1089 : !! SIDE EFFECTS
1090 : !! mat(n,n)=in output, it contains alpha * CONJG(mat^T).
1091 : !!
1092 : !! SOURCE
1093 :
1094 0 : subroutine sqmat_iconjgtrans_dpc(n, mat, alpha)
1095 :
1096 : !Arguments ------------------------------------
1097 : !scalars
1098 : integer,intent(in) :: n
1099 : complex(dp),optional,intent(in) :: alpha
1100 : !arrays
1101 : complex(dp),intent(inout) :: mat(n,n)
1102 : ! *************************************************************************
1103 :
1104 : #ifdef HAVE_LINALG_MKL_IMATCOPY
1105 0 : if (PRESENT(alpha)) then
1106 0 : call mkl_zimatcopy("Column", "C", n, n, alpha, mat, n, n)
1107 : else
1108 0 : call mkl_zimatcopy("Column", "C", n, n, cone_dpc, mat, n, n)
1109 : end if
1110 : #else
1111 : ! Fallback to Fortran.
1112 : if (PRESENT(alpha)) then
1113 : mat = alpha * TRANSPOSE(CONJG(mat))
1114 : else
1115 : mat = TRANSPOSE(CONJG(mat))
1116 : end if
1117 : #endif
1118 :
1119 0 : end subroutine sqmat_iconjgtrans_dpc
1120 : !!***
1121 :
1122 : !----------------------------------------------------------------------
1123 :
1124 : !!****f* m_hide_blas/sqmat_oconjgtrans_spc
1125 : !! NAME
1126 : !! sqmat_oconjgtrans_spc
1127 : !!
1128 : !! FUNCTION
1129 : !! Compute alpha * CONJG(mat^T) out-of-place. target: single precision complex matrix.
1130 : !!
1131 : !! INPUTS
1132 : !! n=size of the matrix
1133 : !! [alpha]=scalar, set to 1.0 if not present
1134 : !! imat(n,n)=Input matrix.
1135 : !!
1136 : !! OUTPUT
1137 : !! omat(n,n)=contains alpha * CONJG(imat^T).
1138 : !!
1139 : !! SOURCE
1140 :
1141 0 : subroutine sqmat_oconjgtrans_spc(n, imat, omat, alpha)
1142 :
1143 : !Arguments ------------------------------------
1144 : !scalars
1145 : integer,intent(in) :: n
1146 : complex(sp),optional,intent(in) :: alpha
1147 : !arrays
1148 : complex(sp),intent(in) :: imat(n,n)
1149 : complex(sp),intent(out) :: omat(n,n)
1150 : ! *************************************************************************
1151 :
1152 : #ifdef HAVE_LINALG_MKL_OMATCOPY
1153 0 : if (PRESENT(alpha)) then
1154 0 : call mkl_comatcopy("Column", "C", n, n, alpha, imat, n, omat, n)
1155 : else
1156 0 : call mkl_comatcopy("Column", "C", n, n, cone_sp, imat, n, omat, n)
1157 : end if
1158 : #else
1159 : ! Fallback to Fortran.
1160 : if (PRESENT(alpha)) then
1161 : omat = alpha * TRANSPOSE(CONJG(imat))
1162 : else
1163 : omat = TRANSPOSE(CONJG(imat))
1164 : end if
1165 : #endif
1166 :
1167 0 : end subroutine sqmat_oconjgtrans_spc
1168 : !!***
1169 :
1170 : !----------------------------------------------------------------------
1171 :
1172 : !!****f* m_hide_blas/sqmat_oconjgtrans_dpc
1173 : !! NAME
1174 : !! sqmat_oconjgtrans_dpc
1175 : !!
1176 : !! FUNCTION
1177 : !! Compute alpha * CONJG(mat^T) out-of-place. target: double precision complex matrix.
1178 : !!
1179 : !! INPUTS
1180 : !! n=size of the matrix
1181 : !! [alpha]=scalar, set to 1.0 if not present
1182 : !! imat(n,n)=Input matrix.
1183 : !!
1184 : !! OUTPUT
1185 : !! omat(n,n)=contains alpha * CONJG(imat^T).
1186 : !!
1187 : !! SOURCE
1188 :
1189 0 : subroutine sqmat_oconjgtrans_dpc(n,imat,omat,alpha)
1190 :
1191 : !Arguments ------------------------------------
1192 : !scalars
1193 : integer,intent(in) :: n
1194 : complex(dp),optional,intent(in) :: alpha
1195 : !arrays
1196 : complex(dp),intent(in) :: imat(n,n)
1197 : complex(dp),intent(out) :: omat(n,n)
1198 : ! *************************************************************************
1199 :
1200 : #ifdef HAVE_LINALG_MKL_OMATCOPY
1201 0 : if (PRESENT(alpha)) then
1202 0 : call mkl_zomatcopy("Column", "C", n, n, alpha, imat, n, omat, n)
1203 : else
1204 0 : call mkl_zomatcopy("Column", "C", n, n, cone_dpc, imat, n, omat, n)
1205 : end if
1206 : #else
1207 : ! Fallback to Fortran.
1208 : if (PRESENT(alpha)) then
1209 : omat = alpha * TRANSPOSE(CONJG(imat))
1210 : else
1211 : omat = TRANSPOSE(CONJG(imat))
1212 : end if
1213 : #endif
1214 :
1215 0 : end subroutine sqmat_oconjgtrans_dpc
1216 : !!***
1217 :
1218 : !----------------------------------------------------------------------
1219 :
1220 : end module m_hide_blas
1221 : !!***
|
