[Fortran] 纯文本查看 复制代码
!$OMP PARALLEL DO PRIVATE(i,j,l,m,pf,th,ll,sr,rl,ls,rs,rg,ag,ipiv,work,lwork,info,fd,tr)&
!$OMP SHARED(po,km,kp,nm,nb,ei,ic,tv,lr,ol,ks,ra,k,ch,te) COLLAPSE(2)&
!$OMP REDUCTION(+:t1,t2,t3)
DO i = 1, nm, 1
DO j = 1, km(1)*km(2)*km(3), 1
DO l = 1, nb, 1
DO m = 1, po(1)*po(2)*po(3), 1
pf = EXP(ic*(kp(j,1)*(tv(1,m)*lr(1,1)+tv(2,m)*lr(2,1)+tv(3,m)*lr(3,1))+kp(j,2)*(tv(1,m)*lr(1,2)+tv(2,m)*lr(2,2)+tv(3,m)*lr(3,2))+kp(j,3)*(tv(1,m)*lr(1,3)+tv(2,m)*lr(2,3)+tv(3,m)*lr(3,3))))
th(:,l) = pf*((ei(i)+ic)*ol(:,l,m)-ks(:,l,m))
END DO
END DO
ll(1,:,:) = th(ra(1)-k:ra(2)-k,ra(1)-k*2:ra(2)-k*2)
ll(2,:,:) = th(ra(1)-k:ra(2)-k,ra(1)-k:ra(2)-k)
ll(3,:,:) = th(ra(1)-k*2:ra(2)-k*2,ra(1)-k:ra(2)-k)
ll(4,:,:) = th(ra(1)-k:ra(2)-k,ra(1):ra(2))
ll(5,:,:) = th(ra(1):ra(2),ra(1)-k:ra(2)-k)
sr = th(ra(1):ra(2),ra(1):ra(2))
rl(1,:,:) = th(ra(1)+k:ra(2)+k,ra(1)+k*2:ra(2)+k*2)
rl(2,:,:) = th(ra(1)+k:ra(2)+k,ra(1)+k:ra(2)+k)
rl(3,:,:) = th(ra(1)+k*2:ra(2)+k*2,ra(1)+k:ra(2)+k)
rl(4,:,:) = th(ra(1):ra(2),ra(1)+k:ra(2)+k)
rl(5,:,:) = th(ra(1)+k:ra(2)+k,ra(1):ra(2))
ls = MATMUL(MATMUL(ll(1,:,:),ll(2,:,:)),ll(3,:,:))
ll(2,:,:) = ll(2,:,:) - ls
rs = MATMUL(MATMUL(rl(1,:,:),rl(2,:,:)),rl(3,:,:))
rl(2,:,:) = rl(2,:,:) - rs
rg = DCMPLX(0.0d0, 0.0d0)
rg = sr-ls-rs
CALL ZGETRF(ra(2)-ra(1),ra(2)-ra(1),rg,l,ipiv,info)
CALL ZGETRI(ra(2)-ra(1),rg,ra(2)-ra(1),ipiv,work,lwork,info)
DO l = 1, k, 1
ag(l,:) = DCONJG(rg(:,l))
END DO
fd(1) = 1.0d0/(EXP((ei(i)-ch(1))/te)+1.0d0)
fd(2) = 1.0d0/(EXP((ei(i)-ch(3))/te)+1.0d0)
IF (fd(1) > fd(2)) THEN
fd(3) = fd(1)
fd(1) = fd(2)
fd(2) = fd(3)
fd(3) = fd(2)-fd(1)
ELSE IF (fd(1) == fd(2)) THEN
fd(3) = 1.0d0
ELSE
fd(3) = fd(2)-fd(1)
END IF
t1(i,j,:,:) = t1(i,j,:,:)+fd(3)*MATMUL(MATMUL(rg,ls+rs),ag)/2.0d0-fd(1)*IMAG(rg)
tr = MATMUL(MATMUL(MATMUL(ls,ag),rs),rg)
DO l = 1, k, 1
t2(i) = t2(i)+tr(l,l)
END DO
t2(i) = t2(i)*fd(3)
tr = rg - ag
DO l = 1, k, 1
t3(i) = t3(i)+tr(l,l)
END DO
END DO
END DO
!$OMP END PARALLEL DO
[Fortran] 纯文本查看 复制代码
PROGRAM TES
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14)
INTEGER :: i, j, k, l, m
INTEGER :: na = 177
INTEGER :: nm = 100
INTEGER :: km(3)
INTEGER :: po(3)
REAL (KIND=dp), ALLOCATABLE :: kp(:,:)
COMPLEX (KIND=dp), ALLOCATABLE :: th(:,:)
COMPLEX (KIND=dp) :: ic = DCMPLX(0.0d0, 1.0d0)
COMPLEX (KIND=dp), ALLOCATABLE :: t1(:,:,:,:)
COMPLEX (KIND=dp), ALLOCATABLE :: t2(:)
COMPLEX (KIND=dp), ALLOCATABLE :: t3(:)
INTEGER :: omp_get_num_threads
INTEGER :: omp_get_thread_num
INTEGER :: num_thread
INTEGER :: id_thread
km(1) = 100
km(2) = 1
km(3) = 75
ALLOCATE (kp(km(1)*km(2)*km(3),3))
DO i = 1, km(1), 1
DO j = 1, km(2), 1
DO k = 1, km(3), 1
kp((i-1)*km(2)*km(3)+(j-1)*km(3)+k,1) = DBLE(i)
kp((i-1)*km(2)*km(3)+(j-1)*km(3)+k,2) = DBLE(j)
kp((i-1)*km(2)*km(3)+(j-1)*km(3)+k,3) = DBLE(k)
END DO
END DO
END DO
po(1) = 1
po(2) = 0
po(3) = 5
ALLOCATE (th(na,na))
ALLOCATE (t1(nm,km(1)*km(2)*km(3),na,na))
ALLOCATE (t2(nm))
ALLOCATE (t3(nm))
!$OMP PARALLEL DEFAULT(NONE) SHARED(num_thread)
num_thread = omp_get_num_threads()
!$OMP END PARALLEL
!$OMP PARALLEL DO PRIVATE(i,j,k,l,m,th)&
!$OMP SHARED(km,kp,nm,po,ic) COLLAPSE(2)&
!$OMP REDUCTION(+:t1,t2,t3)
DO i = 1, nm, 1
DO j = 1, km(1)*km(2)*km(3), 1
th = DCMPLX(0.0d0, 0.0d0)
DO k = 1, na, 1
DO l = 1, na, 1
DO m = 1, (po(1)*2+1)*(po(2)*2+1)*(po(3)*2+1), 1
th(l,k) = th(l,k) + kp(j,1)*k + kp(j,2)*ic - kp(j,3)*l
END DO
END DO
END DO
t1(i,j,:,:) = th*(i-j+ic)
DO k = 1, na, 1
t2(i) = t2(i) + th(k,k)*(i*ic-j)
END DO
DO k = na, 1, -1
t3(i) = t3(i) + th(k,na-k+1)*(i*ic+j)
END DO
END DO
END DO
!$OMP END PARALLEL DO
OPEN (UNIT=3, FILE='output.dat', STATUS='UNKNOWN')
DO i = 1, nm, 1
DO j = 1, km(1)*km(2)*km(3), 1
DO k = 1, na, 1
WRITE (UNIT=3, FMT=*) t1(i,j,k,:)
END DO
END DO
END DO
WRITE (UNIT=3, FMT=*)
DO i = 1, nm, 1
WRITE (UNIT=3, FMT=*) t2(i)
END DO
WRITE (UNIT=3, FMT=*)
DO i = 1, nm, 1
WRITE (UNIT=3, FMT=*) t3(i)
END DO
CLOSE (UNIT=3)
DEALLOCATE (kp)
DEALLOCATE (th)
DEALLOCATE (t1)
DEALLOCATE (t2)
DEALLOCATE (t3)
STOP
END PROGRAM TES
[Fortran] 纯文本查看 复制代码
SUBROUTINE MULTIPLICATION_MATRIX(di,ou,i1,i2)
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14) !Precision of the floating variable
INTEGER :: di !Dimension of the matrix
COMPLEX (KIND=dp), ALLOCATABLE :: ou(:,:) !Output matrix
COMPLEX (KIND=dp), ALLOCATABLE :: i1(:,:) !Input matrix
COMPLEX (KIND=dp), ALLOCATABLE :: i2(:,:) !Input matrix
COMPLEX (KIND=dp) :: i3
INTEGER :: i, j, k
ou = DCMPLX(0.0d0, 0.0d0)
DO i = 1, di, 1
DO j = 1, di, 1
i3 = DCMPLX(0.0d0, 0.0d0)
DO k = 1, di, 1
! ou(i,j) = ou(i,j)+i1(j,k)*i2(k,j)
i3 = i3+i1(j,k)*i2(k,j)
END DO
ou(i,j) = i3
END DO
END DO
RETURN
END SUBROUTINE MULTIPLICATION_MATRIX
PROGRAM CODE
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14) !Precision of the floating variable
INTEGER :: di = 2 !Dimension of the matrix
COMPLEX (KIND=dp), ALLOCATABLE :: ou(:,:) !Output matrix
COMPLEX (KIND=dp), ALLOCATABLE :: i1(:,:) !Input matrix
COMPLEX (KIND=dp), ALLOCATABLE :: i2(:,:) !Input matrix
INTEGER :: i, j
ALLOCATE (ou(di,di))
ALLOCATE (i1(di,di))
ALLOCATE (i2(di,di))
DO i = 1, di, 1
DO j = 1, di, 1
i1(i,j) = DCMPLX(0.0d0, 1.0d0)*i-j
i2(i,j) = DCMPLX(1.0d0, 1.0d0)*j-i
END DO
END DO
CALL MULTIPLICATION_MATRIX(di,ou,i1,i2)
OPEN (UNIT=3, FILE='out.dat', STATUS='UNKNOWN')
DO i = 1, di, 1
WRITE (UNIT=3, FMT=*) ou(i,:)
END DO
DEALLOCATE (ou)
DEALLOCATE (i1)
DEALLOCATE (i2)
CLOSE (UNIT=3)
STOP
END PROGRAM CODE