TABLE OF CONTENTS

ABINIT/m_hide_blas [ Modules ]

[ Top ] [ Modules ]

NAME

  m_hide_blas

FUNCTION

 This module defines interfaces for overloading BLAS routines.
 whose goal is twofold. On one hand, using generic interfaces renders
 the code more readable, especially when the routine can be compiled with
 different precision type (single-precision or double precision as done for example in the GW code)
 On the other hand, the generic interfaces defined here introduce a programming
 layer that can be exploited for interfacing non-standard libraries such as for
 example CUBLAS routines for GPU computations.

COPYRIGHT

 Copyright (C) 1992-2024 ABINIT group (MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

NOTES

 The convention about names of interfaced routine is: x<name>,
 where <name> is usually equal to the name of the standard routine
 without the first character specifying the type and kind.
 The full list of names is reported below.
 BLAS procedures interfaced in this module are marked with an asterisk.
 A complete list of possible overloaded interfaces is provided as guide for future additions.

 ================
 ==== BLAS 1 ====
 ================
 FUNCTION idamax isamax icamax izamax  ---> XIAMAX(n,dx,incx)
 * FUNCTION  snrm2  dnrm2 scnrm2 dznmr2  ---> XNRM2(n,x,incx)
 FUNCTION  sasum  dasum scasum dzasum  ---> XASUM(n,x,incx)
 * FUNCTION               cdotu  zdotu ---> XDOTU(n,x,incx,y,incy)
 * FUNCTION               cdotc  zdotc ---> XDOTC(n,x,incx,y,incy)
 FUNCTION  sdot   ddot                 ---> XDOT(n,x,incx,y,incy)
 FUNCTION  sdsdot sdot                 ---> XDSDOT(n,x,incx,y,incy)
 SUBROUTINE saxpy daxpy caxpy  zaxpy   ---> XAXPY(n,ca,cx,incx,cy,incy)
 * SUBROUTINE scopy dcopy ccopy  zcopy   ---> XCOPY(n,cx,incx,cy,incy)
 SUBROUTINE srotg drotg crotg  zrotg   ---> XROTG(a,b,c,s)
 SUBROUTINE srot  drot  csrot  zdrot   ---> XROT(n,x,incx,y,incy,c,s)
 * SUBROUTINE sscal dscal cscal  zscal
                        csscal zdscal  ---> XSCAL(n,a,x,incx)
 SUBROUTINE sswap dswap cswap  zswap   ---> XSWAP(n,x,incx,y,incy)

 ================
 ==== BLAS 2 ====
 ================
 SUBROUTINE sgbmv dgbmv cgbmv zgbmv    ---> XGBMV(trans,m,kl,ku,n,alpha,A,lda,X,incx,beta,Y,incy)
 * SUBROUTINE sgemv dgemv cgemv zgemv    ---> XGEMV(trans,m,n,alpha,A,lda,X,incx,beta,Y,incy)
 * SUBROUTINE             cgerc zgerc    ---> XGERC(m,n,alpha,x,incx,y,incy,A,lda)
 SUBROUTINE             cgeru zgeru    ---> XGERU(m,n,alpha,x,incx,y,incy,A,lda)
 SUBROUTINE             chbmv zhbmv    ---> XHBMV(uplo,n,k,alpha,A,lda,X,incx,beta,Y,incy)
 SUBROUTINE             chemv zhemv    ---> XHEMV(uplo,n,alpha,A,lda,X,incx,beta,Y,incy)
 * SUBROUTINE             cher  zher     ---> XHER(uplo,n,alpha,x,incx,A,lda)
 SUBROUTINE             cher2 zher2    ---> XHER2(uplo,n,alpha,x,incx,y,incy,A,lda)
 SUBROUTINE             chpr  zhpr     ---> XHPR(uplo,n,alpha,x,incx,AP)
 SUBROUTINE             chpr2 zhpr2    ---> XHPR2(uplo,n,alpha,x,incx,y,incy,AP)
 SUBROUTINE             chpmv zhpmv    ---> XHPMV(uplo,n,alpha,AP,X,incx,beta,Y,incy)
 SUBROUTINE stbmv dtbmv ctbmv ztbmv    ---> XTBMV(uplo,trans,diag,n,k,A,lda,X,incx)
 SUBROUTINE stpmv dtpmv ctpmv ztpmv    ---> XTPMV(uplo,trans,diag,n,AP,X,incx)
 SUBROUTINE strmv dtrmv ctrmv ztrmv    ---> XTRMV(uplo,trans,diag,n,A,lda,X,incx)
 SUBROUTINE ssymv dsymv                ---> XSYMV(uplo,n,alpha,A,lda,X,incx,beta,Y,incy)
 SUBROUTINE ssbmv dsbmv                ---> XSBMV(uplo,n,k,alpha,A,lda,X,incx,beta,Y,incy)
 SUBROUTINE sspmv dspmv                ---> XSPMV(uplo,n,alpha,AP,X,incx,beta,Y,incy)
 SUBROUTINE stbsv dtbsv ctbsv ztbsv    ---> XTBSV(uplo,trans,diag,n,k,A,lda,X,incx)
 SUBROUTINE stpsv dtpsv ctpsv ztpsv    ---> XTPSV(uplo,trans,diag,n,AP,X,incx)
 SUBROUTINE strsv dtrsv ctrsv ztrsv    ---> XTRSV(uplo,trans,diag,n,A,lda,X,incx)
 SUBROUTINE  sger  dger                ---> XGER(m,n,alpha,x,incx,y,incy,A,lda)
 SUBROUTINE  sspr  dspr                ---> XSPR(uplo,n,alpha,x,incx,AP)
 SUBROUTINE sspr2 dspr2                ---> XSPR2(uplo,n,alpha,x,incx,y,incy,AP)
 SUBROUTINE  ssyr  dsyr                ---> XSYR(uplo,n,alpha,x,incx,A,lda)
 SUBROUTINE ssyr2 dsyr2                ---> XSYR2(uplo,n,alpha,x,incx,y,incy,A,lda)

 ================
 ==== BLAS 3 ====
 ================
 * SUBROUTINE sgemm dgemm cgemm zgemm      ---> XGEMM(transA,transB,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc)
 SUBROUTINE             chemm zhemm      ---> XHEMM(side,uplo,m,n,alpha,A,lda,B,ldb,beta,C,ldc)
 SUBROUTINE            cher2k zher2k     ---> XHER2K(uplo,trans,n,k,alpha,A,lda,B,ldb,beta,C,ldc)
 * SUBROUTINE             cherk zherk      ---> XHERK(uplo,trans,n,k,alpha,A,lda,beta,C,ldc)
 SUBROUTINE ssymm dsymm csymm zsymm      ---> XSYMM(side,uplo,m,n,alpha,A,lda,B,ldb,beta,C,ldc)
 SUBROUTINE ssyr2k dsyr2k csyr2k zsyr2k  ---> XSYR2K(uplo,trans,n,k,alpha,A,lda,B,ldb,beta,C,ldc)
 SUBROUTINE ssyrk dsyrk csyrk zsyrk      ---> XSYRK(uplo,trans,n,k,alpha,A,lda,beta,C,ldc)
 SUBROUTINE strmm dtrmm ctrmm ztrmm      ---> XTRMM(side,uplo,transa,diag,m,n,alpha,A,lda,B,ldb)
 SUBROUTINE strsm dtrsm ctrsm ztrsm      ---> XTRSM(side,uplo,transa,diag,m,n,alpha,A,lda,B,ldb)

SOURCE

 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(spc),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(dpc),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     implicit none
175     integer :: incx
176     integer :: n
177     real(sp) :: sa
178     real(sp) :: sx(*)
179   end subroutine sscal
180   !
181   subroutine dscal(n,da,dx,incx)
182     use defs_basis
183     implicit none
184     integer :: incx
185     integer :: n
186     real(dp):: da
187     real(dp):: dx(*)
188   end subroutine dscal
189   !
190   subroutine cscal(n,ca,cx,incx)
191     use defs_basis
192     implicit none
193     integer :: incx
194     integer :: n
195     complex(spc) :: ca
196     complex(spc) :: cx(*)
197   end subroutine cscal
198   !
199   subroutine zscal(n,za,zx,incx)
200     use defs_basis
201     implicit none
202     integer :: incx
203     integer :: n
204     complex(dpc) :: za
205     complex(dpc) :: zx(*)
206   end subroutine zscal
207   !
208   subroutine  csscal(n,sa,cx,incx)
209     use defs_basis
210     implicit none
211     integer :: incx
212     integer :: n
213     real(sp) :: sa
214     complex(spc) :: cx(*)
215   end subroutine csscal
216   !
217   subroutine  zdscal(n,da,zx,incx)
218     use defs_basis
219     implicit none
220     integer :: incx
221     integer :: n
222     real(dp) :: da
223     complex(dpc) :: zx(*)
224   end subroutine zdscal
225   !
226 end interface xscal
227 
228 !-------------------------------------------------------------------------------
229 
230 interface xdotu
231   !
232 #ifdef HAVE_LINALG_ZDOTU_BUG
233   module procedure cdotu
234   module procedure zdotu
235 #else
236   function cdotu(n,cx,incx,cy,incy)
237     use defs_basis
238     complex(spc) :: cdotu
239     complex(spc),intent(in) :: cx(*),cy(*)
240     integer,intent(in) :: incx,incy,n
241   end function cdotu
242   !
243   function zdotu(n,zx,incx,zy,incy)
244     use defs_basis
245     complex(dpc) :: zdotu
246     complex(dpc),intent(in) :: zx(*),zy(*)
247     integer,intent(in) :: incx,incy,n
248   end function zdotu
249 #endif
250   !
251 end interface xdotu
252 
253 !-------------------------------------------------------------------------------
254 
255 
256 ! CDOTC, CDOTU, ZDOTC, and ZDOTU are problematic if Mac OS X's Vec lib is used.
257 ! See http://developer.apple.com/hardwaredrivers/ve/errata.html.
258 ! If needed, we replace them with plain Fortran code.
259 
260 interface xdotc
261   !
262 #ifdef HAVE_LINALG_ZDOTC_BUG
263    module procedure cdotc
264    module procedure zdotc
265 #else
266   function cdotc(n,cx,incx,cy,incy)
267     use defs_basis
268     complex(spc) :: cdotc
269     complex(spc),intent(in) :: cx(*),cy(*)
270     integer,intent(in) :: incx,incy,n
271   end function cdotc
272   !
273   function zdotc(n,zx,incx,zy,incy)
274     use defs_basis
275     complex(dpc) :: zdotc
276     complex(dpc),intent(in) :: zx(*),zy(*)
277     integer,intent(in) :: incx,incy,n
278   end function zdotc
279 #endif
280   !
281 end interface xdotc
282 
283 !-------------------------------------------------------------------------------
284 
285 interface xcopy
286    !module procedure ABI_xcopy
287  !
288  subroutine scopy(n,sx,incx,sy,incy)
289    use defs_basis
290    implicit none
291    integer,intent(in) :: incx
292    integer,intent(in) :: incy
293    integer,intent(in) :: n
294    real(sp),intent(in) ::  sx(*)
295    real(sp),intent(inout) :: sy(*)
296  end subroutine scopy
297  !
298  subroutine  dcopy(n,dx,incx,dy,incy)
299    use defs_basis
300    implicit none
301    integer,intent(in) :: incx
302    integer,intent(in) :: incy
303    integer,intent(in) :: n
304    real(dp),intent(in) :: dx(*)
305    real(dp),intent(inout) :: dy(*)
306  end subroutine dcopy
307  !
308  subroutine  ccopy(n,cx,incx,cy,incy)
309    use defs_basis
310    implicit none
311    integer,intent(in) :: incx
312    integer,intent(in) :: incy
313    integer,intent(in) :: n
314    complex(spc),intent(in) :: cx(*)
315    complex(spc),intent(inout) :: cy(*)
316  end subroutine ccopy
317  !
318  subroutine  zcopy(n,cx,incx,cy,incy)
319    use defs_basis
320    implicit none
321    integer,intent(in) :: incx
322    integer,intent(in) :: incy
323    integer,intent(in) :: n
324    complex(dpc),intent(in) :: cx(*)
325    complex(dpc),intent(inout) :: cy(*)
326  end subroutine zcopy
327  !
328 end interface xcopy
329 
330 !-------------------------------------------------------------------------------
331 
332 interface xgemv
333   !
334   subroutine sgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
335     use defs_basis
336     real(sp),intent(in) :: alpha, beta
337     integer,intent(in) :: incx, incy, lda, m, n
338     character(len=1),intent(in) :: trans
339     real(sp),intent(in) :: a( lda, * ), x( * )
340     real(sp),intent(inout) :: y( * )
341   end subroutine sgemv
342   !
343   subroutine dgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
344     use defs_basis
345     real(dp),intent(in) :: alpha, beta
346     integer,intent(in) :: incx, incy, lda, m, n
347     character(len=1),intent(in) :: trans
348     real(dp),intent(in) :: a( lda, * ), x( * )
349     real(dp),intent(inout) :: y( * )
350   end subroutine dgemv
351   !
352   subroutine cgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
353     use defs_basis
354     complex(spc),intent(in) :: alpha, beta
355     integer,intent(in) :: incx, incy, lda, m, n
356     character(len=1),intent(in) :: trans
357     complex(spc),intent(in) :: a( lda, * ), x( * )
358     complex(spc),intent(inout) :: y( * )
359   end subroutine cgemv
360   !
361   subroutine zgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
362     use defs_basis
363     complex(dpc),intent(in) :: alpha, beta
364     integer,intent(in) :: incx, incy, lda, m, n
365     character(len=1),intent(in) :: trans
366     complex(dpc),intent(in) :: a( lda, * ), x( * )
367     complex(dpc),intent(inout) :: y( * )
368   end subroutine zgemv
369   !
370 end interface xgemv
371 
372 !-------------------------------------------------------------------------------
373 
374 interface xgerc
375   !
376   subroutine cgerc ( m, n, alpha, x, incx, y, incy, a, lda )
377     use defs_basis
378     complex(spc),intent(in) :: alpha
379     integer,intent(in) :: incx, incy, lda, m, n
380     complex(spc),intent(inout) :: a( lda, * )
381     complex(spc),intent(in) :: x( * ), y( * )
382   end subroutine cgerc
383   !
384   subroutine zgerc ( m, n, alpha, x, incx, y, incy, a, lda )
385     use defs_basis
386     complex(dpc),intent(in) :: alpha
387     integer,intent(in) :: incx, incy, lda, m, n
388     complex(dpc),intent(inout) :: a( lda, * )
389     complex(dpc),intent(in) :: x( * ), y( * )
390   end subroutine zgerc
391   !
392 end interface xgerc
393 
394 !-------------------------------------------------------------------------------
395 
396 interface xher
397   !
398   subroutine cher ( uplo, n, alpha, x, incx, a, lda )
399     use defs_basis
400     character(len=1),intent(in) :: uplo
401     real(spc),intent(in) :: alpha
402     integer,intent(in) :: incx, lda, n
403     complex(spc),intent(inout) :: a( lda, * )
404     complex(spc),intent(in) :: x( * )
405   end subroutine cher
406   !
407   subroutine zher ( uplo, n, alpha, x, incx, a, lda )
408     use defs_basis
409     character(len=1),intent(in) :: uplo
410     real(dpc),intent(in) :: alpha
411     integer,intent(in) :: incx, lda, n
412     complex(dpc),intent(inout) :: a( lda, * )
413     complex(dpc),intent(in) :: x( * )
414   end subroutine zher
415   !
416 end interface xher
417 
418 !-------------------------------------------------------------------------------
419 
420 interface xgemm
421   !
422   subroutine sgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
423     use defs_basis
424     character(len=1),intent(in) :: transa, transb
425     integer,intent(in) :: m, n, k, lda, ldb, ldc
426     real(sp),intent(in) :: alpha, beta
427     real(sp),intent(in) :: a( lda, * ), b( ldb, * )
428     real(sp),intent(inout) :: c( ldc, * )
429   end subroutine sgemm
430   !
431   subroutine dgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
432     use defs_basis
433     character(len=1),intent(in) :: transa, transb
434     integer,intent(in) :: m, n, k, lda, ldb, ldc
435     real(dp),intent(in) :: alpha, beta
436     real(dp),intent(in) :: a( lda, * ), b( ldb, * )
437     real(dp),intent(inout) :: c( ldc, * )
438   end subroutine dgemm
439   !
440   subroutine cgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
441     use defs_basis
442     character(len=1),intent(in) :: transa, transb
443     integer,intent(in) :: m, n, k, lda, ldb, ldc
444     complex(spc),intent(in) :: alpha, beta
445     complex(spc),intent(in) :: a( lda, * ), b( ldb, * )
446     complex(spc),intent(inout) :: c( ldc, * )
447   end subroutine cgemm
448   !
449   subroutine zgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
450     use defs_basis
451     character(len=1),intent(in) :: transa, transb
452     integer,intent(in) :: m, n, k, lda, ldb, ldc
453     complex(dpc),intent(in) :: alpha, beta
454     complex(dpc),intent(in) :: a( lda, * ), b( ldb, * )
455     complex(dpc),intent(inout) :: c( ldc, * )
456   end subroutine zgemm
457   !
458 end interface xgemm
459 
460 interface xherk
461   !
462   subroutine cherk( uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
463     use defs_basis
464     character(len=1),intent(in) :: uplo
465     character(len=1),intent(in) :: trans
466     integer,intent(in) :: n,k,lda,ldc
467     real(sp),intent(in) :: alpha
468     complex(spc),intent(in) :: a( lda, * )
469     real(sp),intent(in) :: beta
470     complex(spc),intent(inout) :: c( ldc, * )
471   end subroutine cherk
472   !
473   subroutine zherk( uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
474     use defs_basis
475     character(len=1),intent(in) :: uplo
476     character(len=1),intent(in) :: trans
477     integer,intent(in) :: n,k,lda,ldc
478     real(dp),intent(in) :: alpha
479     complex(dpc),intent(in) :: a( lda, * )
480     real(dp),intent(in) :: beta
481     complex(dpc),intent(inout) :: c( ldc, * )
482   end subroutine zherk
483   !
484 end interface xherk
485 
486 !-------------------------------------------------------------------------------
487 
488 interface blas_cholesky_ortho
489   module procedure blas_cholesky_ortho_spc
490   module procedure blas_cholesky_ortho_dpc
491 end interface blas_cholesky_ortho
492 
493 interface sqmat_itranspose
494   module procedure sqmat_itranspose_sp
495   module procedure sqmat_itranspose_dp
496   module procedure sqmat_itranspose_spc
497   module procedure sqmat_itranspose_dpc
498 end interface sqmat_itranspose
499 
500 interface sqmat_otranspose
501   module procedure sqmat_otranspose_sp
502   module procedure sqmat_otranspose_dp
503   module procedure sqmat_otranspose_spc
504   module procedure sqmat_otranspose_dpc
505 end interface sqmat_otranspose
506 
507 interface sqmat_iconjgtrans
508   module procedure sqmat_iconjgtrans_spc
509   module procedure sqmat_iconjgtrans_dpc
510 end interface sqmat_iconjgtrans
511 
512 interface sqmat_oconjgtrans
513   module procedure sqmat_oconjgtrans_spc
514   module procedure sqmat_oconjgtrans_dpc
515 end interface sqmat_oconjgtrans
516 
517  real(dp),private,parameter ::  zero_dp = 0._dp
518  real(dp),private,parameter ::  one_dp  = 1._dp
519 
520  complex(dpc),private,parameter :: czero_dpc = (0._dp,0._dp)
521  complex(dpc),private,parameter :: cone_dpc  = (1._dp,0._dp)
522 
523 CONTAINS  !========================================================================================
524 
525 ! CDOTC, CDOTU, ZDOTC, and ZDOTU are problematic if Mac OS X's Vec lib is used.
526 ! See http://developer.apple.com/hardwaredrivers/ve/errata.html.
527 ! Here we replace them with plain Fortran code.
528 
529 #ifdef HAVE_LINALG_ZDOTC_BUG
530 !#warning "Using internal replacement for zdotc. External library cannot be used"
531 #include "replacements/cdotc.f"
532 #include "replacements/zdotc.f"
533 #endif
534 
535 #ifdef HAVE_LINALG_ZDOTU_BUG
536 !#warning "Using internal replacement for zdotu. External library cannot be used"
537 #include "replacements/cdotu.f"
538 #include "replacements/zdotu.f"
539 #endif
540 
541 !----------------------------------------------------------------------

m_hide_blas/blas_cholesky_ortho_dpc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  blas_cholesky_ortho_dpc

FUNCTION

  Performs the Cholesky orthonormalization of the vectors stored in iomat.

INPUTS

  vec_size=Size of each vector.
  nvec=Number of vectors in iomat

OUTPUT

  cf_ovlp=Cholesky factorization of the overlap matrix. ovlp = U^H U with U upper triangle matrix returned in cf_ovlp

SIDE EFFECTS

  iomat(vec_size,nvec)
    input: Input set of vectors.
    output: Orthonormalized set.

SOURCE

628 subroutine blas_cholesky_ortho_dpc(vec_size,nvec,iomat,cf_ovlp,use_gemm)
629 
630 !Arguments ------------------------------------
631  integer,intent(in) :: vec_size,nvec
632  logical,optional,intent(in) :: use_gemm
633  complex(dpc),intent(inout) :: iomat(vec_size,nvec)
634  complex(dpc),intent(out) :: cf_ovlp(nvec,nvec)
635 
636 !Local variables ------------------------------
637 !scalars
638  integer :: ierr
639  logical :: my_usegemm
640  character(len=500) :: msg
641 
642 ! *************************************************************************
643 
644  ! 1) Calculate overlap_ij =  <phi_i|phi_j>
645  my_usegemm = .FALSE.; if (PRESENT(use_gemm)) my_usegemm = use_gemm
646 
647  if (my_usegemm) then
648    call xgemm("Conjugate","Normal",nvec,nvec,vec_size,cone_dpc,iomat,vec_size,iomat,vec_size,czero_dpc,cf_ovlp,nvec)
649  else
650    call xherk("U","C", nvec, vec_size, one_dp, iomat, vec_size, zero_dp, cf_ovlp, nvec)
651  end if
652  !
653  ! 2) Cholesky factorization: ovlp = U^H U with U upper triangle matrix.
654  call ZPOTRF('U',nvec,cf_ovlp,nvec,ierr)
655  if (ierr/=0)  then
656    write(msg,'(a,i0)')' ZPOTRF returned info= ',ierr
657    ABI_ERROR(msg)
658  end if
659  !
660  ! 3) Solve X U = io_mat. On exit io_mat is orthonormalized.
661  call ZTRSM('Right','Upper','Normal','Normal',vec_size,nvec,cone_dpc,cf_ovlp,nvec,iomat,vec_size)
662 
663 end subroutine blas_cholesky_ortho_dpc

m_hide_blas/blas_cholesky_ortho_spc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  blas_cholesky_ortho_spc

FUNCTION

  Performs the Cholesky orthonormalization of the vectors stored in iomat.

INPUTS

  vec_size=Size of each vector.
  nvec=Number of vectors in iomat

OUTPUT

  cf_ovlp=Cholesky factorization of the overlap matrix. ovlp = U^H U with U upper triangle matrix returned in cf_ovlp

SIDE EFFECTS

  iomat(vec_size,nvec)
    input: Input set of vectors.
    output: Orthonormalized set.

SOURCE

566 subroutine blas_cholesky_ortho_spc(vec_size,nvec,iomat,cf_ovlp,use_gemm)
567 
568 !Arguments ------------------------------------
569  integer,intent(in) :: vec_size,nvec
570  logical,optional,intent(in) :: use_gemm
571  complex(spc),intent(inout) :: iomat(vec_size,nvec)
572  complex(spc),intent(out) :: cf_ovlp(nvec,nvec)
573 
574 !Local variables ------------------------------
575 !scalars
576  integer :: ierr
577  logical :: my_usegemm
578  character(len=500) :: msg
579 
580 ! *************************************************************************
581 
582  ! 1) Calculate overlap_ij =  <phi_i|phi_j>
583  ! TODO: use dsyrk
584  my_usegemm = .FALSE.; if (PRESENT(use_gemm)) my_usegemm = use_gemm
585 
586  if (my_usegemm) then
587    call xgemm("Conjugate","Normal",nvec,nvec,vec_size,cone_sp,iomat,vec_size,iomat,vec_size,czero_sp,cf_ovlp,nvec)
588  else
589    call xherk("U","C", nvec, vec_size, one_sp, iomat, vec_size, zero_sp, cf_ovlp, nvec)
590  end if
591  !
592  ! 2) Cholesky factorization: ovlp = U^H U with U upper triangle matrix.
593  call CPOTRF('U',nvec,cf_ovlp,nvec,ierr)
594  if (ierr/=0)  then
595    write(msg,'(a,i0)')' ZPOTRF returned info= ',ierr
596    ABI_ERROR(msg)
597  end if
598  !
599  ! 3) Solve X U = io_mat. On exit iomat is orthonormalized.
600  call CTRSM('Right','Upper','Normal','Normal',vec_size,nvec,cone_sp,cf_ovlp,nvec,iomat,vec_size)
601 
602 end subroutine blas_cholesky_ortho_spc

m_hide_blas/sqmat_iconjgtrans_dpc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_iconjgtrans_dpc

FUNCTION

  Compute alpha * CONJG(mat^T) in place. target: double precision complex matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present

SIDE EFFECTS

   mat(n,n)=in output, it contains alpha * CONJG(mat^T).

SOURCE

1115 subroutine sqmat_iconjgtrans_dpc(n, mat, alpha)
1116 
1117 !Arguments ------------------------------------
1118 !scalars
1119  integer,intent(in) :: n
1120  complex(dpc),optional,intent(in) :: alpha
1121 !arrays
1122  complex(dpc),intent(inout) :: mat(n,n)
1123 
1124 ! *************************************************************************
1125 
1126 #ifdef HAVE_LINALG_MKL_IMATCOPY
1127   if (PRESENT(alpha)) then
1128     call mkl_zimatcopy("Column", "C", n, n, alpha, mat, n, n)
1129   else
1130     call mkl_zimatcopy("Column", "C", n, n, cone_dpc, mat, n, n)
1131   end if
1132 #else
1133   ! Fallback to Fortran.
1134   if (PRESENT(alpha)) then
1135     mat = alpha * TRANSPOSE(CONJG(mat))
1136   else
1137     mat = TRANSPOSE(CONJG(mat))
1138   end if
1139 #endif
1140 
1141 end subroutine sqmat_iconjgtrans_dpc

m_hide_blas/sqmat_iconjgtrans_spc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_iconjgtrans_spc

FUNCTION

  Compute alpha * CONJG(mat^T) in place. target: single precision complex matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present

SIDE EFFECTS

   mat(n,n)=in output, it contains alpha * CONJG(mat^T).

SOURCE

1068 subroutine sqmat_iconjgtrans_spc(n,mat,alpha)
1069 
1070 !Arguments ------------------------------------
1071 !scalars
1072  integer,intent(in) :: n
1073  complex(spc),optional,intent(in) :: alpha
1074 !arrays
1075  complex(spc),intent(inout) :: mat(n,n)
1076 
1077 ! *************************************************************************
1078 
1079 #ifdef HAVE_LINALG_MKL_IMATCOPY
1080   if (PRESENT(alpha)) then
1081     call mkl_cimatcopy("Column", "C", n, n, alpha, mat, n, n)
1082   else
1083     call mkl_cimatcopy("Column", "C", n, n, cone_sp, mat, n, n)
1084   end if
1085 #else
1086   ! Fallback to Fortran.
1087   if (PRESENT(alpha)) then
1088     mat = alpha * TRANSPOSE(CONJG(mat))
1089   else
1090     mat = TRANSPOSE(CONJG(mat))
1091   end if
1092 #endif
1093 
1094 end subroutine sqmat_iconjgtrans_spc

m_hide_blas/sqmat_itranspose_dp [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_itranspose_dp

FUNCTION

  Compute alpha * mat^T in place. target: double precision real matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present

SIDE EFFECTS

   mat(n,n)=in output, it contains alpha * mat^T.

SOURCE

731 subroutine sqmat_itranspose_dp(n,mat,alpha)
732 
733 !Arguments ------------------------------------
734 !scalars
735  integer,intent(in) :: n
736  real(dp),optional,intent(in) :: alpha
737 !arrays
738  real(dp),intent(inout) :: mat(n,n)
739 
740 ! *************************************************************************
741 
742 #ifdef HAVE_LINALG_MKL_IMATCOPY
743   if (PRESENT(alpha)) then
744     call mkl_dimatcopy("Column", "Trans", n, n, alpha, mat, n, n)
745   else
746     call mkl_dimatcopy("Column", "Trans", n, n, one_dp, mat, n, n)
747   end if
748 #else
749   ! Fallback to Fortran.
750   if (PRESENT(alpha)) then
751     mat = alpha * TRANSPOSE(mat)
752   else
753     mat = TRANSPOSE(mat)
754   end if
755 #endif
756 
757 end subroutine sqmat_itranspose_dp

m_hide_blas/sqmat_itranspose_dpc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_itranspose_dpc

FUNCTION

  Compute alpha * mat^T in place. target: double precision complex matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present

SIDE EFFECTS

   mat(n,n)=in output, it contains alpha * mat^T.

SOURCE

825 subroutine sqmat_itranspose_dpc(n,mat,alpha)
826 
827 !Arguments ------------------------------------
828 !scalars
829  integer,intent(in) :: n
830  complex(dpc),optional,intent(in) :: alpha
831 !arrays
832  complex(dpc),intent(inout) :: mat(n,n)
833 
834 ! *************************************************************************
835 
836 #ifdef HAVE_LINALG_MKL_IMATCOPY
837   if (PRESENT(alpha)) then
838     call mkl_zimatcopy("Column", "Trans", n, n, alpha, mat, n, n)
839   else
840     call mkl_zimatcopy("Column", "Trans", n, n, cone_dpc, mat, n, n)
841   end if
842 #else
843   ! Fallback to Fortran.
844   if (PRESENT(alpha)) then
845     mat = alpha * TRANSPOSE(mat)
846   else
847     mat = TRANSPOSE(mat)
848   end if
849 #endif
850 
851 end subroutine sqmat_itranspose_dpc

m_hide_blas/sqmat_itranspose_sp [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_itranspose_sp

FUNCTION

  Compute alpha * mat^T in place. target: single precision real matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present

SIDE EFFECTS

   mat(n,n)=in output, it contains alpha * mat^T.

SOURCE

684 subroutine sqmat_itranspose_sp(n,mat,alpha)
685 
686 !Arguments ------------------------------------
687 !scalars
688  integer,intent(in) :: n
689  real(sp),optional,intent(in) :: alpha
690 !arrays
691  real(sp),intent(inout) :: mat(n,n)
692 
693 ! *************************************************************************
694 
695 #ifdef HAVE_LINALG_MKL_IMATCOPY
696   if (PRESENT(alpha)) then
697     call mkl_simatcopy("Column", "Trans", n, n, alpha, mat, n, n)
698   else
699     call mkl_simatcopy("Column", "Trans", n, n, one_sp, mat, n, n)
700   end if
701 #else
702   ! Fallback to Fortran.
703   if (PRESENT(alpha)) then
704     mat = alpha * TRANSPOSE(mat)
705   else
706     mat = TRANSPOSE(mat)
707   end if
708 #endif
709 
710 end subroutine sqmat_itranspose_sp

m_hide_blas/sqmat_itranspose_spc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_itranspose_spc

FUNCTION

  Compute alpha * mat^T in place. target: single precision complex matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present

SIDE EFFECTS

   mat(n,n)=in output, it contains alpha * mat^T.

SOURCE

778 subroutine sqmat_itranspose_spc(n,mat,alpha)
779 
780 !Arguments ------------------------------------
781 !scalars
782  integer,intent(in) :: n
783  complex(spc),optional,intent(in) :: alpha
784 !arrays
785  complex(spc),intent(inout) :: mat(n,n)
786 
787 ! *************************************************************************
788 
789 #ifdef HAVE_LINALG_MKL_IMATCOPY
790   if (PRESENT(alpha)) then
791     call mkl_cimatcopy("Column", "Trans", n, n, alpha, mat, n, n)
792   else
793     call mkl_cimatcopy("Column", "Trans", n, n, cone_sp, mat, n, n)
794   end if
795 #else
796   ! Fallback to Fortran.
797   if (PRESENT(alpha)) then
798     mat = alpha * TRANSPOSE(mat)
799   else
800     mat = TRANSPOSE(mat)
801   end if
802 #endif
803 
804 end subroutine sqmat_itranspose_spc

m_hide_blas/sqmat_oconjgtrans_dpc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_oconjgtrans_dpc

FUNCTION

  Compute alpha * CONJG(mat^T) out-of-place. target: double precision complex matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present
  imat(n,n)=Input matrix.

OUTPUT

  omat(n,n)=contains alpha * CONJG(imat^T).

SOURCE

1212 subroutine sqmat_oconjgtrans_dpc(n,imat,omat,alpha)
1213 
1214 !Arguments ------------------------------------
1215 !scalars
1216  integer,intent(in) :: n
1217  complex(dpc),optional,intent(in) :: alpha
1218 !arrays
1219  complex(dpc),intent(in) :: imat(n,n)
1220  complex(dpc),intent(out) :: omat(n,n)
1221 
1222 ! *************************************************************************
1223 
1224 #ifdef HAVE_LINALG_MKL_OMATCOPY
1225   if (PRESENT(alpha)) then
1226     call mkl_zomatcopy("Column", "C", n, n, alpha, imat, n, omat, n)
1227   else
1228     call mkl_zomatcopy("Column", "C", n, n, cone_dpc, imat, n, omat, n)
1229   end if
1230 #else
1231   ! Fallback to Fortran.
1232   if (PRESENT(alpha)) then
1233     omat = alpha * TRANSPOSE(CONJG(imat))
1234   else
1235     omat = TRANSPOSE(CONJG(imat))
1236   end if
1237 #endif
1238 
1239 end subroutine sqmat_oconjgtrans_dpc

m_hide_blas/sqmat_oconjgtrans_spc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_oconjgtrans_spc

FUNCTION

  Compute alpha * CONJG(mat^T) out-of-place. target: single precision complex matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present
  imat(n,n)=Input matrix.

OUTPUT

  omat(n,n)=contains alpha * CONJG(imat^T).

SOURCE

1163 subroutine sqmat_oconjgtrans_spc(n, imat, omat, alpha)
1164 
1165 !Arguments ------------------------------------
1166 !scalars
1167  integer,intent(in) :: n
1168  complex(spc),optional,intent(in) :: alpha
1169 !arrays
1170  complex(spc),intent(in) :: imat(n,n)
1171  complex(spc),intent(out) :: omat(n,n)
1172 
1173 ! *************************************************************************
1174 
1175 #ifdef HAVE_LINALG_MKL_OMATCOPY
1176   if (PRESENT(alpha)) then
1177     call mkl_comatcopy("Column", "C", n, n, alpha, imat, n, omat, n)
1178   else
1179     call mkl_comatcopy("Column", "C", n, n, cone_sp, imat, n, omat, n)
1180   end if
1181 #else
1182   ! Fallback to Fortran.
1183   if (PRESENT(alpha)) then
1184     omat = alpha * TRANSPOSE(CONJG(imat))
1185   else
1186     omat = TRANSPOSE(CONJG(imat))
1187   end if
1188 #endif
1189 
1190 end subroutine sqmat_oconjgtrans_spc

m_hide_blas/sqmat_otranspose_dp [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_otranspose_dp

FUNCTION

  Compute alpha * mat^T out-of-place. target: double precision real matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present
  imat(n,n)=Input matrix.

OUTPUT

  omat(n,n)=contains alpha * imat^T.

SOURCE

922 subroutine sqmat_otranspose_dp(n,imat,omat,alpha)
923 
924 !Arguments ------------------------------------
925 !scalars
926  integer,intent(in) :: n
927  real(dp),optional,intent(in) :: alpha
928 !arrays
929  real(dp),intent(in) :: imat(n,n)
930  real(dp),intent(out) :: omat(n,n)
931 
932 ! *************************************************************************
933 
934 #ifdef HAVE_LINALG_MKL_OMATCOPY
935   if (PRESENT(alpha)) then
936     call mkl_domatcopy("Column", "Transpose", n, n, alpha, imat, n, omat, n)
937   else
938     call mkl_domatcopy("Column", "Transpose", n, n, one_dp, imat, n, omat, n)
939   end if
940 #else
941   ! Fallback to Fortran.
942   if (PRESENT(alpha)) then
943     omat = alpha * TRANSPOSE(imat)
944   else
945     omat = TRANSPOSE(imat)
946   end if
947 #endif
948 
949 end subroutine sqmat_otranspose_dp

m_hide_blas/sqmat_otranspose_dpc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_otranspose_dpc

FUNCTION

  Compute alpha * mat^T out-of-place. target: double precision complex matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present
  imat(n,n)=Input matrix.

OUTPUT

  omat(n,n)=contains alpha * imat^T.

SOURCE

1020 subroutine sqmat_otranspose_dpc(n,imat,omat,alpha)
1021 
1022 !Arguments ------------------------------------
1023 !scalars
1024  integer,intent(in) :: n
1025  complex(dpc),optional,intent(in) :: alpha
1026 !arrays
1027  complex(dpc),intent(in) :: imat(n,n)
1028  complex(dpc),intent(out) :: omat(n,n)
1029 
1030 ! *************************************************************************
1031 
1032 #ifdef HAVE_LINALG_MKL_OMATCOPY
1033   if (PRESENT(alpha)) then
1034     call mkl_zomatcopy("Column", "Transpose", n, n, alpha, imat, n, omat, n)
1035   else
1036     call mkl_zomatcopy("Column", "Transpose", n, n, cone_dpc, imat, n, omat, n)
1037   end if
1038 #else
1039   ! Fallback to Fortran.
1040   if (PRESENT(alpha)) then
1041     omat = alpha * TRANSPOSE(imat)
1042   else
1043     omat = TRANSPOSE(imat)
1044   end if
1045 #endif
1046 
1047 end subroutine sqmat_otranspose_dpc

m_hide_blas/sqmat_otranspose_sp [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_otranspose_sp

FUNCTION

  Compute alpha * mat^T out-of-place. target: single precision real matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present
  imat(n,n)=Input matrix.

OUTPUT

  omat(n,n)=contains alpha * imat^T.

SOURCE

873 subroutine sqmat_otranspose_sp(n,imat,omat,alpha)
874 
875 !Arguments ------------------------------------
876 !scalars
877  integer,intent(in) :: n
878  real(sp),optional,intent(in) :: alpha
879 !arrays
880  real(sp),intent(in) :: imat(n,n)
881  real(sp),intent(out) :: omat(n,n)
882 
883 ! *************************************************************************
884 
885 #ifdef HAVE_LINALG_MKL_OMATCOPY
886   if (PRESENT(alpha)) then
887     call mkl_somatcopy("Column", "Transpose", n, n, alpha, imat, n, omat, n)
888   else
889     call mkl_somatcopy("Column", "Transpose", n, n, one_sp, imat, n, omat, n)
890   end if
891 #else
892   ! Fallback to Fortran.
893   if (PRESENT(alpha)) then
894     omat = alpha * TRANSPOSE(imat)
895   else
896     omat = TRANSPOSE(imat)
897   end if
898 #endif
899 
900 end subroutine sqmat_otranspose_sp

m_hide_blas/sqmat_otranspose_spc [ Functions ]

[ Top ] [ m_hide_blas ] [ Functions ]

NAME

  sqmat_otranspose_spc

FUNCTION

  Compute alpha * mat^T out-of-place. target: single precision complex matrix.

INPUTS

  n=size of the matrix
  [alpha]=scalar, set to 1.0 if not present
  imat(n,n)=Input matrix.

OUTPUT

  omat(n,n)=contains alpha * imat^T.

SOURCE

971 subroutine sqmat_otranspose_spc(n,imat,omat,alpha)
972 
973 !Arguments ------------------------------------
974 !scalars
975  integer,intent(in) :: n
976  complex(spc),optional,intent(in) :: alpha
977 !arrays
978  complex(spc),intent(in) :: imat(n,n)
979  complex(spc),intent(out) :: omat(n,n)
980 
981 ! *************************************************************************
982 
983 #ifdef HAVE_LINALG_MKL_OMATCOPY
984   if (PRESENT(alpha)) then
985     call mkl_comatcopy("Column", "Transpose", n, n, alpha, imat, n, omat, n)
986   else
987     call mkl_comatcopy("Column", "Transpose", n, n, cone_sp, imat, n, omat, n)
988   end if
989 #else
990   ! Fallback to Fortran.
991   if (PRESENT(alpha)) then
992     omat = alpha * TRANSPOSE(imat)
993   else
994     omat = TRANSPOSE(imat)
995   end if
996 #endif
997 
998 end subroutine sqmat_otranspose_spc