!!matvec with matmul !#define __MATVEC_N(y,Mat,Vec) y=MATMUL(Mat,Vec) !#define __MATVEC_T(y,Mat,Vec) y=MATMUL(Vec,Mat) !#define __PMATVEC_N(fy,y,Mat,Vec) y=fyy+MATMUL(Mat,Vec) !#define __PMATVEC_T(fy,y,Mat,Vec) y=fyy+MATMUL(Vec,Mat) !#define __AMATVEC_N(y,fMat,Mat,Vec) y=fMatMATMUL(Mat,Vec) !#define __AMATVEC_T(y,fMat,Mat,Vec) y=fMatMATMUL(Vec,Mat) !#define __PAMATVEC_N(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Mat,Vec) !#define __PAMATVEC_T(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Vec,Mat)
!!#define __GENERICMATVEC(NT,fy,y,fMat,Mat,Vec) CALL DGEMV(NT,SIZE(Mat,1),SIZE(Mat,2),fMat,Mat,SIZE(Mat,1),Vec,1,fy,y,1)
!!matmat with matmul !#define __MATMAT_NN(Y,A,B) Y=MATMUL(A,B) !#define __MATMAT_TN(Y,A,B) Y=MATMUL(TRANSPOSE(A),B) !#define __MATMAT_NT(Y,A,B) Y=MATMUL(A,TRANSPOSE(B)) !#define __MATMAT_TT(Y,A,B) Y=TRANSPOSE(MATMUL(B,A))
!#define __PMATMAT_NN(fy,Y,A,B) Y=fyY+MATMUL(A,B) !#define __PMATMAT_TN(fy,Y,A,B) Y=fyY+MATMUL(TRANSPOSE(A),B) !#define __PMATMAT_NT(fy,Y,A,B) Y=fyY+MATMUL(A,TRANSPOSE(B)) !#define __PMATMAT_TT(fy,Y,A,B) Y=fyY+TRANSPOSE(MATMUL(B,A))
!#define __AMATMAT_NN(Y,fa,A,B) Y=faMATMUL(A,B) !#define __AMATMAT_TN(Y,fa,A,B) Y=faMATMUL(TRANSPOSE(A),B) !#define __AMATMAT_NT(Y,fa,A,B) Y=faMATMUL(A,TRANSPOSE(B)) !#define __AMATMAT_TT(Y,fa,A,B) Y=faTRANSPOSE(MATMUL(B,A))
!#define __PAMATMAT_NN(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,B) !#define __PAMATMAT_TN(fy,Y,fa,A,B) Y=fyY+faMATMUL(TRANSPOSE(A),B) !#define __PAMATMAT_NT(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,TRANSPOSE(B)) !#define __PAMATMAT_TT(fy,Y,fa,A,B) Y=fyY+faTRANSPOSE(MATMUL(B,A))
! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)
!#define __GENERICMATMAT_NN(fy,Y,fa,A,B) CALL DGEMM('N','N',SIZE(A,1),SIZE(B,2),SIZE(B,1),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TN(fy,Y,fa,A,B) CALL DGEMM('T','N',SIZE(A,2),SIZE(B,2),SIZE(B,1),fa,A,SIZE(B,1),B,SIZE(B,1),fy,Y,SIZE(A,2)) !#define __GENERICMATMAT_NT(fy,Y,fa,A,B) CALL DGEMM('N','T',SIZE(A,1),SIZE(B,1),SIZE(B,2),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TT(fy,Y,fa,A,B) CALL DGEMM('T','T',SIZE(A,2),SIZE(B,1),SIZE(B,2),fa,A,SIZE(B,2),B,SIZE(B,1),fy,Y,SIZE(A,2))
! SIMPLE INTERFACE FOR DGEMM, specifying nrows/ncols of mat A and nrows/ncols of mat B (for any transpose!) ! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)
!=================================================================================================================================== ! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics ! License: MIT !=================================================================================================================================== #include "defines.FPP" !=================================================================================================================================== !> !!# **GVEC** Driver program !! !=================================================================================================================================== PROGRAM TEST_GVEC_TO_HOPR USE MODgvec_Globals USE MODgvec_gvec_to_hopr IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- !local variables INTEGER :: i,nArgs,SFL CHARACTER(LEN=255) :: filename REAL(wp) :: StartTime,EndTime REAL(wp) :: xin(3,4),xout(3,4),data_out(9,4) REAL(wp) :: phi_edge_axis(2) REAL(wp) :: chi_edge_axis(2) !=================================================================================================================================== CALL CPU_TIME(StartTime) nArgs=COMMAND_ARGUMENT_COUNT() IF(nArgs.GE.1)THEN CALL GET_COMMAND_ARGUMENT(1,filename) ELSE STOP ' TEST GVEC TO HOPR: gvec filename not given, usage: "./executable gvec_file.dat"' END IF !header WRITE(Unit_stdOut,'(132("="))') WRITE(Unit_stdOut,'(5(("*",A128,2X,"*",:,"\n")))')& ' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '& ,' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '& ,' - - - - - - - - - - TEST GVEC => HOPR - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '& ,' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '& ,' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ' WRITE(Unit_stdOut,'(132("="))') DO SFL=0,2 !initialization phase CALL Init_gvec_to_hopr(filename,SFLcoord_in=SFL,factorSFL_in=2) WRITE(UNIT_stdOut,'(A,I4)')'===> SFLcoord: ',SFL xin(:,1)=(/0.0,0.5,0.3/) xin(:,2)=(/0.3,0.13,0.65/) xin(:,3)=(/0.6,0.43,0.15/) xin(:,4)=(/1.0,-0.33,-0.45/) CALL gvec_to_hopr(4,xin,xout,data_out,phi_edge_axis,chi_edge_axis) WRITE(UNIT_stdOut,'(A,2E21.13)')'phi_edge_axis: ',phi_edge_axis WRITE(UNIT_stdOut,'(A,2E21.13)')'chi_edge_axis: ',chi_edge_axis DO i=1,4 WRITE(UNIT_stdOut,'(A,3E21.13)')'s,thet,zeta: ',xin(:,i) WRITE(UNIT_stdOut,'(A,3E21.13)')'x,y,z : ',xout(:,i) WRITE(UNIT_stdOut,'(A, E21.13)')'pressure : ',data_out(1,i) WRITE(UNIT_stdOut,'(A,3E21.13)')'Bcart : ',data_out(2:4,i) WRITE(UNIT_stdOut,'(A, E21.13)')'|B| : ',SQRT(SUM(data_out(2:4,i)**2)) WRITE(UNIT_stdOut,'(A,2E21.13)')'chi,phi : ',data_out(5:6,i) WRITE(UNIT_stdOut,'(A,3E21.13)')'Acart : ',data_out(7:9,i) WRITE(UNIT_stdOut,*)'-----------------------' END DO CALL Finalize_gvec_to_hopr() END DO CALL CPU_TIME(EndTime) WRITE(Unit_stdOut,fmt_sep) WRITE(Unit_stdOut,'(A,F8.2,A)') ' TEST GVEC TO HOPR FINISHED! [',EndTime-StartTime,' sec ]' WRITE(Unit_stdOut,fmt_sep) END PROGRAM TEST_GVEC_TO_HOPR