globals.F90

!!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)

Source Code

!===================================================================================================================================
! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics
! License: MIT
!===================================================================================================================================
#include "defines.FPP"

!===================================================================================================================================
!>
!!# Module **Globals**
!!
!! Here globally used variables /functions are defined
!!
!===================================================================================================================================
MODULE MODgvec_Globals

#ifndef NOISOENV
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY : INPUT_UNIT, OUTPUT_UNIT, ERROR_UNIT
#endif
USE_MPI

! USE MODgvec_py_abort, ONLY: PY_ABORT

IMPLICIT NONE

PUBLIC

!-----------------------------------------------------------------------------------------------------------------------------------
! Select here the working precision wp
!INTEGER, PARAMETER :: wp = selected_real_kind(6,35)   !! single precision
INTEGER, PARAMETER :: wp = selected_real_kind(15,307)  !! double precision
!INTEGER, PARAMETER :: wp = selected_real_kind(33,307) !! quadruple precision
!-----------------------------------------------------------------------------------------------------------------------------------
CHARACTER(LEN=20)   :: fmt_sep ='(132("="))'             !! formatting of separator line: WRITE(*,fmt_sep)
REAL(wp),PARAMETER  :: PI   =ACOS(-1.0_wp)               !! pi parameter
REAL(wp),PARAMETER  :: TWOPI=2.0_wp*PI                   !! 2*pi parameter
INTEGER             :: n_warnings_occured=0              !! for final line in screen output: 0 no warnings occured
!-----------------------------------------------------------------------------------------------------------------------------------
!for testing
LOGICAL                     :: testdbg=.FALSE.           !! for debugging the tests, set true for implementing tests, false to run
INTEGER                     :: testlevel =-1             !! flag for testing routines in code: -1: off
INTEGER                     :: ntestCalled=0             !! counter for called tests
INTEGER                     :: nfailedMsg=0              !! counter for messages on failed tests
INTEGER                     :: testUnit                  !! unit for out.test file
!MPI--------------------------------------------------------------------------------------------------------------------------------
LOGICAL                     :: MPIRoot=.TRUE.            !! flag whether process is MPI root process
INTEGER                     :: myRank=0                  !! rank of the MPI task
INTEGER                     :: nRanks=1                  !! total number of MPI tasks
!-----------------------------------------------------------------------------------------------------------------------------------
CHARACTER(LEN=20)           :: active_region(5)=(/"", "", "", "", ""/) !! for abort messages, to identify which (sub-)region was currently
INTEGER                     :: iregion=0                 !! which active_region to fill
INTEGER                     :: ProgressBar_oldpercent    !! for progressBar
REAL(wp)                    :: ProgressBar_starttime     !! for progressBar
!-----------------------------------------------------------------------------------------------------------------------------------
#ifndef NOISOENV
INTEGER, PARAMETER          :: UNIT_stdIn  = input_unit  !! Terminal input
INTEGER, PARAMETER          :: UNIT_stdOut = output_unit !! Terminal output
INTEGER, PARAMETER          :: UNIT_errOut = error_unit  !! For error output
#else
INTEGER, PARAMETER          :: UNIT_stdIn  = 5           !! Terminal input
INTEGER, PARAMETER          :: UNIT_stdOut = 6           !! Terminal output
INTEGER, PARAMETER          :: UNIT_errOut = 0           !! For error output
#endif
LOGICAL                     :: print_backtrace=.TRUE.  !! print backtrace on abort if compiled with GNU compiler
INTEGER, PARAMETER          :: MAXLEN  = 16384           !! max length of strings

INTERFACE reset_subregion
  MODULE PROCEDURE reset_subregion
END INTERFACE

INTERFACE enter_subregion
  MODULE PROCEDURE enter_subregion
END INTERFACE

INTERFACE exit_subregion
  MODULE PROCEDURE exit_subregion
END INTERFACE

INTERFACE Abort
   MODULE PROCEDURE Abort
END INTERFACE

ABSTRACT INTERFACE
  SUBROUTINE RaiseException(ErrorMessage)
    CHARACTER(LEN=*), INTENT(IN) :: ErrorMessage
  END SUBROUTINE RaiseException
END INTERFACE

PROCEDURE(RaiseException), POINTER :: RaiseExceptionPtr  => NULL()

INTERFACE GetTime
  MODULE PROCEDURE GetTime
END INTERFACE GetTime

INTERFACE GetTimeSerial
  MODULE PROCEDURE GetTimeSerial
END INTERFACE GetTimeSerial

INTERFACE ProgressBar
   MODULE PROCEDURE ProgressBar
END INTERFACE

INTERFACE GETFREEUNIT
   MODULE PROCEDURE GETFREEUNIT
END INTERFACE

INTERFACE Eval1DPoly
  MODULE PROCEDURE Eval1DPoly
END INTERFACE

INTERFACE CROSS
   MODULE PROCEDURE CROSS
END INTERFACE

INTERFACE NORMALIZE
   MODULE PROCEDURE NORMALIZE
END INTERFACE

INTERFACE DET33
   MODULE PROCEDURE DET33
END INTERFACE

INTERFACE INV33
   MODULE PROCEDURE INV33
END INTERFACE

CONTAINS


!==================================================================================================================================
!> reset global variables for the subregion output to default
!==================================================================================================================================
SUBROUTINE reset_subregion
  ! MODULES
  IMPLICIT NONE
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  iregion=0
  active_region=""
END SUBROUTINE reset_subregion

!==================================================================================================================================
!> add the current subregion to the active_regions (maximum depth is 5)
!! This information is collected uniquely for the abort error message
!!
!==================================================================================================================================
SUBROUTINE enter_subregion(subregion_name)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  CHARACTER(LEN=*), INTENT(IN) :: subregion_name
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
#if DEBUG
  CHARACTER(LEN=MAXLEN) :: regions
  INTEGER :: i
#endif
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  IF(MPIroot)THEN
    IF(iregion>4) CALL Abort(__STAMP__,&
                         "active subregion reached maximum depth of 5")
    iregion=iregion+1
    active_region(iregion)=subregion_name
#if DEBUG
    regions=active_region(1)
    DO i=2,iregion
      regions=TRIM(regions)//"."//TRIM(active_region(i))
    END DO
    SWRITE(Unit_stdOut,'(A)') '==> entering '//TRIM(regions)
#endif
  END IF
END SUBROUTINE enter_subregion

!==================================================================================================================================
!> remove the current subregion from the active subregions
!!
!==================================================================================================================================
SUBROUTINE exit_subregion(subregion_name)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  CHARACTER(LEN=*), INTENT(IN) :: subregion_name
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
#if DEBUG
  CHARACTER(LEN=MAXLEN) :: regions
  INTEGER :: i
#endif
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  IF(MPIroot)THEN
#if DEBUG
    regions=active_region(1)
    DO i=2,iregion
      regions=TRIM(regions)//"."//TRIM(active_region(i))
    END DO
    SWRITE(Unit_stdOut,'(A)') '<==  exiting '//TRIM(regions)
#endif
    IF(TRIM(subregion_name).NE.TRIM(active_region(iregion))) &
      CALL Abort(__STAMP__,&
                "trying to exit subregion '"//TRIM(subregion_name)// &
                "', but currently active subregion is '"//TRIM(active_region(iregion))//"'")
    active_region(iregion)=""
    iregion=iregion-1
  END IF
END SUBROUTINE exit_subregion



!==================================================================================================================================
!> Terminate program correctly if an error has occurred (important in MPI mode!).
!! Uses a MPI_ABORT which terminates FLUXO if a single proc calls this routine.
!!
!==================================================================================================================================
SUBROUTINE Abort(SourceFile,SourceLine,CompDate,CompTime,ErrorMessage,IntInfo,RealInfo,ErrorCode,TypeInfo)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  CHARACTER(LEN=*)                  :: SourceFile      !! Source file where error has occurred
  INTEGER                           :: SourceLine      !! Line in source file
  CHARACTER(LEN=*)                  :: CompDate        !! Compilation date
  CHARACTER(LEN=*)                  :: CompTime        !! Compilation time
  CHARACTER(LEN=*)                  :: ErrorMessage    !! Error message
  INTEGER,OPTIONAL                  :: IntInfo         !! additional integer value for error message
  REAL(wp),OPTIONAL                 :: RealInfo        !! additional real value for error message
  INTEGER,OPTIONAL                  :: ErrorCode       !! used for MPI
  CHARACTER(LEN=*),OPTIONAL         :: TypeInfo        !! Error type, default is "RuntimeError". Or e.g.
                                                       !! "MissingParameterError","InvalidParameterError","FileNotFoundError","InitializationError"
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  CHARACTER(LEN=50)                 :: IntString,RealString,errtype
#if MPI
  INTEGER                           :: errOut          ! Output of MPI_ABORT
  INTEGER                           :: signalout       ! Output errorcode
#endif
  CHARACTER(LEN=MAXLEN)             :: errmsg
  INTEGER                           :: i
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  IntString = ""
  RealString = ""
  errtype="RuntimeError"
  errmsg=""
  IF(MPIroot)THEN
    errmsg=TRIM(active_region(1))
    DO i=2,iregion
      errmsg=TRIM(errmsg)//"."//TRIM(active_region(i))
    END DO
    CALL reset_subregion()
  END IF
  IF(PRESENT(TypeInfo)) errtype = TRIM(TypeInfo)
  errmsg=TRIM(errmsg) // " | "//TRIM(errtype)

  IF (PRESENT(IntInfo))  THEN
    WRITE(IntString,"(I8)")  IntInfo
    IntString=",IntInfo="//TRIM(IntString)
  END IF
  IF (PRESENT(RealInfo)) THEN
     WRITE(RealString,"(F24.19)") RealInfo
     RealString=",RealInfo="//TRIM(RealString)
  END IF
  errmsg=TRIM(errmsg)//" | "//TRIM(ErrorMessage)//TRIM(IntString)//TRIM(RealString)

  WRITE(UNIT_stdOut,*) '_____________________________________________________________________________\n', &
                       'Program abort caused on Proc ',myRank, '\n', &
                       '  in File : ',TRIM(SourceFile),' Line ',SourceLine, '\n', &
                       '  This file was compiled at ',TRIM(CompDate),'  ',TRIM(CompTime), '\n', &
                       'Message: ',TRIM(errmsg)

  CALL FLUSH(UNIT_stdOut)

#if MPI
signalout=2 ! MPI_ABORT requires an output error-code /=0
IF(PRESENT(ErrorCode)) signalout=ErrorCode
CALL MPI_ABORT(MPI_COMM_WORLD,signalout,errOut)
#endif

#if GNU
  IF(print_backtrace) CALL BACKTRACE
#endif

  IF (ASSOCIATED(RaiseExceptionPtr)) THEN
    CALL RaiseExceptionPtr(errmsg)
  END IF
  ERROR STOP 2
END SUBROUTINE Abort

!==================================================================================================================================
!> Calculates current time (serial / OpenMP /MPI)
!!
!==================================================================================================================================
FUNCTION GetTime() RESULT(t)
  ! MODULES
!$ USE omp_lib
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  REAL(wp) :: t   !< output time
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
#if MPI
  LOGICAL :: barr
  INTEGER :: ierr
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)  ! not possible to 'CALL parBarrier()' because MODgvec_MPI uses MODgvec_Globals!
  t = MPI_WTIME()
#else
  CALL CPU_TIME(t)
!$ t=OMP_GET_WTIME()
#endif
END FUNCTION GetTime

!==================================================================================================================================
!> Calculates current time locally on a MPIrank (no MPI Barrier)
!!
!==================================================================================================================================
FUNCTION GetTimeSerial() RESULT(t)
  ! MODULES
!$ USE omp_lib
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  REAL(wp) :: t   !< output time
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  CALL CPU_TIME(t)
!$ t=OMP_GET_WTIME()
END FUNCTION GetTimeSerial

!==================================================================================================================================
!> Print a progress bar to screen, call either with init=T or init=F
!!
!==================================================================================================================================
SUBROUTINE ProgressBar(iter,n_iter)
  ! MODULES
!$ USE omp_lib
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  INTEGER,INTENT(IN) :: iter,n_iter  !! iter ranges from 1...n_iter
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  CHARACTER(LEN=8)  :: fmtstr
  INTEGER           :: newpercent
  REAL(wp)          :: endTime
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  IF(.NOT.MPIroot)RETURN
  IF(iter.LE.0)THEN !INIT
    ProgressBar_oldpercent=0
    ProgressBar_StartTime=GetTimeSerial()
    WRITE(UNIT_StdOut,'(4X,A,I8)') &
    '|       10%       20%       30%       40%       50%       60%       70%       80%       90%      100%| ... of ',n_iter
    WRITE(UNIT_StdOut,'(4X,A1)',ADVANCE='NO')'|'
    CALL FLUSH(UNIT_stdOut)
  ELSE
    newpercent=FLOOR(REAL(iter,wp)/REAL(n_iter,wp)*(100.0_wp+1.0e-12_wp))
    WRITE(fmtstr,'(I4)')newpercent-ProgressBar_oldpercent
    IF(newpercent-ProgressBar_oldpercent.GT.0)THEN
      WRITE(UNIT_StdOut,'('//TRIM(fmtstr)//'("."))',ADVANCE='NO')
      CALL FLUSH(UNIT_stdOut)
    END IF
    ProgressBar_oldPercent=newPercent
    IF(newpercent.EQ.100)THEN
      EndTime=GetTimeSerial()
      WRITE(Unit_stdOut,'(A3,F8.2,A)') '| [',EndTime-ProgressBar_StartTime,' sec ]'
    END IF
  END IF
END SUBROUTINE ProgressBar

!==================================================================================================================================
!> Get unused file unit number
!==================================================================================================================================
FUNCTION GETFREEUNIT()
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  INTEGER :: GetFreeUnit !! File unit number
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  LOGICAL :: connected
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  GetFreeUnit=55
  INQUIRE(UNIT=GetFreeUnit, OPENED=connected)
  IF(connected)THEN
    DO
      GetFreeUnit=GetFreeUnit+1
      INQUIRE(UNIT=GetFreeUnit, OPENED=connected)
      IF(.NOT.connected)EXIT
    END DO
  END IF
END FUNCTION GETFREEUNIT


!===================================================================================================================================
!> evalute monomial polynomial c_1+c_2*x+c_3*x^2 ...
!!
!===================================================================================================================================
PURE FUNCTION Eval1DPoly(nCoefs,Coefs,x)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  INTEGER,  INTENT(IN)  :: nCoefs                   !! number of coefficients
  REAL(wp), INTENT(IN)  :: Coefs(nCoefs)            !! coefficients
  REAL(wp), INTENT(IN)  :: x                        !! evaluation position
  REAL(wp)              :: Eval1DPoly
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  INTEGER               :: i
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  Eval1DPoly=0.
  DO i=nCoefs,1,-1
    Eval1DPoly=Eval1DPoly*x+Coefs(i)
  END DO

END FUNCTION Eval1DPoly


!===================================================================================================================================
!> evalute first derivative monomial polynomial (c_1+c_2*x+c_3*x^2) -> (c_2+2*c_3*x+3*c_4*x^2 ...
!!
!===================================================================================================================================
PURE FUNCTION Eval1DPoly_deriv(nCoefs,Coefs,x)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  INTEGER,  INTENT(IN)  :: nCoefs                   !! number of coefficients
  REAL(wp), INTENT(IN)  :: Coefs(nCoefs)            !! coefficients
  REAL(wp), INTENT(IN)  :: x                        !! evaluation position
  REAL(wp)              :: Eval1DPoly_deriv
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  INTEGER               :: i
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  Eval1DPoly_deriv=0.
  DO i=nCoefs,2,-1
    Eval1DPoly_deriv=Eval1DPoly_deriv*x+REAL(i-1,wp)*Coefs(i)
  END DO

END FUNCTION Eval1DPoly_deriv


!===================================================================================================================================
!> normalizes a nDim vector with respect to the eucledian norm
!!
!===================================================================================================================================
PURE FUNCTION NORMALIZE(v1,nVal)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  INTEGER,INTENT(IN)  :: nVal     !! vector size
  REAL(wp),INTENT(IN) :: v1(nVal) !! vector
  REAL(wp)            :: normalize(nVal) !! result, normalized vector
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  normalize=v1/SQRT(SUM(v1*v1))
END FUNCTION NORMALIZE


!===================================================================================================================================
!> computes the cross product of to 3 dimensional vectors: cross=v1 x v2
!!
!===================================================================================================================================
PURE FUNCTION CROSS(v1,v2)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  REAL(wp),INTENT(IN) :: v1(3) !! first input vector
  REAL(wp),INTENT(IN) :: v2(3) !! second input vector
  REAL(wp)            :: cross(3)  !! result v1 x v2
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  cross=(/v1(2)*v2(3)-v1(3)*v2(2),v1(3)*v2(1)-v1(1)*v2(3),v1(1)*v2(2)-v1(2)*v2(1)/)
END FUNCTION CROSS


!===================================================================================================================================
!> compute determinant of 3x3 matrix
!!
!===================================================================================================================================
PURE FUNCTION DET33(Mat)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  REAL(wp),INTENT(IN)  :: Mat(3,3) !! input matrix
  REAL(wp)             :: DET33 !! determinant of the input matrix
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  DET33=   ( Mat(1,1) * Mat(2,2) - Mat(1,2) * Mat(2,1) ) * Mat(3,3) &
           + ( Mat(1,2) * Mat(2,3) - Mat(1,3) * Mat(2,2) ) * Mat(3,1) &
           + ( Mat(1,3) * Mat(2,1) - Mat(1,1) * Mat(2,3) ) * Mat(3,2)
END FUNCTION DET33


!===================================================================================================================================
!> compute inverse of 3x3 matrix, needs sDet=1/det(Mat)
!!
!===================================================================================================================================
PURE FUNCTION INV33(Mat, Det_in)
  ! MODULES
  IMPLICIT NONE
  ! INPUT/OUTPUT VARIABLES ------------------------------------------------------------------------------------------------------!
  REAL(wp),INTENT(IN)             :: Mat(3,3) !! input matrix
  REAL(wp),INTENT(IN),OPTIONAL    ::  Det_in  !! determinant of input matrix (otherwise computed here)
  REAL(wp)                        :: INV33(3,3) !! inverse matrix
  ! LOCAL VARIABLES -------------------------------------------------------------------------------------------------------------!
  REAL(wp)                        :: sDet
  ! CODE ------------------------------------------------------------------------------------------------------------------------!
  IF(PRESENT(Det_in))THEN
    sDet=1.0_wp/Det_in
  ELSE
    sDet=1.0_wp/DET33(Mat)
  END IF
  INV33(1,1) = ( Mat(2,2) * Mat(3,3) - Mat(2,3) * Mat(3,2) ) * sDet
  INV33(1,2) = ( Mat(1,3) * Mat(3,2) - Mat(1,2) * Mat(3,3) ) * sDet
  INV33(1,3) = ( Mat(1,2) * Mat(2,3) - Mat(1,3) * Mat(2,2) ) * sDet
  INV33(2,1) = ( Mat(2,3) * Mat(3,1) - Mat(2,1) * Mat(3,3) ) * sDet
  INV33(2,2) = ( Mat(1,1) * Mat(3,3) - Mat(1,3) * Mat(3,1) ) * sDet
  INV33(2,3) = ( Mat(1,3) * Mat(2,1) - Mat(1,1) * Mat(2,3) ) * sDet
  INV33(3,1) = ( Mat(2,1) * Mat(3,2) - Mat(2,2) * Mat(3,1) ) * sDet
  INV33(3,2) = ( Mat(1,2) * Mat(3,1) - Mat(1,1) * Mat(3,2) ) * sDet
  INV33(3,3) = ( Mat(1,1) * Mat(2,2) - Mat(1,2) * Mat(2,1) ) * sDet

END FUNCTION INV33


END MODULE MODgvec_Globals

globals.F90 Source File

Contents

Modules

Source Code

Files dependent on this one

Source Code