LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-03 20:54:00 +01:00
Code Issues Releases Wiki Activity

OpenMP

Browse Source
This commit is contained in:
Anthony Scemama 2017-09-25 22:34:56 +02:00
parent 987fc65984
commit 63af3aa6a2
1 changed files with 29 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
plugins/FourIdx/four_index.irp.f
View File

@ -30,7 +30,10 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB, &
integer(key_kind), allocatable :: key(:) integer(key_kind), allocatable :: key(:)
real(integral_kind), allocatable :: value(:) real(integral_kind), allocatable :: value(:)
ASSERT (k_start == i_start)
ASSERT (l_start == j_start)
ASSERT (a_start == c_start)
ASSERT (b_start == d_start)
i_min = min(i_start,a_start) i_min = min(i_start,a_start)
i_max = max(i_end ,a_end ) i_max = max(i_end ,a_end )
@ -50,42 +53,40 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB, &
ASSERT (LDB >= k_max) ASSERT (LDB >= k_max)
ASSERT (LDB >= l_max) ASSERT (LDB >= l_max)
!$OMP PARALLEL DEFAULT(PRIVATE) SHARED( &
!$OMP a_start,a_end,b_start,b_end,c_start,c_end,d_start,d_end,&
!$OMP i_start,i_end,j_start,j_end,k_start,k_end,l_start,l_end,&
!$OMP i_min,i_max,j_min,j_max,k_min,k_max,l_min,l_max, &
!$OMP map_a,map_c,matrix_B)
allocate( key(i_max*j_max*k_max), value(i_max*j_max*k_max) ) allocate( key(i_max*j_max*k_max), value(i_max*j_max*k_max) )
allocate( U(a_start:a_end, b_start:b_end, c_start:c_end) ) allocate( U(a_start:a_end, b_start:b_end, c_start:c_end) )
!$OMP DO
do d=d_start,d_end do d=d_start,d_end
U = 0.d0 U = 0.d0
print *, d print *, d
do l=l_start,l_end do l=1,l_end
if (dabs(matrix_B(l,d)) < 1.d-10) then if (dabs(matrix_B(l,d)) < 1.d-10) then
cycle cycle
endif endif
allocate( T(i_start:i_end, k_start:k_end, j_start:j_end) ) allocate( T(i_start:i_end, k_start:k_end, j_start:j_end), &
V(a_start:a_end, k_start:k_end, j_start:j_end) )
do k=k_start,k_end do k=k_start,k_end
do j=j_start,j_end do j=j_start,j_end
do i=i_start,i_end do i=i_start,k
call bielec_integrals_index(i,j,k,l,idx) call bielec_integrals_index(i,j,k,l,idx)
call map_get(map_a,idx,tmp) call map_get(map_a,idx,tmp)
T(i, k,j) = tmp T(i, k,j) = tmp
T(k, i,j) = tmp
enddo enddo
enddo enddo
enddo enddo
allocate( V(a_start:a_end, k_start:k_end, j_start:j_end) )
! V = 0.d0 call DGEMM('T','N', (a_end-a_start+1), &
! do a=a_start,a_end (k_end-k_start+1)*(j_end-j_start+1), &
! do k=k_start,k_end
! do j=j_start,j_end
! do i=i_start,i_end
! V(a, k,j) = V(a, k,j) + T(i, k,j)*matrix_B(i, a)
! enddo
! enddo
! enddo
! enddo
call DGEMM('T','N', (a_end-a_start+1), (k_end-k_start+1)*(j_end-j_start+1),&
(i_end-i_start+1), 1.d0, & (i_end-i_start+1), 1.d0, &
matrix_B(i_start,a_start), size(matrix_B,1), & matrix_B(i_start,a_start), size(matrix_B,1), &
T(i_start,k_start,j_start), size(T,1), 0.d0, & T(i_start,k_start,j_start), size(T,1), 0.d0, &
@ -94,50 +95,16 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB, &
deallocate(T) deallocate(T)
allocate( T(a_start:a_end, k_start:k_end, b_start:b_end) ) allocate( T(a_start:a_end, k_start:k_end, b_start:b_end) )
! V = 0.d0 call DGEMM('N','N', (a_end-a_start+1)*(k_end-k_start+1), &
! do a=a_start,a_end (b_end-b_start+1), &
! do k=k_start,k_end
! do b=b_start,b_end
! do j=j_start,j_end
! V(a,k, b) = V(a,k, b) + T(a,k, j)*matrix_B(j, b)
! enddo
! enddo
! enddo
! enddo
call DGEMM('N','N', (a_end-a_start+1)*(k_end-k_start+1),(b_end-b_start+1),&
(j_end-j_start+1), 1.d0, & (j_end-j_start+1), 1.d0, &
V(a_start,k_start,j_start), size(V,1)*size(V,2), & V(a_start,k_start,j_start), size(V,1)*size(V,2), &
matrix_B(j_start,b_start), size(matrix_B,1),0.d0, & matrix_B(j_start,b_start), size(matrix_B,1),0.d0, &
T(a_start,k_start,b_start), size(T,1)*size(T,2) ) T(a_start,k_start,b_start), size(T,1)*size(T,2) )
deallocate(V) deallocate(V)
allocate( V(a_start:a_end, k_start:k_end, b_start:b_end) )
V = T
deallocate(T)
allocate( T(a_start:a_end, k_start:k_end, b_start:b_end) )
do b=b_start,b_end
do a=a_start,a_end
do k=k_start,k_end
T(a, k,b) = V(a, k,b)
enddo
enddo
enddo
deallocate(V)
allocate( V(a_start:a_end, b_start:b_end, c_start:c_end) ) allocate( V(a_start:a_end, b_start:b_end, c_start:c_end) )
! V = 0.d0
! do b=b_start,b_end
! do c=c_start,c_end
! do a=a_start,a_end
! do k=k_start,k_end
! V(a,b,c) = V(a,b,c) + T(a,k ,b)*matrix_B(k, c)
! enddo
! enddo
! enddo
! enddo
V = 0.d0 V = 0.d0
do b=b_start,b_end do b=b_start,b_end
call DGEMM('N','N', (a_end-a_start+1), (c_end-c_start+1), & call DGEMM('N','N', (a_end-a_start+1), (c_end-c_start+1), &
@ -147,22 +114,8 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB, &
V(a_start,c_start,b), size(V,1) ) V(a_start,c_start,b), size(V,1) )
enddo enddo
deallocate(T)
U = U + V*matrix_B(l, d) U = U + V*matrix_B(l, d)
deallocate(T,V)
! do a=a_start,a_end
! do b=b_start,b_end
! do c=c_start,c_end
! do c=c_start,c_end
! do b=b_start,d
! do a=a_start,min(b,c)
! U(a,b,c) = U(a,b,c) + V(a,b,c) * matrix_B(l, d)
! enddo
! enddo
! enddo
deallocate(V)
enddo enddo
@ -182,11 +135,16 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB, &
enddo enddo
enddo enddo
enddo enddo
!$OMP CRITICAL
call map_append(map_c, key, value, idx) call map_append(map_c, key, value, idx)
call map_sort(map_c) call map_sort(map_c)
call map_unique(map_c) call map_unique(map_c)
!$OMP END CRITICAL
enddo enddo
!$OMP END DO
deallocate(key,value) deallocate(key,value)
!$OMP END PARALLEL
end end