Skip to content

Commit

Permalink
Added intent(in, out, inout) and use :: only. to several modules
Browse files Browse the repository at this point in the history
  • Loading branch information
kousuke-nakano committed Jul 17, 2023
1 parent 87bca0b commit 49b72c8
Show file tree
Hide file tree
Showing 7 changed files with 134 additions and 44 deletions.
22 changes: 18 additions & 4 deletions src/m_common/upwinv.f90
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,15 @@
subroutine upwinv(nel, jel, indt, nelorb, winv, v, psi)
use constants, only: yes_ontarget
implicit none
integer nel, indt, i, j, jel, nelorb
real*8 winv(nel, indt), psi(indt, nel), v(nel)

! argument parameters
integer, intent(in) :: nel, indt, jel, nelorb
real*8, intent(in) :: psi(indt, nel), v(nel)
real*8, intent(inout) :: winv(nel, indt)

! local variables
integer i, j

#ifdef _OFFLOAD
!$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
#endif
Expand All @@ -43,11 +50,18 @@ subroutine upwinv(nel, jel, indt, nelorb, winv, v, psi)
end do
return
end

subroutine upwinv_complex(nel, jel, indt, nelorb, winv, v, psi)
use constants, only: yes_ontarget
implicit none
integer nel, indt, i, j, jel, nelorb
complex*16 winv(nel, indt), psi(indt, nel), v(nel)

! argument parameters
integer, intent(in) :: nel, indt, jel, nelorb
complex*16, intent(in) :: psi(indt, nel), v(nel)
complex*16, intent(inout) :: winv(nel, indt)

! local variables
integer i, j
#ifdef _OFFLOAD
!$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
#endif
Expand Down
82 changes: 58 additions & 24 deletions src/m_common/upwinvp.f90
Original file line number Diff line number Diff line change
Expand Up @@ -13,11 +13,19 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.

!> subroutine to update the matrix winv, for real(8)
subroutine upwinvp(nel, indt, winv, ainv, ainvn, psi)
use constants, only: yes_ontarget
implicit none
integer nel, indt, i, j
real*8 winv(nel, indt), psi(indt, nel, 2), ainv(nel), ainvn(nel)

! argument parameters
integer, intent(in) :: nel, indt
real*8, intent(in) :: psi(indt, nel, 2), ainv(nel), ainvn(nel)
real*8, intent(inout) :: winv(nel, indt)

! local variables
integer :: i, j

#ifdef _OFFLOAD
if (yes_ontarget) then
!$omp target teams distribute parallel do collapse(2)
Expand All @@ -42,12 +50,48 @@ subroutine upwinvp(nel, indt, winv, ainv, ainvn, psi)
#endif
return
end

!> subroutine to update the matrix winv, for complex(16)
subroutine upwinvp_complex(nel, indt, winv, ainv, ainvn, psi)
use constants, only: yes_ontarget
implicit none

! argument parameters
integer, intent(in) :: nel, indt
complex*16, intent(in) :: psi(indt, nel, 2), ainv(nel), ainvn(nel)
complex*16, intent(inout) :: winv(nel, indt)

! local variables
integer i, j

#ifdef _OFFLOAD
!$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
#endif
do i = 1, indt
do j = 1, nel
winv(j, i) = winv(j, i) + ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2)
end do
end do
#ifdef _OFFLOAD
!$omp end target teams distribute parallel do
#endif
return
end

!> subroutine to update the matrix winv of Pfaffian, for real(8)
subroutine upwinvp_pfaff(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, winvup, winvdo, psi, ainv)
use constants, only: yes_ontarget
implicit none
integer nelc, nelcdo, nelup, neldo, nmol, nmolipf, nmolshift, indt, i, j
real*8 winvup(nelup, indt), winvdo(neldo, indt), psi(nmolipf, indt), ainv(nmol)

! argument parameters
integer, intent(in) :: nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt
real*8, intent(in) :: psi(nmolipf, indt), ainv(nmol)
real*8, intent(inout) :: winvup(nelup, indt), winvdo(neldo, indt)

! local variables
integer i, j, nelcdo
real*8 csum

if (yes_ontarget) then
if (nelc .le. nelup) then
#ifdef _OFFLOAD
Expand Down Expand Up @@ -116,11 +160,19 @@ subroutine upwinvp_pfaff(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, win
end if
return
end

!> subroutine to update the matrix winv of Pfaffian, for complex(16)
subroutine upwinvp_pfaff_complex(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, winvup, winvdo, psi, ainv)
use constants, only: yes_ontarget
implicit none
integer nelc, nelcdo, nelup, neldo, nmol, nmolipf, nmolshift, indt, i, j
complex*16 winvup(nelup, indt), winvdo(neldo, indt), psi(nmolipf, indt), ainv(nmol)

! argument parameters
integer, intent(in) :: nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt
complex*16, intent(in) :: psi(nmolipf, indt), ainv(nmol)
complex*16, intent(inout) :: winvup(nelup, indt), winvdo(neldo, indt)

! local variables
integer i, j, nelcdo
complex*16 csum
if (yes_ontarget) then
if (nelc .le. nelup) then
Expand Down Expand Up @@ -194,21 +246,3 @@ subroutine upwinvp_pfaff_complex(nelc, nelup, neldo, nmol, nmolipf, nmolshift, i
end if
return
end
subroutine upwinvp_complex(nel, indt, winv, ainv, ainvn, psi)
use constants, only: yes_ontarget
implicit none
integer nel, indt, i, j
complex*16 winv(nel, indt), psi(indt, nel, 2), ainv(nel), ainvn(nel)
#ifdef _OFFLOAD
!$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
#endif
do i = 1, indt
do j = 1, nel
winv(j, i) = winv(j, i) + ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2)
end do
end do
#ifdef _OFFLOAD
!$omp end target teams distribute parallel do
#endif
return
end
8 changes: 7 additions & 1 deletion src/m_common/write_type_orb.f90
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,15 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.

!> This subroutine returns typeorb for jastrow orbitals
subroutine write_type_orb(nshellj, multij, ioccj, typeorb)
implicit none
integer nshellj, multij(*), ioccj(*), typeorb(*)

! argument parameters
integer, intent(in) :: nshellj, multij(*), ioccj(*)
integer, intent(out) :: typeorb(*)

! local variables
integer i, ii, j, ind_type

ii = 0
Expand Down
1 change: 1 addition & 0 deletions src/m_common/zgemm_my.f90
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.

!> This subroutine is a wrapper for the LAPACK zgemm routine
subroutine ZGEMM_MY(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA&
&, B, LDB, BETA, C, LDC, nproc, rank, comm_mpi)
implicit none
Expand Down
25 changes: 17 additions & 8 deletions src/m_common/zsktri.f90
Original file line number Diff line number Diff line change
Expand Up @@ -13,19 +13,29 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.

!> This subroutine manipulates a skewsymmetric matrix A.
subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
implicit none
character uplo
integer n, i, j, lda, ldinv, info, lwork, ipiv(*)
complex*16 a(lda, *), ainv(ldinv, *), work(*)

! argument parameters
character, intent(in) :: uplo
integer, intent(in) :: n, lda, ldinv, ipiv(*)
complex*16, intent(in out) :: a(lda, *), work(*)
complex*16, intent(out) :: ainv(ldinv, *)
integer, intent(out) :: info

! local variables
integer i, j
real*8 rcond
logical yeslap

! First factorization assumed
! lwork=3*n
! lwork=3*n
! CALL ZSKTRF( UPLO, 'N', N, A, LDA, IPIV, WORK, LWORK, INFO)
! Identity matrix
! ipiv dimension required 3n
! work dimension required n^2+12*n-2
info = 0
yeslap = .false. ! if false the homemade algorithm is done.
do i = 1, n
do j = 1, i - 1
Expand All @@ -36,13 +46,13 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
ainv(j, i) = dcmplx(0.d0, 0.d0)
end do
end do
! fisrt permutation of ainv
! fisrt permutation of ainv
do i = n, 1, -1
work(1:n) = ainv(i, 1:n)
ainv(i, 1:n) = ainv(ipiv(i), 1:n)
ainv(ipiv(i), 1:n) = work(1:n)
end do
! SECOND
! second
if (UPLO .eq. 'u' .or. UPLO .eq. 'U') then
do i = 1, N - 1
work(i) = a(i, i + 1)
Expand All @@ -63,7 +73,7 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
a(2:n, 1) = dcmplx(0.d0, 0.d0)
end if
work(2*N - 1:3*N - 2) = dcmplx(0.d0, 0.d0) ! diagonal elements of skew matrix , obviously set to zero.
call ZTRTRS(UPLO, 'N', 'U', N, N, A, LDA, ainv, LDINV, INFO)
call ZTRTRS(UPLO, 'N', 'U', N, N, A, LDA, ainv, LDINV, INFO) ! ainv = A^-1
! USE STANDARD LAPACK DGTSVX or the simples DGTSV?
if (yeslap) then
call ZGTSVX('N', 'N', N, N, work(N), work(2*N - 1), work, work(3*N)&
Expand All @@ -83,7 +93,6 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
ainv(i, 1:n) = ainv(ipiv(i), 1:n)
ainv(ipiv(i), 1:n) = work(1:n)
end do

return
end

19 changes: 15 additions & 4 deletions src/m_common/zsktrs.f90
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
! Copyright (C) 2022 TurboRVB group
! Copyright (C) 2023- TurboRVB group
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
Expand All @@ -13,11 +13,22 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.

!> This subroutine solves A X = B
!> To be deleted!
subroutine zsktrs(uplo, n, nhrs, a, b, ldb, x, info)
implicit none
integer n, nhrs, ldb, info, i, j
complex*16 a(*), b(ldb, nhrs), x(nhrs, *)
character uplo

! argument parameters
integer, intent(in) :: n, nhrs, ldb
integer, intent(out) :: info
complex*16, intent(in out) :: a(*), b(ldb, nhrs)
complex*16, intent(out) :: x(nhrs, *)
character, intent(in) :: uplo

! local variables
integer i, j

info = 0
if (uplo .eq. 'l' .or. uplo .eq. 'L') then
a(1:n - 1) = -a(1:n - 1)
end if
Expand Down
21 changes: 18 additions & 3 deletions src/m_common/zsktrs_qp.f90
Original file line number Diff line number Diff line change
Expand Up @@ -13,16 +13,31 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.

!> This subroutine solves A * X = B.
!> Our homemade algorithm to solve A * X = B (No transpose)
!> where A is a skewsymmetric matrix, X and B are rectanglar matrices.
!> The input vector a stores the elements of A such that
!> A(i+1,i)=a(i) and A(i,i+1)=-a(i) for i=1,..,n-1
!> a(n), b(ldb, n), x(n)
subroutine zsktrs(uplo, n, nhrs, a, b, ldb, x, info)
implicit none
integer n, nhrs, ldb, info, i, j
complex*16 a(*), b(ldb, *), x(*)

! argument variables
integer, intent(in) :: n, nhrs, ldb
integer, intent(out) :: info
complex*16, intent(inout) :: a(*), b(ldb, *)
complex*16, intent(out) :: x(*)
character, intent(in) :: uplo

! local variables
integer i, j
#ifdef __PORT
complex*16 aq, a2q, xo
#else
complex*32 aq, a2q, xo
#endif
character uplo

info = 0
if (uplo .eq. 'l' .or. uplo .eq. 'L') then
a(1:n - 1) = -a(1:n - 1)
end if
Expand Down

0 comments on commit 49b72c8

Please sign in to comment.