diff --git a/Documentation/fassaad.bib b/Documentation/fassaad.bib index ac06a5618a715293489acccac8f6ebe18cebb9ca..1d8484d818f1b82b25615dab4f882e9dc0ad6b7b 100644 --- a/Documentation/fassaad.bib +++ b/Documentation/fassaad.bib @@ -3260,7 +3260,7 @@ Title = {Statistical Mechanics: Algorithms and Computations}, Year = {2006}} -@article{neal1993, +@book{neal1993, Author = {Neal, Radford M}, Date-Added = {2018-05-14 18:37:55 +0000}, Date-Modified = {2018-05-14 18:37:55 +0000}, diff --git a/Documentation/interaction_vertices.tex b/Documentation/interaction_vertices.tex index b330184752dccec7ea2bf6decad056bb460f49cd..7ed97576724d91079f2fbb8abcd5f940bb56fd01 100644 --- a/Documentation/interaction_vertices.tex +++ b/Documentation/interaction_vertices.tex @@ -181,6 +181,7 @@ Another predefined vertex is: - \frac{|J_z|}{2} \left( S^{z}_i - \sgn|J_z| S^{z}_j \right)^2 = J_z S^{z}_i S^{z}_j - \frac{|J_z|}{2} (S^{z}_i)^2 - \frac{|J_z|}{2}(S^{z}_j)^2 \end{align} + which, if particle fluctuations are frozen on the $i$ and $j$ sites, then $(S^{z}_i)^2 = 1/4$ and the interaction corresponds to a $J_z$-$J_z$ ferromagnetic or antiferromagnetic coupling. The implementation of the interaction in \texttt{Predefined\_Int\_Jz} defines two operators: diff --git a/Libraries/Modules/Mat_subroutines.F90 b/Libraries/Modules/Mat_subroutines.F90 index 1d2731740e94799651222b94465d8a946e288059..177051ad1fd4573cd7845299122f5b73893a971d 100644 --- a/Libraries/Modules/Mat_subroutines.F90 +++ b/Libraries/Modules/Mat_subroutines.F90 @@ -957,3 +957,464 @@ subroutine ZSLHEMM(side, uplo, N, M1, M2, A, P, Mat) ENDIF end subroutine ZSLHEMM + + +subroutine ZDSLSYMM(side, uplo, N, M1, M2, A, P, Mat) +! Mixed Complex-Real Small Large symmetric matrix multiplication + +!-------------------------------------------------------------------- +!> @author +!> ALF-project +! +!> @brief +!> P^T A P is symmetric +!> !!!!! Side = L +!> M = op( P^T A P) * M +!> Side = R +!> M = M * op( P^T A P) +!> On input: P = Op%P and A = Op%O +!> !!!!! type +!> op = N --> None +!> op = T --> Transposed +!> op = C --> Transposed + Complex conjugation. Same as N. +!> !!!!! Mat has dimensions M1,M2 +!> +!-------------------------------------------------------------------- + use iso_fortran_env, only: output_unit, error_unit + + IMPLICIT NONE + CHARACTER (1) , INTENT(IN) :: side, uplo + INTEGER , INTENT(IN) :: N, M1, M2 + REAL (KIND=KIND(0.D0)), INTENT(IN) , DIMENSION(N,N) :: A + COMPLEX (KIND=KIND(0.D0)), INTENT(INOUT), DIMENSION(M1,M2) :: Mat + INTEGER , INTENT(IN) , DIMENSION(N) :: P + + REAL (KIND=KIND(0.D0)), DIMENSION(:,:), ALLOCATABLE :: WORK + REAL (KIND=KIND(0.D0)), DIMENSION(:), ALLOCATABLE :: RWORK ! required for the ZLACRM calls + COMPLEX (KIND=KIND(0.D0)), DIMENSION(:,:), ALLOCATABLE :: WORK2, WORKCMPLX + Complex (Kind = Kind(0.D0)) :: alpha, beta, Z(8) + INTEGER :: I,L,IDX, NUMBLOCKS + INTEGER, DIMENSION(:), ALLOCATABLE :: IDXLIST, DIMLIST + LOGICAL :: COMPACT, ALLOC + + alpha = 1.D0 + beta = 0.D0 + + !identify possible block structure + !only used in default case for n>4 + IF(N > 8) THEN + IF(uplo=='L' .or. uplo=='l') THEN ! uplo == l unimplemented and never used in this case + ERROR STOP + ENDIF + COMPACT = .TRUE. + L = 1 + IDX = 1 + ALLOCATE(IDXLIST(N),DIMLIST(N)) + ALLOC = .TRUE. + NUMBLOCKS=0 + DO I=1,N-1 + IF ( P(I)+1 .ne. P(I+1) ) THEN + COMPACT = .FALSE. + NUMBLOCKS=NUMBLOCKS+1 + IDXLIST(NUMBLOCKS)=IDX + DIMLIST(NUMBLOCKS)=L + IDX=IDX+L + L=1 + ELSE + L=L+1 + ENDIF + ENDDO + IF(IDX1) THEN + ALLOCATE(WORK(N,N)) + CALL DLACPY(uplo, N, N, A(1,1), N, WORK(1,1), N) + ! Fill the rest of WORK (thereby ignoring uplo) + IF(uplo=='U' .or. uplo=='u') THEN + DO L=1,N + DO I=1,L-1 + WORK(L,I)=WORK(I,L) + ENDDO + ENDDO + ELSE + DO L=1,N + DO I=L+1,N + WORK(L,I)=WORK(I,L) + ENDDO + ENDDO + ENDIF + ENDIF + + IF ( side == 'L' .or. side == 'l' ) THEN + ! multiply op(A) from the left [ Mat = op(A)*Mat ] + + SELECT CASE(N) + CASE (1) + ! Here only one row is rescaled + ! uplo and transpositions have no effect for 1x1 herm. Matrices + CALL ZDSCAL(M2,A(1,1),Mat(P(1),1),M1) + CASE (2) + ! perform inplace matmult +! L=M2/4 +! DO I=1,4*L-1,4 +! Z(1)=Mat(P(1),I) +! Z(2)=Mat(P(2),I) +! Z(3)=Mat(P(1),I+1) +! Z(4)=Mat(P(2),I+1) +! Z(5)=Mat(P(1),I+2) +! Z(6)=Mat(P(2),I+2) +! Z(7)=Mat(P(1),I+3) +! Z(8)=Mat(P(2),I+3) +! Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2) +! Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2) +! Mat(P(1),I+1)=WORK(1,1)*Z(3)+WORK(1,2)*Z(4) +! Mat(P(2),I+1)=WORK(2,1)*Z(3)+WORK(2,2)*Z(4) +! Mat(P(1),I+2)=WORK(1,1)*Z(5)+WORK(1,2)*Z(6) +! Mat(P(2),I+2)=WORK(2,1)*Z(5)+WORK(2,2)*Z(6) +! Mat(P(1),I+3)=WORK(1,1)*Z(7)+WORK(1,2)*Z(8) +! Mat(P(2),I+3)=WORK(2,1)*Z(7)+WORK(2,2)*Z(8) +! ENDDO + DO I=1,M2!4*L,M2 + Z(1)=Mat(P(1),I) + Z(2)=Mat(P(2),I) + Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2) + Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2) + ENDDO + DEALLOCATE(WORK) + CASE (3) + ! perform inplace matmult +! L=M2/2 +! DO I=1,2*L-1,2 +! Z(1)=Mat(P(1),I) +! Z(2)=Mat(P(2),I) +! Z(3)=Mat(P(3),I) +! Z(4)=Mat(P(1),I+1) +! Z(5)=Mat(P(2),I+1) +! Z(6)=Mat(P(3),I+1) +! Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2)+WORK(1,3)*Z(3) +! Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2)+WORK(2,3)*Z(3) +! Mat(P(3),I)=WORK(3,1)*Z(1)+WORK(3,2)*Z(2)+WORK(3,3)*Z(3) +! Mat(P(1),I+1)=WORK(1,1)*Z(4)+WORK(1,2)*Z(5)+WORK(1,3)*Z(6) +! Mat(P(2),I+1)=WORK(2,1)*Z(4)+WORK(2,2)*Z(5)+WORK(2,3)*Z(6) +! Mat(P(3),I+1)=WORK(3,1)*Z(4)+WORK(3,2)*Z(5)+WORK(3,3)*Z(6) +! ENDDO + DO I=1,M2!2*L,M2 + Z(1)=Mat(P(1),I) + Z(2)=Mat(P(2),I) + Z(3)=Mat(P(3),I) + Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2)+WORK(1,3)*Z(3) + Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2)+WORK(2,3)*Z(3) + Mat(P(3),I)=WORK(3,1)*Z(1)+WORK(3,2)*Z(2)+WORK(3,3)*Z(3) + ENDDO + DEALLOCATE(WORK) + CASE (4) + ! perform inplace matmult + DO I=1,M2 + Z(1)=Mat(P(1),I) + Z(2)=Mat(P(2),I) + Z(3)=Mat(P(3),I) + Z(4)=Mat(P(4),I) + Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2)+WORK(1,3)*Z(3)+WORK(1,4)*Z(4) + Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2)+WORK(2,3)*Z(3)+WORK(2,4)*Z(4) + Mat(P(3),I)=WORK(3,1)*Z(1)+WORK(3,2)*Z(2)+WORK(3,3)*Z(3)+WORK(3,4)*Z(4) + Mat(P(4),I)=WORK(4,1)*Z(1)+WORK(4,2)*Z(2)+WORK(4,3)*Z(3)+WORK(4,4)*Z(4) + ENDDO + DEALLOCATE(WORK) + CASE (5) + ! perform inplace matmult + DO I=1,M2 + Z(1)=Mat(P(1),I) + Z(2)=Mat(P(2),I) + Z(3)=Mat(P(3),I) + Z(4)=Mat(P(4),I) + Z(5)=Mat(P(5),I) + Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2)+WORK(1,3)*Z(3)+WORK(1,4)*Z(4)+WORK(1,5)*Z(5) + Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2)+WORK(2,3)*Z(3)+WORK(2,4)*Z(4)+WORK(2,5)*Z(5) + Mat(P(3),I)=WORK(3,1)*Z(1)+WORK(3,2)*Z(2)+WORK(3,3)*Z(3)+WORK(3,4)*Z(4)+WORK(3,5)*Z(5) + Mat(P(4),I)=WORK(4,1)*Z(1)+WORK(4,2)*Z(2)+WORK(4,3)*Z(3)+WORK(4,4)*Z(4)+WORK(4,5)*Z(5) + Mat(P(5),I)=WORK(5,1)*Z(1)+WORK(5,2)*Z(2)+WORK(5,3)*Z(3)+WORK(5,4)*Z(4)+WORK(5,5)*Z(5) + ENDDO + DEALLOCATE(WORK) + CASE (6) + ! perform inplace matmult + DO I=1,M2 + Z(1)=Mat(P(1),I) + Z(2)=Mat(P(2),I) + Z(3)=Mat(P(3),I) + Z(4)=Mat(P(4),I) + Z(5)=Mat(P(5),I) + Z(6)=Mat(P(6),I) + Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2)+WORK(1,3)*Z(3)& + &+WORK(1,4)*Z(4)+WORK(1,5)*Z(5)+WORK(1,6)*Z(6) + Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2)+WORK(2,3)*Z(3)& + &+WORK(2,4)*Z(4)+WORK(2,5)*Z(5)+WORK(2,6)*Z(6) + Mat(P(3),I)=WORK(3,1)*Z(1)+WORK(3,2)*Z(2)+WORK(3,3)*Z(3)& + &+WORK(3,4)*Z(4)+WORK(3,5)*Z(5)+WORK(3,6)*Z(6) + Mat(P(4),I)=WORK(4,1)*Z(1)+WORK(4,2)*Z(2)+WORK(4,3)*Z(3)& + &+WORK(4,4)*Z(4)+WORK(4,5)*Z(5)+WORK(4,6)*Z(6) + Mat(P(5),I)=WORK(5,1)*Z(1)+WORK(5,2)*Z(2)+WORK(5,3)*Z(3)& + &+WORK(5,4)*Z(4)+WORK(5,5)*Z(5)+WORK(5,6)*Z(6) + Mat(P(6),I)=WORK(6,1)*Z(1)+WORK(6,2)*Z(2)+WORK(6,3)*Z(3)& + &+WORK(6,4)*Z(4)+WORK(6,5)*Z(5)+WORK(6,6)*Z(6) + ENDDO + DEALLOCATE(WORK) + CASE (7) + ! perform inplace matmult + DO I=1,M2 + Z(1)=Mat(P(1),I) + Z(2)=Mat(P(2),I) + Z(3)=Mat(P(3),I) + Z(4)=Mat(P(4),I) + Z(5)=Mat(P(5),I) + Z(6)=Mat(P(6),I) + Z(7)=Mat(P(7),I) + Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2)+WORK(1,3)*Z(3)& + &+WORK(1,4)*Z(4)+WORK(1,5)*Z(5)+WORK(1,6)*Z(6)+WORK(1,7)*Z(7) + Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2)+WORK(2,3)*Z(3)& + &+WORK(2,4)*Z(4)+WORK(2,5)*Z(5)+WORK(2,6)*Z(6)+WORK(2,7)*Z(7) + Mat(P(3),I)=WORK(3,1)*Z(1)+WORK(3,2)*Z(2)+WORK(3,3)*Z(3)& + &+WORK(3,4)*Z(4)+WORK(3,5)*Z(5)+WORK(3,6)*Z(6)+WORK(3,7)*Z(7) + Mat(P(4),I)=WORK(4,1)*Z(1)+WORK(4,2)*Z(2)+WORK(4,3)*Z(3)& + &+WORK(4,4)*Z(4)+WORK(4,5)*Z(5)+WORK(4,6)*Z(6)+WORK(4,7)*Z(7) + Mat(P(5),I)=WORK(5,1)*Z(1)+WORK(5,2)*Z(2)+WORK(5,3)*Z(3)& + &+WORK(5,4)*Z(4)+WORK(5,5)*Z(5)+WORK(5,6)*Z(6)+WORK(5,7)*Z(7) + Mat(P(6),I)=WORK(6,1)*Z(1)+WORK(6,2)*Z(2)+WORK(6,3)*Z(3)& + &+WORK(6,4)*Z(4)+WORK(6,5)*Z(5)+WORK(6,6)*Z(6)+WORK(6,7)*Z(7) + Mat(P(7),I)=WORK(7,1)*Z(1)+WORK(7,2)*Z(2)+WORK(7,3)*Z(3)& + &+WORK(7,4)*Z(4)+WORK(7,5)*Z(5)+WORK(7,6)*Z(6)+WORK(7,7)*Z(7) + ENDDO + DEALLOCATE(WORK) + CASE (8) + ! perform inplace matmult + DO I=1,M2 + Z(1)=Mat(P(1),I) + Z(2)=Mat(P(2),I) + Z(3)=Mat(P(3),I) + Z(4)=Mat(P(4),I) + Z(5)=Mat(P(5),I) + Z(6)=Mat(P(6),I) + Z(7)=Mat(P(7),I) + Z(8)=Mat(P(8),I) + Mat(P(1),I)=WORK(1,1)*Z(1)+WORK(1,2)*Z(2)+WORK(1,3)*Z(3)+WORK(1,4)*Z(4)& + &+WORK(1,5)*Z(5)+WORK(1,6)*Z(6)+WORK(1,7)*Z(7)+WORK(1,8)*Z(8) + Mat(P(2),I)=WORK(2,1)*Z(1)+WORK(2,2)*Z(2)+WORK(2,3)*Z(3)+WORK(2,4)*Z(4)& + &+WORK(2,5)*Z(5)+WORK(2,6)*Z(6)+WORK(2,7)*Z(7)+WORK(2,8)*Z(8) + Mat(P(3),I)=WORK(3,1)*Z(1)+WORK(3,2)*Z(2)+WORK(3,3)*Z(3)+WORK(3,4)*Z(4)& + &+WORK(3,5)*Z(5)+WORK(3,6)*Z(6)+WORK(3,7)*Z(7)+WORK(3,8)*Z(8) + Mat(P(4),I)=WORK(4,1)*Z(1)+WORK(4,2)*Z(2)+WORK(4,3)*Z(3)+WORK(4,4)*Z(4)& + &+WORK(4,5)*Z(5)+WORK(4,6)*Z(6)+WORK(4,7)*Z(7)+WORK(4,8)*Z(8) + Mat(P(5),I)=WORK(5,1)*Z(1)+WORK(5,2)*Z(2)+WORK(5,3)*Z(3)+WORK(5,4)*Z(4)& + &+WORK(5,5)*Z(5)+WORK(5,6)*Z(6)+WORK(5,7)*Z(7)+WORK(5,8)*Z(8) + Mat(P(6),I)=WORK(6,1)*Z(1)+WORK(6,2)*Z(2)+WORK(6,3)*Z(3)+WORK(6,4)*Z(4)& + &+WORK(6,5)*Z(5)+WORK(6,6)*Z(6)+WORK(6,7)*Z(7)+WORK(6,8)*Z(8) + Mat(P(7),I)=WORK(7,1)*Z(1)+WORK(7,2)*Z(2)+WORK(7,3)*Z(3)+WORK(7,4)*Z(4)& + &+WORK(7,5)*Z(5)+WORK(7,6)*Z(6)+WORK(7,7)*Z(7)+WORK(7,8)*Z(8) + Mat(P(8),I)=WORK(8,1)*Z(1)+WORK(8,2)*Z(2)+WORK(8,3)*Z(3)+WORK(8,4)*Z(4)& + &+WORK(8,5)*Z(5)+WORK(8,6)*Z(6)+WORK(8,7)*Z(7)+WORK(8,8)*Z(8) + ENDDO + DEALLOCATE(WORK) + CASE DEFAULT + ! allocate memory and copy blocks of Mat to work + ALLOCATE(WORKCMPLX(N,M2)) + DO I=1,NUMBLOCKS + CALL ZLACPY('A', DIMLIST(I), M2, Mat(P(IDXLIST(I)),1), M1, WORKCMPLX(IDXLIST(I),1), N) + ENDDO + + ! Perform Mat multiplication + IF(COMPACT) THEN + !write result directly into mat + allocate (RWORK(2*N*M2)) + call ZLARCM(N, M2, A(1, 1), N, WORKCMPLX(1,1), N, Mat(P(1), 1), M1, RWORK) + +! CALL ZHEMM(side, uplo, N, M2, alpha, A(1, 1), N, WORKCMPLX(1, 1), N, beta, Mat(P(1), 1), M1) + DEALLOCATE(RWORK) + ELSE + !additional space for result + ALLOCATE(WORK2(N,M2), RWORK(2*N*M2)) + call ZLARCM(N, M2, A(1, 1), N, WORKCMPLX(1,1), N, WORK2(1,1), N, RWORK) +! CALL ZHEMM(side, uplo, N, M2, alpha, A(1, 1), N, WORKCMPLX(1, 1), N, beta, WORK2(1, 1), N) + !distribute result back into mat using blocks + DO I=1,NUMBLOCKS + CALL ZLACPY('A', DIMLIST(I), M2, WORK2(IDXLIST(I),1), N, Mat(P(IDXLIST(I)),1), M1) + ENDDO + !free result memory + DEALLOCATE(WORK2, RWORK) + ENDIF + !free memory of first mat copy + DEALLOCATE(WORKCMPLX,IDXLIST,DIMLIST) + END SELECT + + ELSEIF ( side == 'R' .or. side == 'r' ) THEN + ! multiply op(A) from the right [ Mat = Mat*op(A) ] + + SELECT CASE(N) + CASE (1) + ! Here only one column is rescaled + ! uplo and transpositions have no effect for 1x1 matrices + CALL ZDSCAL(M1,A(1,1),Mat(1,P(1)),1) + CASE (2) + ! perform inplace matmult + DO I=1,M1 + Z(1)=Mat(I,P(1)) + Z(2)=Mat(I,P(2)) + Mat(I,P(1))=WORK(1,1)*Z(1)+WORK(2,1)*Z(2) + Mat(I,P(2))=WORK(1,2)*Z(1)+WORK(2,2)*Z(2) + ENDDO + DEALLOCATE(WORK) + CASE (3) + ! perform inplace matmult + DO I=1,M1 + Z(1)=Mat(I,P(1)) + Z(2)=Mat(I,P(2)) + Z(3)=Mat(I,P(3)) + Mat(I,P(1))=WORK(1,1)*Z(1)+WORK(2,1)*Z(2)+WORK(3,1)*Z(3) + Mat(I,P(2))=WORK(1,2)*Z(1)+WORK(2,2)*Z(2)+WORK(3,2)*Z(3) + Mat(I,P(3))=WORK(1,3)*Z(1)+WORK(2,3)*Z(2)+WORK(3,3)*Z(3) + ENDDO + DEALLOCATE(WORK) + CASE (4) + ! perform inplace matmult + DO I=1,M1 + Z(1)=Mat(I,P(1)) + Z(2)=Mat(I,P(2)) + Z(3)=Mat(I,P(3)) + Z(4)=Mat(I,P(4)) + Mat(I,P(1))=WORK(1,1)*Z(1)+WORK(2,1)*Z(2)+WORK(3,1)*Z(3)+WORK(4,1)*Z(4) + Mat(I,P(2))=WORK(1,2)*Z(1)+WORK(2,2)*Z(2)+WORK(3,2)*Z(3)+WORK(4,2)*Z(4) + Mat(I,P(3))=WORK(1,3)*Z(1)+WORK(2,3)*Z(2)+WORK(3,3)*Z(3)+WORK(4,3)*Z(4) + Mat(I,P(4))=WORK(1,4)*Z(1)+WORK(2,4)*Z(2)+WORK(3,4)*Z(3)+WORK(4,4)*Z(4) + ENDDO + DEALLOCATE(WORK) + CASE (5) + ! perform inplace matmult + DO I=1,M1 + Z(1)=Mat(I,P(1)) + Z(2)=Mat(I,P(2)) + Z(3)=Mat(I,P(3)) + Z(4)=Mat(I,P(4)) + Z(5)=Mat(I,P(5)) + Mat(I,P(1))=WORK(1,1)*Z(1)+WORK(2,1)*Z(2)+WORK(3,1)*Z(3)+WORK(4,1)*Z(4)+WORK(5,1)*Z(5) + Mat(I,P(2))=WORK(1,2)*Z(1)+WORK(2,2)*Z(2)+WORK(3,2)*Z(3)+WORK(4,2)*Z(4)+WORK(5,2)*Z(5) + Mat(I,P(3))=WORK(1,3)*Z(1)+WORK(2,3)*Z(2)+WORK(3,3)*Z(3)+WORK(4,3)*Z(4)+WORK(5,3)*Z(5) + Mat(I,P(4))=WORK(1,4)*Z(1)+WORK(2,4)*Z(2)+WORK(3,4)*Z(3)+WORK(4,4)*Z(4)+WORK(5,4)*Z(5) + Mat(I,P(5))=WORK(1,5)*Z(1)+WORK(2,5)*Z(2)+WORK(3,5)*Z(3)+WORK(4,5)*Z(4)+WORK(5,5)*Z(5) + ENDDO + DEALLOCATE(WORK) + CASE (6) + ! perform inplace matmult + DO I=1,M1 + Z(1)=Mat(I,P(1)) + Z(2)=Mat(I,P(2)) + Z(3)=Mat(I,P(3)) + Z(4)=Mat(I,P(4)) + Z(5)=Mat(I,P(5)) + Z(6)=Mat(I,P(6)) + Mat(I,P(1))=WORK(1,1)*Z(1)+WORK(2,1)*Z(2)+WORK(3,1)*Z(3)& + &+WORK(4,1)*Z(4)+WORK(5,1)*Z(5)+WORK(6,1)*Z(6) + Mat(I,P(2))=WORK(1,2)*Z(1)+WORK(2,2)*Z(2)+WORK(3,2)*Z(3)& + &+WORK(4,2)*Z(4)+WORK(5,2)*Z(5)+WORK(6,2)*Z(6) + Mat(I,P(3))=WORK(1,3)*Z(1)+WORK(2,3)*Z(2)+WORK(3,3)*Z(3)& + &+WORK(4,3)*Z(4)+WORK(5,3)*Z(5)+WORK(6,3)*Z(6) + Mat(I,P(4))=WORK(1,4)*Z(1)+WORK(2,4)*Z(2)+WORK(3,4)*Z(3)& + &+WORK(4,4)*Z(4)+WORK(5,4)*Z(5)+WORK(6,4)*Z(6) + Mat(I,P(5))=WORK(1,5)*Z(1)+WORK(2,5)*Z(2)+WORK(3,5)*Z(3)& + &+WORK(4,5)*Z(4)+WORK(5,5)*Z(5)+WORK(6,5)*Z(6) + Mat(I,P(6))=WORK(1,6)*Z(1)+WORK(2,6)*Z(2)+WORK(3,6)*Z(3)& + &+WORK(4,6)*Z(4)+WORK(5,6)*Z(5)+WORK(6,6)*Z(6) + ENDDO + DEALLOCATE(WORK) + CASE (7) + ! perform inplace matmult + DO I=1,M1 + Z(1)=Mat(I,P(1)) + Z(2)=Mat(I,P(2)) + Z(3)=Mat(I,P(3)) + Z(4)=Mat(I,P(4)) + Z(5)=Mat(I,P(5)) + Z(6)=Mat(I,P(6)) + Z(7)=Mat(I,P(7)) + Mat(I,P(1))=WORK(1,1)*Z(1)+WORK(2,1)*Z(2)+WORK(3,1)*Z(3)& + &+WORK(4,1)*Z(4)+WORK(5,1)*Z(5)+WORK(6,1)*Z(6)+WORK(7,1)*Z(7) + Mat(I,P(2))=WORK(1,2)*Z(1)+WORK(2,2)*Z(2)+WORK(3,2)*Z(3)& + &+WORK(4,2)*Z(4)+WORK(5,2)*Z(5)+WORK(6,2)*Z(6)+WORK(7,2)*Z(7) + Mat(I,P(3))=WORK(1,3)*Z(1)+WORK(2,3)*Z(2)+WORK(3,3)*Z(3)& + &+WORK(4,3)*Z(4)+WORK(5,3)*Z(5)+WORK(6,3)*Z(6)+WORK(7,3)*Z(7) + Mat(I,P(4))=WORK(1,4)*Z(1)+WORK(2,4)*Z(2)+WORK(3,4)*Z(3)& + &+WORK(4,4)*Z(4)+WORK(5,4)*Z(5)+WORK(6,4)*Z(6)+WORK(7,4)*Z(7) + Mat(I,P(5))=WORK(1,5)*Z(1)+WORK(2,5)*Z(2)+WORK(3,5)*Z(3)& + &+WORK(4,5)*Z(4)+WORK(5,5)*Z(5)+WORK(6,5)*Z(6)+WORK(7,5)*Z(7) + Mat(I,P(6))=WORK(1,6)*Z(1)+WORK(2,6)*Z(2)+WORK(3,6)*Z(3)& + &+WORK(4,6)*Z(4)+WORK(5,6)*Z(5)+WORK(6,6)*Z(6)+WORK(7,6)*Z(7) + Mat(I,P(7))=WORK(1,7)*Z(1)+WORK(2,7)*Z(2)+WORK(3,7)*Z(3)& + &+WORK(4,7)*Z(4)+WORK(5,7)*Z(5)+WORK(6,7)*Z(6)+WORK(7,7)*Z(7) + ENDDO + DEALLOCATE(WORK) + CASE (8) + ! perform inplace matmult + DO I=1,M1 + Z(1)=Mat(I,P(1)) + Z(2)=Mat(I,P(2)) + Z(3)=Mat(I,P(3)) + Z(4)=Mat(I,P(4)) + Z(5)=Mat(I,P(5)) + Z(6)=Mat(I,P(6)) + Z(7)=Mat(I,P(7)) + Z(8)=Mat(I,P(8)) + Mat(I,P(1))=WORK(1,1)*Z(1)+WORK(2,1)*Z(2)+WORK(3,1)*Z(3)+WORK(4,1)*Z(4)& + &+WORK(5,1)*Z(5)+WORK(6,1)*Z(6)+WORK(7,1)*Z(7)+WORK(8,1)*Z(8) + Mat(I,P(2))=WORK(1,2)*Z(1)+WORK(2,2)*Z(2)+WORK(3,2)*Z(3)+WORK(4,2)*Z(4)& + &+WORK(5,2)*Z(5)+WORK(6,2)*Z(6)+WORK(7,2)*Z(7)+WORK(8,2)*Z(8) + Mat(I,P(3))=WORK(1,3)*Z(1)+WORK(2,3)*Z(2)+WORK(3,3)*Z(3)+WORK(4,3)*Z(4)& + &+WORK(5,3)*Z(5)+WORK(6,3)*Z(6)+WORK(7,3)*Z(7)+WORK(8,3)*Z(8) + Mat(I,P(4))=WORK(1,4)*Z(1)+WORK(2,4)*Z(2)+WORK(3,4)*Z(3)+WORK(4,4)*Z(4)& + &+WORK(5,4)*Z(5)+WORK(6,4)*Z(6)+WORK(7,4)*Z(7)+WORK(8,4)*Z(8) + Mat(I,P(5))=WORK(1,5)*Z(1)+WORK(2,5)*Z(2)+WORK(3,5)*Z(3)+WORK(4,5)*Z(4)& + &+WORK(5,5)*Z(5)+WORK(6,5)*Z(6)+WORK(7,5)*Z(7)+WORK(8,5)*Z(8) + Mat(I,P(6))=WORK(1,6)*Z(1)+WORK(2,6)*Z(2)+WORK(3,6)*Z(3)+WORK(4,6)*Z(4)& + &+WORK(5,6)*Z(5)+WORK(6,6)*Z(6)+WORK(7,6)*Z(7)+WORK(8,6)*Z(8) + Mat(I,P(7))=WORK(1,7)*Z(1)+WORK(2,7)*Z(2)+WORK(3,7)*Z(3)+WORK(4,7)*Z(4)& + &+WORK(5,7)*Z(5)+WORK(6,7)*Z(6)+WORK(7,7)*Z(7)+WORK(8,7)*Z(8) + Mat(I,P(8))=WORK(1,8)*Z(1)+WORK(2,8)*Z(2)+WORK(3,8)*Z(3)+WORK(4,8)*Z(4)& + &+WORK(5,8)*Z(5)+WORK(6,8)*Z(6)+WORK(7,8)*Z(7)+WORK(8,8)*Z(8) + ENDDO + DEALLOCATE(WORK) + CASE DEFAULT + ! allocate memory and copy blocks of Mat to work + ALLOCATE(WORKCMPLX(M1,N)) + DO I=1,NUMBLOCKS + CALL ZLACPY('A', M1, DIMLIST(I), Mat(1,P(IDXLIST(I))), M1, WORKCMPLX(1,IDXLIST(I)), M1) + ENDDO + + ! Perform Mat multiplication + IF(COMPACT) THEN ! B*A + C + ALLOCATE(RWORK(2*M1*N)) + !write result directly into mat + call ZLACRM(M1, N, WORKCMPLX(1,1), M1, A(1, 1), N, Mat(1, P(1)), M1, RWORK) +! CALL ZHEMM(side, uplo, M1, N, alpha, A(1, 1), N, WORKCMPLX(1, 1), M1, beta, Mat(1, P(1)), M1) + DEALLOCATE(RWORK) + ELSE + !additional space for result + ALLOCATE(WORK2(M1,N), RWORK(2*M1*N)) + call ZLACRM(M1, N, WORKCMPLX(1,1), M1, A(1, 1), N, WORK2(1, 1), M1, RWORK) +! CALL ZHEMM(side, uplo, M1, N, alpha, A(1, 1), N, WORKCMPLX(1, 1), M1, beta, WORK2(1, 1), M1) + !distribute result back into mat using blocks + DO I=1,NUMBLOCKS + CALL ZLACPY('A', M1, DIMLIST(I), WORK2(1,IDXLIST(I)), M1, Mat(1,P(IDXLIST(I))), M1) + ENDDO + !free result memory + DEALLOCATE(WORK2, RWORK) + ENDIF + !free memory of first mat copy + DEALLOCATE(WORKCMPLX,IDXLIST,DIMLIST) + END SELECT + + ELSE + write(error_unit,*) 'ZDSLSYMM: Illegal argument for side: It is not one of [R,r,L,l] !' + error stop 1 + ENDIF + +end subroutine ZDSLSYMM diff --git a/Prog/ContainerElementBase_mod.F90 b/Prog/ContainerElementBase_mod.F90 new file mode 100644 index 0000000000000000000000000000000000000000..decd5b21183e90cc62c3dfcfe80a11fc1c176d82 --- /dev/null +++ b/Prog/ContainerElementBase_mod.F90 @@ -0,0 +1,138 @@ +! Copyright (C) 2020 The ALF project +! +! The ALF project is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! The ALF project is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with ALF. If not, see http://www.gnu.org/licenses/. +! +! Under Section 7 of GPL version 3 we require you to fulfill the following additional terms: +! +! - It is our hope that this program makes a contribution to the scientific community. Being +! part of that community we feel that it is reasonable to require you to give an attribution +! back to the original authors if you have benefitted from this program. +! Guidelines for a proper citation can be found on the project's homepage +! https://alf.physik.uni-wuerzburg.de . +! +! - We require the preservation of the above copyright notice and this license in all original files. +! +! - We prohibit the misrepresentation of the origin of the original source files. To obtain +! the original source files please visit the homepage https://alf.physik.uni-wuerzburg.de . +! +! - If you make substantial changes to the program we require you to either consider contributing +! to the ALF project or to mark your material in a reasonable way as different from the original version. + + +! Declare a common base class for the interface: it multiplies with complex matrices +module ContainerElementBase_mod + implicit none + + private + public :: ContainerElementBase + + ! Base for defining the interface + type, abstract :: ContainerElementBase + contains + procedure(rmultinterface), deferred :: rmult + procedure(lmultinterface), deferred :: lmult + procedure(rmultinvinterface), deferred :: rmultinv + procedure(lmultinvinterface), deferred :: lmultinv + procedure(adjointactioninterface), deferred :: adjointaction + procedure(dump), deferred :: dump + procedure(dealloc), deferred :: dealloc + end type ContainerElementBase + + abstract interface + + + !-------------------------------------------------------------------- + !> @brief + !> multiplies this with arg from the right. + ! + !> @param[in] this + !-------------------------------------------------------------------- + subroutine rmultinterface(this, arg) + import ContainerElementBase + class(ContainerElementBase), intent(in) :: this + Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: arg + end subroutine + + !-------------------------------------------------------------------- + !> @brief + !> multiplies this with arg from the left. + ! + !> @param[in] this + !-------------------------------------------------------------------- + subroutine lmultinterface(this, arg) + import ContainerElementBase + class(ContainerElementBase), intent(in) :: this + Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: arg + end subroutine + + !-------------------------------------------------------------------- + !> @brief + !> multiplies this^-1 with arg from the right. + ! + !> @param[in] this + !-------------------------------------------------------------------- + subroutine rmultinvinterface(this, arg) + import ContainerElementBase + class(ContainerElementBase), intent(in) :: this + Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: arg + end subroutine + + !-------------------------------------------------------------------- + !> @brief + !> multiplies this^-1 with arg from the left. + ! + !> @param[in] this + !-------------------------------------------------------------------- + subroutine lmultinvinterface(this, arg) + import ContainerElementBase + class(ContainerElementBase), intent(in) :: this + Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: arg + end subroutine + + !-------------------------------------------------------------------- + !> @brief + !> This dumps the content to the screen. + ! + !> @param[in] this + !-------------------------------------------------------------------- + subroutine dump(this) + import ContainerElementBase + class(ContainerElementBase), intent(in) :: this + end subroutine + + !-------------------------------------------------------------------- + !> @brief + !> Free the used memory + ! + !> @param[in] this + !-------------------------------------------------------------------- + subroutine dealloc(this) + import ContainerElementBase + class(ContainerElementBase), intent(inout) :: this + end subroutine + + !-------------------------------------------------------------------- + !> @brief + !> Perform the similarity transform e^{-T/2} arg e^{T/2} + ! + !> @param[in] this + !> @param[inout] the matrix that we intend to transform. + !-------------------------------------------------------------------- + subroutine adjointactioninterface(this, arg) + import ContainerElementBase + class(ContainerElementBase), intent(in) :: this + Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: arg + end subroutine + end interface +end module ContainerElementBase_mod diff --git a/Prog/DynamicMatrixArray_mod.F90 b/Prog/DynamicMatrixArray_mod.F90 new file mode 100644 index 0000000000000000000000000000000000000000..577f20b894e48744fdfcda6e5b6d9f57fa291948 --- /dev/null +++ b/Prog/DynamicMatrixArray_mod.F90 @@ -0,0 +1,179 @@ +! Copyright (C) 2020 The ALF project +! +! The ALF project is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! The ALF project is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with ALF. If not, see http://www.gnu.org/licenses/. +! +! Under Section 7 of GPL version 3 we require you to fulfill the following additional terms: +! +! - It is our hope that this program makes a contribution to the scientific community. Being +! part of that community we feel that it is reasonable to require you to give an attribution +! back to the original authors if you have benefitted from this program. +! Guidelines for a proper citation can be found on the project's homepage +! https://alf.physik.uni-wuerzburg.de . +! +! - We require the preservation of the above copyright notice and this license in all original files. +! +! - We prohibit the misrepresentation of the origin of the original source files. To obtain +! the original source files please visit the homepage https://alf.physik.uni-wuerzburg.de . +! +! - If you make substantial changes to the program we require you to either consider contributing +! to the ALF project or to mark your material in a reasonable way as different from the original version. + +! helpful Fortran docs: +! https://materials.prace-ri.eu/400/1/advFortranIntro.pdf +! http://www.chem.helsinki.fi/~manninen/fortran2014/7_Object_oriented_features.pdf + +module DynamicMatrixArray_mod + Use ContainerElementBase_mod + implicit none + + private + public :: DynamicMatrixArray + + type :: OpTBasePtrWrapper + class(ContainerElementBase), pointer :: dat => null() + end type + + type :: DynamicMatrixArray + integer :: avamem ! amount of available space + integer :: tail ! last valid Fortran index + Type(OpTbasePtrWrapper), allocatable, dimension(:) :: data ! actual effective array of the pointers + contains + procedure :: init => DynamicMatrixArray_init + procedure :: dealloc => DynamicMatrixArray_dealloc + procedure :: pushback => DynamicMatrixArray_pushback + procedure :: at => DynamicMatrixArray_at + procedure :: back => DynamicMatrixArray_back + procedure :: length => DynamicMatrixArray_length + ! FIXME: do we need insert? + end type DynamicMatrixArray + +contains + +!-------------------------------------------------------------------- +!> @author +!> The ALF Project contributors +! +!> @brief +!> set up initial state of the vector. +!> +!> @param[inout] this the vector +!-------------------------------------------------------------------- +subroutine DynamicMatrixArray_init(this) + class(DynamicMatrixArray) :: this + type(OpTbasePtrWrapper) :: temp + this%tail = 0 !when the vector has no content this is invalid memory + this%avamem = 4096/(STORAGE_SIZE(temp)/8) ! allocate a page of memory ! Note STORAGE_SIZE: F2008, SIZEOF: GCC Extension + allocate(this%data(this%avamem)) +end subroutine DynamicMatrixArray_init + +!-------------------------------------------------------------------- +!> @author +!> The ALF Project contributors +! +!> @brief +!> Deallocates internal storage. Pointees are not deleted! +!> +!> @param[inout] this the vector +!-------------------------------------------------------------------- +subroutine DynamicMatrixArray_dealloc(this) + class(DynamicMatrixArray) :: this + deallocate(this%data) +end subroutine + +!-------------------------------------------------------------------- +!> @author +!> The ALF Project contributors +! +!> @brief +!> Attach a pointer to the object given by itm at the end of the vector. +!> If out of space the vector grows. +!> +!> @param[inout] this the vector +!> @param[in] itm the object that we like to store a pointer to. +! +!-------------------------------------------------------------------- +subroutine DynamicMatrixArray_pushback(this, itm) + class(DynamicMatrixArray) :: this + class(ContainerElementBase), intent(in), target :: itm !Type(...) has to match exactly, class(...) allows for polymorphism + type(OpTbasePtrWrapper), allocatable, dimension(:) :: temp + integer :: i + + if (this%tail == this%avamem) then + ! reallocate the memory + write (*,*) "not enough space -> growing." + call MOVE_ALLOC(this%data, temp) + allocate(this%data(2*this%avamem)) + do i = 1, this%avamem + this%data(i) = temp(i) + enddo + deallocate(temp) + this%avamem = 2*this%avamem + endif + this%tail = this%tail + 1 + this%data(this%tail)%dat => itm ! let the pointer point to the object +end subroutine + +!-------------------------------------------------------------------- +!> @author +!> The ALF Project contributors +! +!> @brief +!> return a pointer to the object stored at position i. +!> +!> @param[in] i the index. +!> @param[out] the content stored at the position i. +! +!-------------------------------------------------------------------- +function DynamicMatrixArray_at(this, pos) result(itm) + class(DynamicMatrixArray), intent(in) :: this + integer, intent(in) :: pos + class(ContainerElementBase), pointer :: itm + itm => this%data(pos)%dat +end function + +!-------------------------------------------------------------------- +!> @author +!> The ALF Project contributors +! +!> @brief +!> returns the pointer to the last element +!> +!> @param[inout] this the vector. +!> @return itm the element at the end of the vector. +! +!-------------------------------------------------------------------- +function DynamicMatrixArray_back(this) result(itm) + class(DynamicMatrixArray), intent(in) :: this + class(ContainerElementBase), pointer :: itm + itm => this%data(this%tail)%dat +end function + +!-------------------------------------------------------------------- +!> @author +!> The ALF Project contributors +! +!> @brief +!> Inquire the length of the vector. +!> +!> @param[inout] this the vector +!> @return the current length of the vector +! +!-------------------------------------------------------------------- +function DynamicMatrixArray_length(this) result(l) + class(DynamicMatrixArray) :: this + integer :: l + l = this%tail +end function + +end module DynamicMatrixArray_mod diff --git a/Prog/Hop_mod.F90 b/Prog/Hop_mod.F90 index 4ab1b75e18338055549c404139051f73fa0ba982..3da8911bc66dba5a7b83a52df43b223c3569ad6f 100644 --- a/Prog/Hop_mod.F90 +++ b/Prog/Hop_mod.F90 @@ -48,16 +48,55 @@ Use Hamiltonian Use Random_wrap + Use DynamicMatrixArray_mod + Use ContainerElementBase_mod + Use OpTTypes_mod use iso_fortran_env, only: output_unit, error_unit ! Private variables - Complex (Kind=Kind(0.d0)), allocatable, private :: Exp_T(:,:,:,:), Exp_T_M1(:,:,:,:) - Complex (Kind=Kind(0.d0)), allocatable, private :: Exp_T_1D2(:,:,:,:), Exp_T_M1_1D2(:,:,:,:) - Complex (Kind=Kind(0.d0)), allocatable, private :: U_HLP(:,:), U_HLP1(:,:), V_HLP(:,:), V_HLP1(:,:) - Integer, private, save :: Ncheck, Ndim_hop + Type(DynamicMatrixArray), private, allocatable :: ExpOpT_vec(:) ! for now we have for simplicity for each flavour a vector + Integer, private, save :: Ncheck Real (Kind=Kind(0.d0)), private, save :: Zero Contains + +!-------------------------------------------------------------------- +!> @author +!> ALF-project +! +!> @brief +!> This function serves as a central entry point to collect the +!> processing that occurs in mapping an OpT input matrix to the internal +!> matrix-like data structure. +! +!> @param ExpOpT_vec[inout] a DynamicMatrixArray structure to which we append new elements. +!> @param op[in] an Operator that describes an OpT hopping matrix. +! +!-------------------------------------------------------------------- + subroutine OpT_postprocess(ExpOpT_vec, op) + use Operator_mod + implicit none + + Type(DynamicMatrixArray), intent(inout) :: ExpOpT_vec + Type(Operator), intent(in) :: op + + Class(CmplxExpOpT), pointer :: cmplxexp => null() + Class(RealExpOpT), pointer :: realexp => null() + + if (Op_is_real(op)) then + ! branch for real operators + allocate(realexp) ! Yep, this is a manifest memory leak. Using the ptr we can allocate onto the same variable + call realexp%init(op) + call ExpOpT_vec%pushback(realexp) + else + ! branch for complex operators + allocate(cmplxexp) + call cmplxexp%init(op) + call ExpOpT_vec%pushback(cmplxexp) + endif + end subroutine + + !-------------------------------------------------------------------- !> @author !> ALF-project @@ -68,62 +107,26 @@ ! !-------------------------------------------------------------------- subroutine Hop_mod_init - Implicit none - Integer :: nc, nf, i,j - Complex (Kind=Kind(0.d0)) :: g + Integer :: nc, nf Ncheck = size(Op_T,1) If ( size(Op_T,2) /= N_FL ) then Write(error_unit,*) 'Hop_mod_init: Error in the number of flavors.' error stop 1 Endif - Ndim_hop = Op_T(1,1)%N - !Write(6,*) 'In Hop_mod: ', Ndim, Ndim_hop, Ncheck - Do nc = 1, Ncheck - do nf = 1,N_FL - if ( Ndim_hop /= Op_T(nc,nf)%N ) Then - Write(error_unit,*) 'Hop_mod_init: Different size of Hoppings not implemented ' - error stop 1 - endif - enddo - enddo - - Allocate ( Exp_T (Ndim_hop,Ndim_hop,Ncheck,N_FL) ) - Allocate ( Exp_T_M1 (Ndim_hop,Ndim_hop,Ncheck,N_FL) ) - Allocate ( Exp_T_1D2 (Ndim_hop,Ndim_hop,Ncheck,N_FL) ) - Allocate ( Exp_T_M1_1D2(Ndim_hop,Ndim_hop,Ncheck,N_FL) ) - Allocate ( V_Hlp(Ndim_hop,Ndim) ) - Allocate ( V_Hlp1(Ndim_hop,Ndim) ) - Allocate ( U_Hlp (Ndim, Ndim_hop) ) - Allocate ( U_Hlp1(Ndim, Ndim_hop) ) + allocate(ExpOpT_vec(N_FL)) - Exp_T = cmplx(0.d0, 0.d0, kind(0.D0)) - Exp_T_M1 = cmplx(0.d0, 0.d0, kind(0.D0)) do nf = 1,N_FL + call ExpOpT_vec(nf)%init() do nc = 1,Ncheck - g = Op_T(nc,nf)%g - Call Op_exp(g,Op_T(nc,nf),Exp_T(:,:,nc,nf)) - g = -Op_T(nc,nf)%g - Call Op_exp(g,Op_T(nc,nf),Exp_T_M1(:,:,nc,nf)) - g = Op_T(nc,nf)%g/2.d0 - Call Op_exp(g,Op_T(nc,nf),Exp_T_1D2(:,:,nc,nf)) - g = -Op_T(nc,nf)%g/2.d0 - Call Op_exp(g,Op_T(nc,nf),Exp_T_M1_1D2(:,:,nc,nf)) - ! symmetrize the upper part of Exp_T and Exp_T_M1 - DO i = 1, Ndim_hop - DO j = i, Ndim_hop - Exp_T(i, j, nc, nf) = (Exp_T(i, j, nc, nf) + Conjg(Exp_T(j, i, nc, nf)))/2.D0 - Exp_T_M1(i, j, nc, nf) = (Exp_T_M1(i, j, nc, nf) + Conjg(Exp_T_M1(j, i, nc, nf)))/2.D0 - ENDDO - ENDDO + call OpT_postprocess(ExpOpT_vec(nf), Op_T(nc, nf)) enddo enddo Zero = 1.E-12 - end subroutine Hop_mod_init !-------------------------------------------------------------------- @@ -138,43 +141,15 @@ Integer, intent(IN) :: nf !Local - Integer :: nc, N1, N2 - - N1=size(In,1) - N2=size(In,2) + Integer :: nc + class(ContainerElementBase), pointer :: dummy do nc = Ncheck,1,-1 - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - call ZSLHEMM('L','U',Ndim_hop,N1,N2,Exp_T(:,:,nc,nf),Op_T(nc,nf)%P,In) - Endif + dummy => ExpOpT_vec(nf)%at(nc) + call dummy%lmult(In) Enddo end Subroutine Hop_mod_mmthr -!-------------------------------------------------------------------- - Subroutine Hop_mod_mmthr_1D2(In,nf) - - - ! InOut: In = e^{ -dtau T /2 }.IN - Implicit none - - Complex (Kind=Kind(0.d0)), intent(INOUT) :: IN(:,:) - Integer, intent(IN) :: nf - - !Local - Integer :: nc, N1, N2 - - N1=size(In,1) - N2=size(In,2) - - do nc = Ncheck,1,-1 - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - call ZSLHEMM('L','U',Ndim_hop,N1,N2,Exp_T_1D2(:,:,nc,nf),Op_T(nc,nf)%P,In) - Endif - Enddo - end Subroutine Hop_mod_mmthr_1D2 - -!-------------------------------------------------------------------- - Subroutine Hop_mod_mmthr_m1(In,nf) @@ -185,15 +160,12 @@ Integer :: nf !Local - Integer :: nc , N1, N2 - - N1=size(In,1) - N2=size(In,2) + Integer :: nc + class(ContainerElementBase), pointer :: dummy do nc = 1,Ncheck - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - call ZSLHEMM('L','U',Ndim_hop,N1,N2,Exp_T_m1(:,:,nc,nf),Op_T(nc,nf)%P,In) - Endif + dummy => ExpOpT_vec(nf)%at(nc) + call dummy%lmultinv(In) Enddo end Subroutine Hop_mod_mmthr_m1 @@ -210,15 +182,12 @@ Integer :: nf !Local - Integer :: nc, N1, N2 - - N1=size(In,1) - N2=size(In,2) + Integer :: nc + class(ContainerElementBase), pointer :: dummy do nc = 1, Ncheck - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - call ZSLHEMM('R','U',Ndim_hop,N1,N2,Exp_T(:,:,nc,nf),Op_T(nc,nf)%P,In) - Endif + dummy => ExpOpT_vec(nf)%at(nc) + call dummy%rmult(In) Enddo end Subroutine Hop_mod_mmthl @@ -235,15 +204,12 @@ Integer :: nf !Local - Integer :: nc, N1, N2 - - N1=size(In,1) - N2=size(In,2) + Integer :: nc + class(ContainerElementBase), pointer :: dummy do nc = 1, Ncheck - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - call ZSLHEMM('L','U',Ndim_hop,N1,N2,Exp_T(:,:,nc,nf),Op_T(nc,nf)%P,In) - Endif + dummy => ExpOpT_vec(nf)%at(nc) + call dummy%lmult(In) Enddo end Subroutine Hop_mod_mmthlc @@ -260,46 +226,16 @@ Integer :: nf !Local - Integer :: nc, N1, N2 - - N1=size(In,1) - N2=size(In,2) + Integer :: nc + class(ContainerElementBase), pointer :: dummy do nc = Ncheck,1,-1 - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - call ZSLHEMM('R','U',Ndim_hop,N1,N2,Exp_T_m1(:,:,nc,nf),Op_T(nc,nf)%P,In) - Endif + dummy => ExpOpT_vec(nf)%at(nc) + call dummy%rmultinv(In) Enddo end Subroutine Hop_mod_mmthl_m1 - -!-------------------------------------------------------------------- - - Subroutine Hop_mod_mmthl_m1_1D2(In, nf) - - - ! InOut: In = IN * e^{ dtau T/2 } - Implicit none - - Complex (Kind=Kind(0.d0)), intent(INOUT) :: IN(:,:) - Integer :: nf - - !Local - Integer :: nc, N1, N2 - - N1=size(In,1) - N2=size(In,2) - - do nc = Ncheck,1,-1 - If ( dble( Op_T(nc,nf)%g*conjg(Op_T(nc,nf)%g) ) > Zero ) then - call ZSLHEMM('R','U',Ndim_hop,N1,N2,Exp_T_m1_1D2(:,:,nc,nf),Op_T(nc,nf)%P,In) - Endif - Enddo - - end Subroutine Hop_mod_mmthl_m1_1D2 - - !!$Subroutine Hop_mod_test !!$ !!$Implicit none @@ -333,12 +269,15 @@ COMPLEX (Kind=Kind(0.d0)), Dimension(:,:,:), Intent(Out):: Out COMPLEX (Kind=Kind(0.d0)), Dimension(:,:,:), Intent(IN):: In - Integer :: nf + Integer :: nf, nc + class(ContainerElementBase), pointer :: dummy Out = In Do nf = 1, size(In,3) - Call Hop_mod_mmthr_1D2 (Out(:,:,nf), nf ) - Call Hop_mod_mmthl_m1_1D2(Out(:,:,nf), nf ) + do nc = Ncheck,1,-1 + dummy => ExpOpT_vec(nf)%at(nc) + call dummy%adjointaction(Out(:, :, nf)) + enddo enddo End Subroutine Hop_mod_Symm diff --git a/Prog/Makefile b/Prog/Makefile index dd1a2b00302aca4923d1802860907de7e2afb33d..1a5e4d41946a06727cc5026329fed60614868ccc 100644 --- a/Prog/Makefile +++ b/Prog/Makefile @@ -1,6 +1,6 @@ .PHONY : Compile tidy clean Z2_Matter Kondo Hubbard tV Hubbard_Plain_Vanilla LRC -OBJS= Hamiltonians/LRC_mod.o Set_random.o control_mod.o Fields_mod.o Operator_mod.o WaveFunction_mod.o observables_mod.o \ +OBJS= Hamiltonians/LRC_mod.o Set_random.o control_mod.o Fields_mod.o Operator_mod.o WaveFunction_mod.o observables_mod.o ContainerElementBase_mod.o DynamicMatrixArray_mod.o OpTTypes_mod.o\ Predefined_Int_mod.o Predefined_Obs_mod.o Predefined_Latt_mod.o Predefined_Hop_mod.o Predefined_Trial_mod.o \$(HAMILTONIAN) QDRP_decompose_mod.o udv_state_mod.o Hop_mod.o UDV_WRAP_mod.o \ Global_mod.o Langevin_HMC_mod.o Wrapgr_mod.o tau_m_mod.o tau_p_mod.o main.o wrapul.o cgr1.o wrapur.o cgr2_2.o upgrade.o @@ -8,7 +8,7 @@ MODS= control.mod fields_mod.mod global_mod.mod hop_mod.mod lrc_mod.mod observab operator_mod.mod predefined_int.mod predefined_lattices.mod predefined_obs.mod \ predefined_hoppings.mod predefined_trial.mod qdrp_mod.mod tau_m_mod.mod tau_p_mod.mod \ udv_state_mod.mod udv_wrap_mod.mod wavefunction_mod.mod wrapgr_mod.mod hamiltonian.mod \ - Langevin_HMC_mod.mod + containerelementbase_mod.mod opttypes_mod.mod dynamicmatrixarray_mod.mod OHAM= Hamiltonians/Hamiltonian_Z2_Matter_mod.o \ Hamiltonians/Hamiltonian_Kondo_mod.o \ Hamiltonians/Hamiltonian_Hubbard_mod.o \ diff --git a/Prog/OpTTypes_mod.F90 b/Prog/OpTTypes_mod.F90 new file mode 100644 index 0000000000000000000000000000000000000000..0d6610349bfc40519a71be140f0abfdd48cc4691 --- /dev/null +++ b/Prog/OpTTypes_mod.F90 @@ -0,0 +1,331 @@ +! Copyright (C) 2020 The ALF project +! +! The ALF project is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! The ALF project is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with ALF. If not, see http://www.gnu.org/licenses/. +! +! Under Section 7 of GPL version 3 we require you to fulfill the following additional terms: +! +! - It is our hope that this program makes a contribution to the scientific community. Being +! part of that community we feel that it is reasonable to require you to give an attribution +! back to the original authors if you have benefitted from this program. +! Guidelines for a proper citation can be found on the project's homepage +! https://alf.physik.uni-wuerzburg.de . +! +! - We require the preservation of the above copyright notice and this license in all original files. +! +! - We prohibit the misrepresentation of the origin of the original source files. To obtain +! the original source files please visit the homepage https://alf.physik.uni-wuerzburg.de . +! +! - If you make substantial changes to the program we require you to either consider contributing +! to the ALF project or to mark your material in a reasonable way as different from the original version. + + +module OpTTypes_mod + use ContainerElementBase_mod + use Operator_mod + implicit none + + private + public :: RealExpOpT, CmplxExpOpT + + !-------------------------------------------------------------------- + !> @author + !> ALF-project + !> @brief + !> Encapsulates operations for real exponentiated OpTs. + !> + !-------------------------------------------------------------------- + type, extends(ContainerElementBase) :: RealExpOpT + Real(kind=kind(0.d0)), allocatable, dimension(:,:) :: mat, invmat, mat_1D2, invmat_1D2 !>We store the matrix in the class + Real(kind=kind(0.d0)) :: g, Zero + integer, pointer :: P(:) + Integer :: m, n, Ndim_hop + + contains + procedure :: init => RealExpOpT_init ! initialize and allocate matrices + procedure :: dealloc => RealExpOpT_dealloc ! dealloc matrices + procedure :: rmult => RealExpOpT_rmult ! right multiplication with Op_T + procedure :: lmult => RealExpOpT_lmult + procedure :: rmultinv => RealExpOpT_rmultinv ! right multiplication with Op_T inverse + procedure :: lmultinv => RealExpOpT_lmultinv + procedure :: adjointaction => RealExpOpT_adjointaction + procedure :: dump => RealExpOpT_dump ! dump matrices for debugging to screen + end type RealExpOpT + + !-------------------------------------------------------------------- + !> @author + !> ALF-project + !> @brief + !> Encapsulates operations for complex exponentiated OpTs. + !> + !-------------------------------------------------------------------- + type, extends(ContainerElementBase) :: CmplxExpOpT + Complex(kind=kind(0.d0)), allocatable, dimension(:,:) :: mat, invmat, mat_1D2, invmat_1D2 !>We store the matrix inclass + Complex(kind=kind(0.d0)) :: g + Real(kind=kind(0.d0)) :: Zero + integer, pointer :: P(:) + Integer :: m, n, Ndim_hop + contains + procedure :: init => CmplxExpOpT_init ! initialize and allocate matrices + procedure :: dealloc => CmplxExpOpT_dealloc ! dealloc matrices + procedure :: rmult => CmplxExpOpT_rmult ! right multiplication with Op_T + procedure :: lmult => CmplxExpOpT_lmult + procedure :: rmultinv => CmplxExpOpT_rmultinv ! right multiplication with Op_T inverse + procedure :: lmultinv => CmplxExpOpT_lmultinv + procedure :: adjointaction => CmplxExpOpT_adjointaction + procedure :: dump => CmplxExpOpT_dump ! dump matrices for debugging to screen + end type CmplxExpOpT + +contains + subroutine RealExpOpT_init(this, Op_T) + class(RealExpOpT) :: this + Type(Operator), intent(in) :: Op_T + Complex(kind=kind(0.d0)), allocatable, dimension(:,:) :: cmat, cinvmat + Complex(kind=kind(0.d0)) :: cg + Integer :: i, j + + this%Zero = 1.E-12 + this%Ndim_hop = Op_T%N + allocate(this%mat(this%Ndim_hop, this%Ndim_hop), this%invmat(this%Ndim_hop, this%Ndim_hop)) + allocate(this%mat_1D2(this%Ndim_hop, this%Ndim_hop), this%invmat_1D2(this%Ndim_hop, this%Ndim_hop)) + allocate(cmat(this%Ndim_hop, this%Ndim_hop), cinvmat(this%Ndim_hop, this%Ndim_hop)) + + cg = -Op_T%g + Call Op_exp(cg, Op_T, cinvmat) + cg = Op_T%g + Call Op_exp(cg, Op_T, cmat) + ! copy over the data to the real storage + this%mat = DBLE(cmat) + this%invmat = DBLE(cinvmat) + + cg = -Op_T%g/2.0 + Call Op_exp(cg, Op_T, cinvmat) + cg = Op_T%g/2.0 + Call Op_exp(cg, Op_T, cmat) + ! copy over the data to the real storage + this%mat_1D2 = DBLE(cmat) + this%invmat_1D2 = DBLE(cinvmat) + + DO i = 1, this%Ndim_hop + DO j = i, this%Ndim_hop + this%mat(i, j) = (this%mat(i, j) + this%mat(j, i))/2.D0 + this%invmat(i, j) = (this%invmat(i, j) + this%invmat(j, i))/2.D0 + + this%mat_1D2(i, j) = (this%mat_1D2(i, j) + this%mat_1D2(j, i))/2.D0 + this%invmat_1D2(i, j) = (this%invmat_1D2(i, j) + this%invmat_1D2(j, i))/2.D0 + ENDDO + ENDDO + this%P => Op_T%P + this%g = DBLE(Op_T%g) + deallocate(cmat, cinvmat) + end subroutine + + subroutine RealExpOpT_adjointaction(this, arg) + class(RealExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Integer :: n1, n2 + + n1 = size(arg,1) + n2 = size(arg,2) + If ( this%g*this%g > this%Zero ) then + call ZDSLSYMM('L', 'U', this%Ndim_hop, n1, n2, this%mat_1D2, this%P, arg) + call ZDSLSYMM('R', 'U', this%Ndim_hop, n1, n2, this%invmat_1D2, this%P, arg) + Endif + end subroutine + + subroutine RealExpOpT_rmult(this, arg) + class(RealExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Integer :: n1, n2 + + n1 = size(arg,1) + n2 = size(arg,2) + If ( this%g*this%g > this%Zero ) then + call ZDSLSYMM('R', 'U', this%Ndim_hop, n1, n2, this%mat, this%P, arg) + Endif + end subroutine + + subroutine RealExpOpT_rmultinv(this, arg) + class(RealExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Integer :: n1, n2 + + n1 = size(arg,1) + n2 = size(arg,2) + If ( this%g*this%g > this%Zero ) then + call ZDSLSYMM('R', 'U', this%Ndim_hop, n1, n2, this%invmat, this%P, arg) + Endif + end subroutine + + subroutine RealExpOpT_lmult(this, arg) + class(RealExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + integer :: n1, n2 + + ! taken from mmthr + n1 = size(arg,1) + n2 = size(arg,2) + If ( this%g*this%g > this%Zero ) then + call ZDSLSYMM('L', 'U', this%Ndim_hop, n1, n2, this%mat, this%P, arg) + Endif + end subroutine + + subroutine RealExpOpT_lmultinv(this, arg) + class(RealExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + integer :: n1, n2 + + n1 = size(arg,1) + n2 = size(arg,2) + If ( this%g*this%g > this%Zero ) then + call ZDSLSYMM('L', 'U', this%Ndim_hop, n1, n2, this%invmat, this%P, arg) + Endif + end subroutine + + subroutine CmplxExpOpT_init(this, Op_T) + class(CmplxExpOpT) :: this + Type(Operator), intent(in) :: Op_T + Integer :: i, j + + this%Zero = 1.E-12 + this%Ndim_hop = Op_T%N + + allocate(this%mat(this%Ndim_hop, this%Ndim_hop), this%invmat(this%Ndim_hop, this%Ndim_hop)) + allocate(this%mat_1D2(this%Ndim_hop, this%Ndim_hop), this%invmat_1D2(this%Ndim_hop, this%Ndim_hop)) + + this%g = -Op_T%g/2.0 + Call Op_exp(this%g, Op_T, this%invmat_1D2) + this%g = Op_T%g/2.0 + Call Op_exp(this%g, Op_T, this%mat_1D2) + + this%g = -Op_T%g + Call Op_exp(this%g, Op_T, this%invmat) + this%g = Op_T%g + Call Op_exp(this%g, Op_T, this%mat) + + DO i = 1, this%Ndim_hop + DO j = i, this%Ndim_hop + this%mat(i, j) = (this%mat(i, j) + Conjg(this%mat(j, i)))/2.D0 + this%invmat(i, j) = (this%invmat(i, j) + Conjg(this%invmat(j, i)))/2.D0 + + this%mat_1D2(i, j) = (this%mat_1D2(i, j) + Conjg(this%mat_1D2(j, i)))/2.D0 + this%invmat_1D2(i, j) = (this%invmat_1D2(i, j) + Conjg(this%invmat_1D2(j, i)))/2.D0 + ENDDO + ENDDO + this%P => Op_T%P + + end subroutine + + subroutine CmplxExpOpT_adjointaction(this, arg) + class(CmplxExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Integer :: n1, n2 + + n1 = size(arg,1) + n2 = size(arg,2) + If ( dble(this%g*conjg(this%g)) > this%Zero ) then + call ZSLHEMM('L', 'U', this%Ndim_hop, n1, n2, this%mat_1D2, this%P, arg) + call ZSLHEMM('R', 'U', this%Ndim_hop, n1, n2, this%invmat_1D2, this%P, arg) + Endif + + end subroutine + + subroutine CmplxExpOpT_rmult(this, arg) + class(CmplxExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Integer :: n1, n2 + + ! taken from mmthl + n1 = size(arg,1) + n2 = size(arg,2) + If ( dble(this%g*conjg(this%g)) > this%Zero ) then + call ZSLHEMM('R', 'U', this%Ndim_hop, n1, n2, this%mat, this%P, arg) + Endif + end subroutine + + subroutine CmplxExpOpT_rmultinv(this, arg) + class(CmplxExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Integer :: n1, n2 + + ! taken from mmthl_m1 + n1 = size(arg,1) + n2 = size(arg,2) + If ( dble(this%g*conjg(this%g)) > this%Zero ) then + call ZSLHEMM('R', 'U', this%Ndim_hop, n1, n2, this%invmat, this%P, arg) + Endif + end subroutine + + subroutine CmplxExpOpT_lmult(this, arg) + class(CmplxExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + integer :: n1, n2 + + ! taken from mmthr + n1 = size(arg,1) + n2 = size(arg,2) + If ( dble(this%g*conjg(this%g)) > this%Zero ) then + call ZSLHEMM('L', 'U', this%Ndim_hop, n1, n2, this%mat, this%P, arg) + Endif + end subroutine + + subroutine CmplxExpOpT_lmultinv(this, arg) + class(CmplxExpOpT), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + integer :: n1, n2 + n1 = size(arg,1) + n2 = size(arg,2) + If ( dble(this%g*conjg(this%g)) > this%Zero ) then + call ZSLHEMM('L', 'U', this%Ndim_hop, n1, n2, this%invmat, this%P, arg) + Endif + end subroutine + + subroutine CmplxExpOpT_dump(this) + class(CmplxExpOpT), intent(in) :: this + integer :: i,j + + do i = 1, size(this%mat, 1) + write (*,*) (dble(this%mat(i,j)), j = 1,size(this%mat,2) ) + enddo + write (*,*) "---------------" + do i = 1, size(this%mat, 1) + write (*,*) (dble(this%invmat(i,j)), j = 1,size(this%mat,2) ) + enddo + end subroutine + + subroutine RealExpOpT_dump(this) + class(RealExpOpT), intent(in) :: this + integer :: i,j + + do i = 1, size(this%mat, 1) + write (*,*) (dble(this%mat(i,j)), j = 1,size(this%mat,2) ) + enddo + write (*,*) "---------------" + do i = 1, size(this%mat, 1) + write (*,*) (dble(this%invmat(i,j)), j = 1,size(this%mat,2) ) + enddo + end subroutine + + subroutine CmplxExpOpT_dealloc(this) + class(CmplxExpOpT), intent(inout) :: this + + deallocate(this%mat, this%invmat, this%mat_1D2, this%invmat_1D2) + end subroutine + + subroutine RealExpOpT_dealloc(this) + class(RealExpOpT), intent(inout) :: this + + deallocate(this%mat, this%invmat, this%mat_1D2, this%invmat_1D2) + end subroutine + +end module OpTTypes_mod diff --git a/Prog/Operator_mod.F90 b/Prog/Operator_mod.F90 index 0f940cb7b27b42684736139bb5067456bddbb93e..c8b86a96d2d74d74bca56e25e88082e0f9a39d44 100644 --- a/Prog/Operator_mod.F90 +++ b/Prog/Operator_mod.F90 @@ -420,7 +420,6 @@ Contains !-------------------------------------------------------------------- subroutine Op_mmultL(Mat,Op,spin,cop) Implicit none - Integer :: Ndim Type (Operator) , INTENT(IN) :: Op Complex (Kind=Kind(0.d0)), INTENT(INOUT) :: Mat (:,:) Real (Kind=Kind(0.d0)), INTENT(IN) :: spin @@ -654,7 +653,7 @@ Contains Implicit none Integer :: Ndim - Type (Operator) , INTENT(IN ) :: Op + Type (Operator) , INTENT(IN) :: Op Complex (Kind = Kind(0.D0)), INTENT(INOUT) :: Mat (Ndim,Ndim) Real (Kind=Kind(0.d0)), INTENT(IN ) :: spin Integer, INTENT(IN) :: N_Type @@ -712,5 +711,23 @@ Contains endif endif end Subroutine Op_Wrapdo - + + function Op_is_real(Op) result(retval) + Implicit None + + Type (Operator) , INTENT(IN) :: Op + Logical :: retval + Real (Kind=Kind(0.d0)) :: myzero + integer :: i,j + + retval = (Abs(aimag(Op%g)) < Abs(Op%g)*epsilon(1.D0)) + ! calculate a matrix scale + myzero = maxval(abs(Op%E))*epsilon(Op%E) + + do i = 1, Op%N + do j = 1, Op%N + retval = retval .and. (Abs(aimag(Op%O(i,j))) < myzero) + enddo + enddo + end function Op_is_real end Module Operator_mod diff --git a/Prog/matTypes_mod.F90 b/Prog/matTypes_mod.F90 new file mode 100644 index 0000000000000000000000000000000000000000..3c1ae70e3a74d1114fc973bb4e3406a9ab96e319 --- /dev/null +++ b/Prog/matTypes_mod.F90 @@ -0,0 +1,183 @@ +! Copyright (C) 2020 The ALF project +! +! The ALF project is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! The ALF project is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with ALF. If not, see http://www.gnu.org/licenses/. +! +! Under Section 7 of GPL version 3 we require you to fulfill the following additional terms: +! +! - It is our hope that this program makes a contribution to the scientific community. Being +! part of that community we feel that it is reasonable to require you to give an attribution +! back to the original authors if you have benefitted from this program. +! Guidelines for a proper citation can be found on the project's homepage +! https://alf.physik.uni-wuerzburg.de . +! +! - We require the preservation of the above copyright notice and this license in all original files. +! +! - We prohibit the misrepresentation of the origin of the original source files. To obtain +! the original source files please visit the homepage https://alf.physik.uni-wuerzburg.de . +! +! - If you make substantial changes to the program we require you to either consider contributing +! to the ALF project or to mark your material in a reasonable way as different from the original version. + + +module matTypes_mod + use ContainerElementBase_mod + implicit none + + private + public :: RealMat, CmplxMat + + type, extends(ContainerElementBase) :: RealMat + Real(kind=kind(0.d0)), allocatable, dimension(:,:) :: mat + Integer :: m, n + contains + procedure :: init => RealMat_init + procedure :: simt => RealMat_simt + procedure :: rmult => RealMat_rmult + procedure :: lmult => RealMat_lmult + procedure :: rmultinv => RealMat_rmultinv + procedure :: lmultinv => RealMat_lmultinv + end type RealMat + + type, extends(ContainerElementBase) :: CmplxMat + Complex(kind=kind(0.d0)),allocatable, dimension(:,:):: mat + Integer :: m, n + contains + procedure :: init => CmplxMat_init + procedure :: simt => CmplxMat_simt + procedure :: rmult => CmplxMat_rmult + procedure :: lmult => CmplxMat_lmult + procedure :: rmultinv => CmplxMat_rmultinv + procedure :: lmultinv => CmplxMat_lmultinv + end type CmplxMat + +contains + subroutine RealMat_init(this, arg) + class(RealMat) :: this + Real(kind=kind(0.D0)), intent(inout), allocatable, dimension(:,:) :: arg + Integer :: i,j + this%mat = arg !copy argument to local storage + this%m = size(arg, 1) + this%n = size(arg, 2) + +! do i = 1, size(arg,1) +! write (*,*) (this%mat(i,j), j = 1, size(arg,2) ) +! enddo + end subroutine + + subroutine RealMat_simt(this, arg) + class(RealMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), allocatable, dimension(:,:) :: arg + Complex(kind=kind(0.D0)), allocatable, dimension(:,:) :: temp + + end subroutine + + subroutine RealMat_rmult(this, arg) + class(RealMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Complex(kind=kind(0.D0)), allocatable, dimension(:,:) :: out + Real(kind=kind(0.D0)), allocatable, dimension(:) :: rwork + Integer :: i, j, sz1, sz2 + + sz1 = size(arg, 1) + sz2 = size(arg, 2) + allocate(out(sz1, sz2), rwork(2*sz1*sz2)) + call zlacrm(sz1, sz2, arg, sz1, this%mat, this%m, out, this%m, rwork) ! zlarcm assumes mat to be square + arg = out + deallocate(out, rwork) + end subroutine + + subroutine RealMat_rmultinv(this, arg) + class(RealMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Complex(kind=kind(0.D0)), allocatable, dimension(:,:) :: out + Real(kind=kind(0.D0)), allocatable, dimension(:) :: rwork + Integer :: i, j, sz1, sz2 + + end subroutine + + subroutine RealMat_lmult(this, arg) + class(RealMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Complex(kind=kind(0.D0)), allocatable, dimension(:,:) :: temp + + end subroutine + + subroutine RealMat_lmultinv(this, arg) + class(RealMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Complex(kind=kind(0.D0)), allocatable, dimension(:,:) :: temp + + end subroutine + + subroutine CmplxMat_init(this, arg) + class(CmplxMat) :: this + Complex(kind=kind(0.D0)), intent(inout), allocatable, dimension(:,:) :: arg + Integer :: i,j +! do i = 1, size(arg,1) +! write (*,*) (this%mat(i,j), j = 1, size(arg,2) ) +! enddo + this%m = size(arg, 1) + this%n = size(arg, 2) + if (allocated(this%mat)) deallocate(this%mat) + allocate(this%mat(this%m, this%n)) + do i = 1, this%m + do j = 1, this%n + this%mat(i, j) = arg(i, j) + enddo + enddo +! this%mat = arg !copy argument to local storage + + end subroutine + + subroutine CmplxMat_simt(this, arg) + class(CmplxMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), allocatable, dimension(:,:) :: arg + end subroutine + + subroutine CmplxMat_rmult(this, arg) + class(CmplxMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Complex(kind=kind(0.D0)), allocatable, dimension(:,:) :: out + Complex(kind=kind(0.d0)) :: alpha, zero + Integer :: i, j, sz1, sz2 + + alpha = 1.0 + zero = 0 + sz1 = size(arg, 1) + sz2 = size(arg, 2) + allocate(out(sz1, sz2)) + call zhemm('R', 'U', sz1, sz2, alpha, arg, sz1, this%mat, this%m, zero, out, sz1) + arg = out + deallocate(out) + end subroutine + + subroutine CmplxMat_rmultinv(this, arg) + class(CmplxMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + Complex(kind=kind(0.D0)), allocatable, dimension(:,:) :: out + Complex(kind=kind(0.d0)) :: alpha, zero + Integer :: i, j, sz1, sz2 + + end subroutine + + subroutine CmplxMat_lmult(this, arg) + class(CmplxMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + end subroutine + + subroutine CmplxMat_lmultinv(this, arg) + class(CmplxMat), intent(in) :: this + Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg + end subroutine +end module matTypes_mod diff --git a/testsuite/CMakeLists.txt b/testsuite/CMakeLists.txt index c245fcc48af0bc9de7f2947e566f3c3911e7783a..88c8381b7ab2950d5da347120689a67681381810 100644 --- a/testsuite/CMakeLists.txt +++ b/testsuite/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 2.8.5) +cmake_minimum_required(VERSION 2.8.11) Project(matmod.tests C Fortran) add_subdirectory(matmod.tests) diff --git a/testsuite/Prog.tests/26-Test-Polymorphic-Fortran.F90 b/testsuite/Prog.tests/26-Test-Polymorphic-Fortran.F90 new file mode 100644 index 0000000000000000000000000000000000000000..d7b29e1e8217160810379f40929ca8021e268cb6 --- /dev/null +++ b/testsuite/Prog.tests/26-Test-Polymorphic-Fortran.F90 @@ -0,0 +1,97 @@ +program test + Use DynamicMatrixArray_mod + Use ContainerElementBase_mod + Use OpTTypes_mod + Use Operator_mod + implicit none + + Type(DynamicMatrixArray) :: vec + Type(Operator), dimension(:), allocatable :: Op_T + Class(RealExpOpT), pointer :: reopt => null() + Class(CmplxExpOpT), pointer :: cmplxopt => null() + Class(ContainerElementBase), pointer :: dummy + Complex(kind=kind(0.d0)), allocatable, dimension(:,:) :: res + Complex(kind=kind(0.d0)) :: alpha, zero + Integer :: i, j, nmax, ndimmax + + nmax = 5 + ndimmax = 5 + allocate (res(ndimmax, ndimmax)) + call vec%init() + + alpha = 1.0 + zero = 0.0 + ! initialize res as unit matrix + call zlaset('A', ndimmax, ndimmax, zero, alpha, res, ndimmax) + + allocate(Op_T(2*nmax)) + + do i = 1, nmax + Call Op_make(Op_T(i), Ndimmax) + Call Op_make(Op_T(nmax + i), Ndimmax) + Op_T(i)%O = 0 + Op_T(nmax + i)%O = 0 + Do j = 1, ndimmax + Op_T(i)%P(j) = j + Op_T(nmax + i)%P(j) = j + Enddo + Op_T(i)%g = 0.2 + Op_T(i)%type = 2 + Op_T(i)%alpha = cmplx(0.d0, 0.d0, kind(0.D0)) + Op_T(nmax + i)%g = 0.2 + Op_T(nmax + i)%type = 2 + Op_T(nmax + i)%alpha = cmplx(0.d0, 0.d0, kind(0.D0)) + ! fill with some data + do j = 1, ndimmax + Op_T(i)%O(j, j) = 0!j + if (j+1 <= ndimmax) then + Op_T(nmax + i)%O(j, j+1) = cmplx(j, j, kind(0.D0)) + Op_T(nmax + i)%O(j+1, j) = cmplx(j, -j, kind(0.D0)) + endif + enddo + + Call Op_set(Op_T(i)) + Call Op_set(Op_T(nmax + i)) + + allocate(reopt, cmplxopt) + call reopt%init(Op_T(i)) + call cmplxopt%init(Op_T(nmax + i)) + call vec%pushback(reopt) + call vec%pushback(cmplxopt) + enddo + + + ! execute a loop over all stored objects + do i= 1, vec%length() + dummy => vec%at(i) ! get object + call dummy%rmult(res) ! polymorphic dispatch to rmult + call dummy%lmult(res) ! polymorphic dispatch to lmult +! do k = 1, ndimmax +! write (*,*) (aimag(res(k,l)), l = 1,ndimmax ) +! enddo +! write (*,*) "============" + enddo + +! do i = 1, ndimmax +! write (*,*) (res(i,j), j = 1,ndimmax ) +! enddo + if (abs(dble(res(ndimmax,ndimmax-1)) - 559995.58637168515D00) > 559995.58637168515D00*1D-15) then + write (*,*) "error in OpT mult.", abs(dble(res(ndimmax,ndimmax-1))-559995.58637168515D00), 559995.58637168515D00*1D-15 + stop 1 + endif + + if (abs(aimag(res(ndimmax,ndimmax-1)) + 559995.58637168526D00 ) > 559995.58637168526D00 * 1D-15) then ! ref is negative + write (*,*) "error in OpT multiplication",abs(aimag(res(ndimmax, ndimmax-1)) + 559995.58637168526D00 ) + stop 2 + endif +! write (*,*) res(ndimmax, ndimmax) + + ! tidy up + do i = 1, vec%length() + dummy => vec%at(i) ! Fortran doesn't want chaining + deallocate(dummy) + call Op_clear(Op_T(i), ndimmax) + enddo + call vec%dealloc() + deallocate(res, Op_T) +end program diff --git a/testsuite/Prog.tests/CMakeLists.txt b/testsuite/Prog.tests/CMakeLists.txt index f44604ffef7d7db5d648a731b6f9ad30929e1355..20c6938623a83eeb97c1760a500e96dfba5d07f4 100644 --- a/testsuite/Prog.tests/CMakeLists.txt +++ b/testsuite/Prog.tests/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 2.8.5) +cmake_minimum_required(VERSION 2.8.11) Project(matmod.tests C Fortran) find_package(LAPACK) @@ -28,8 +28,9 @@ add_executable(22-sl-mat-mults 21-sl-mat-mults.F90) add_executable(23-cgrp 23-cgrp.F90) add_executable(24-udv 24-udv.F90) add_executable(25-assign-UDV-state 25-assign-UDV-state.F90) +add_executable(26-Test-Polymorphic-Fortran 26-Test-Polymorphic-Fortran.F90) # set(operatorparts 1-copy-sel-rows 2-copy-sel-columns 3-opmult 4-opmultct 5-FillExpOps 6-opexpmult-1 6-opexpmult-2 7-opexpmultct-1 7-opexpmultct-2 9-Op-Phase 10-Op-mmultL 11-Op-mmultR 13-Op-Wrapup 14-Op-Wrapdo 21-Op-make 22-sl-mat-mults) -set(operatorparts 9-Op-Phase 10-Op-mmultL 11-Op-mmultR 13-Op-Wrapup 14-Op-Wrapdo 21-Op-make 22-sl-mat-mults) +set(operatorparts 9-Op-Phase 10-Op-mmultL 11-Op-mmultR 13-Op-Wrapup 14-Op-Wrapdo 21-Op-make 22-sl-mat-mults 26-Test-Polymorphic-Fortran) set(cgr2_2tests 16-get-blocks 17-solve-extended-system) #set(cgr2_1tests 8-scalematrix-1 8-scalematrix-2) set(cgr1tests 15-cgr 20-qdrp 23-cgrp) @@ -75,6 +76,8 @@ target_link_libraries(${iter}${CMAKE_CURRENT_SOURCE_DIR}/../../Libraries/libqrr target_link_libraries(${iter}${LAPACK_LIBRARIES}) endforeach(iter) +target_link_libraries(26-Test-Polymorphic-Fortran ${CMAKE_CURRENT_SOURCE_DIR}/../../Prog/OpTTypes_mod.o${CMAKE_CURRENT_SOURCE_DIR}/../../Prog/DynamicMatrixArray_mod.o \${CMAKE_CURRENT_SOURCE_DIR}/../../Prog/ContainerElementBase_mod.o) + enable_testing() # add_test(1-copy-sel-rows 1-copy-sel-rows) # add_test(2-copy-sel-columns 2-copy-sel-columns) @@ -103,3 +106,4 @@ add_test(22-sl-mat-mults 22-sl-mat-mults) add_test(23-cgrp 23-cgrp) add_test(24-udv 24-udv) add_test(25-assign-UDV-state 25-assign-UDV-state) +add_test(26-Test-Polymorphic-Fortran 26-Test-Polymorphic-Fortran) diff --git a/testsuite/matmod.tests/CMakeLists.txt b/testsuite/matmod.tests/CMakeLists.txt index 1a249602eb01c1f2c3164b6733f6e303fbb7519b..928ff4ca5d0babc283b5ee15057046619319a840 100644 --- a/testsuite/matmod.tests/CMakeLists.txt +++ b/testsuite/matmod.tests/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 2.8.5) +cmake_minimum_required(VERSION 2.8.11) Project(matmod.tests C Fortran) find_package(LAPACK)