!!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" !=================================================================================================================================== !> !!# Module **Linear Algebra** !! !! Provides the linear algebra wrapper routines using LAPACK. !! !!- matrix inverse !!- solve linear system !! !=================================================================================================================================== MODULE MODgvec_LinAlg USE MODgvec_Globals, ONLY:wp,abort IMPLICIT NONE PUBLIC CONTAINS !=================================================================================================================================== !> Computes matrix inverse using LAPACK !! Input matrix should be a square matrix !! !=================================================================================================================================== FUNCTION INV(A) RESULT(Ainv) ! MODULES IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT/OUTPUT VARIABLES REAL(wp),INTENT(IN) :: A(:,:) !! input matrix REAL(wp) :: Ainv(SIZE(A,1),SIZE(A,2)) !! result: inverse of A !----------------------------------------------------------------------------------------------------------------------------------- ! External procedures defined in LAPACK EXTERNAL DGETRF EXTERNAL DGETRI ! LOCAL VARIABLES REAL(wp):: work(SIZE(A,1)) ! work array for LAPACK INTEGER :: ipiv(SIZE(A,1)) ! pivot indices INTEGER :: n,info !=================================================================================================================================== ! Store A in Ainv to prevent it from being overwritten by LAPACK Ainv = A n = size(A,1) ! DGETRF computes an LU factorization of a general M-by-N matrix A ! using partial pivoting with row interchanges. CALL DGETRF(n, n, Ainv, n, ipiv, info) IF(info.NE.0)THEN CALL abort(__STAMP__,& 'Matrix is numerically singular!') END IF ! DGETRI computes the inverse of a matrix using the LU factorization ! computed by DGETRF. CALL DGETRI(n, Ainv, n, ipiv, work, n, info) IF(info.NE.0)THEN CALL abort(__STAMP__,& 'Matrix inversion failed!') END IF END FUNCTION INV !=================================================================================================================================== !> Solve linear system of dimension dims and multiple RHS !! !=================================================================================================================================== FUNCTION SOLVE(A,RHS) RESULT(X) ! MODULES IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES REAL(wp),INTENT(IN) :: A(:,:) !! matrix REAL(wp),INTENT(IN) :: RHS(:) !! RHS, sorting: (dimA,nRHS), two dimensions can be used in input !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES REAL(wp) :: X(SIZE(RHS,1)) !! result: solution of A X=RHS !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES ! External procedures defined in LAPACK EXTERNAL DGETRF EXTERNAL DGETRS ! LOCAL VARIABLES REAL(wp) :: Atmp(SIZE(A,1), SIZE(A,1)) INTEGER :: ipiv(SIZE(A,1)) ! pivot indices INTEGER :: nRHS,n,info !=================================================================================================================================== Atmp=A X = RHS n = SIZE(A,1) nRHS=SIZE(RHS,1)/SIZE(A,1) CALL DGETRF(n, n, Atmp, n, ipiv, info) IF(info.NE.0)THEN CALL abort(__STAMP__,& 'Matrix is numerically singular!') END IF CALL DGETRS('N',n, nRHS,Atmp, n, ipiv,X,n, info) IF(info.NE.0)THEN CALL abort(__STAMP__,& 'Matrix solve does not work!') END IF END FUNCTION SOLVE !=================================================================================================================================== !> Solve linear system of dimension dims and multiple RHS !! !=================================================================================================================================== FUNCTION SOLVEMAT(A,RHS) RESULT(X) ! MODULES IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES REAL(wp),INTENT(IN) :: A(:,:) !! matrix REAL(wp),INTENT(IN) :: RHS(:,:) !! RHS, sorting: (dimA,nRHS), two dimensions can be used in input !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES REAL(wp) :: X(SIZE(A,1),SIZE(RHS,2)) !! result: solution of A X=RHS !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES ! External procedures defined in LAPACK EXTERNAL DGETRF EXTERNAL DGETRS ! LOCAL VARIABLES REAL(wp) :: Atmp(SIZE(A,1), SIZE(A,1)) INTEGER :: ipiv(SIZE(A,1)) ! pivot indices INTEGER :: nRHS,n,info !=================================================================================================================================== Atmp=A X = RHS n = SIZE(A,1) nRHS=SIZE(RHS,2) CALL DGETRF(n, n, Atmp, n, ipiv, info) IF(info.NE.0)THEN CALL abort(__STAMP__,& 'Matrix is numerically singular!') END IF CALL DGETRS('N',n, nRHS,Atmp, n, ipiv,X,n, info) IF(info.NE.0)THEN CALL abort(__STAMP__,& 'Matrix solve does not work!') END IF END FUNCTION SOLVEMAT !=================================================================================================================================== !> Return P L U matrices of the LU decomposition, cmoputed from LAPACK Routine (if P is not passed, L=P*L) !! !=================================================================================================================================== SUBROUTINE getLU(dimA,A,L,U,P) ! MODULES IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES INTEGER, INTENT(IN) :: dimA REAL(wp),INTENT(IN) :: A(1:dimA,1:dimA) !! matrix !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES REAL(wp),INTENT(OUT) :: L(1:dimA,1:dimA) !! L or P*L (if P omitted) REAL(wp),INTENT(OUT) :: U(1:dimA,1:dimA) REAL(wp),INTENT(OUT),OPTIONAL :: P(1:dimA,1:dimA) !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES ! External procedures defined in LAPACK EXTERNAL DGETRF ! LOCAL VARIABLES REAL(wp) :: tmp,Atmp(dimA,dimA) INTEGER :: ipiv(dimA) ! pivot indices INTEGER :: i,j,info !=================================================================================================================================== Atmp=A CALL DGETRF(dimA, dimA, Atmp, dimA, ipiv, info) IF(info.NE.0)THEN CALL abort(__STAMP__,& 'Matrix is numerically singular!') END IF !upper diagonal part U=0. DO i=1,dimA DO j=i,dimA U(i,j)=Atmp(i,j) END DO END DO !lower diagonal part, known: 1 on the diagonal L=0. DO i=1,dimA L(i,i)=1. DO j=1,i-1 L(i,j)=Atmp(i,j) END DO END DO IF(PRESENT(P))THEN P=0. DO i=1,dimA P(i,i)=1. END DO !pivoting of rows in L is pivoting of columns in P DO i=1,dimA DO j=1,dimA tmp=P(j,i) P(j,i)=P(j,ipiv(i)) P(j,ipiv(i))=tmp END DO END DO !CHECK IF(SUM(ABS(MATMUL(P,MATMUL(L,U))-A)).GT.1e-12*dimA*dimA) THEN CALL abort(__STAMP__,& 'A=P*L*U decomposition not correct') END IF ELSE !pivoting of rows in L, backwards DO i=dimA,1,-1 DO j=1,dimA tmp=L(i,j) L(i,j)=L(ipiv(i),j) L(ipiv(i),j)=tmp END DO END DO !CHECK IF(SUM(ABS(MATMUL(L,U)-A)).GT.1e-12*dimA*dimA) THEN CALL abort(__STAMP__,& 'A=L*U decomposition not correct') END IF END IF END SUBROUTINE getLU END MODULE MODgvec_LinAlg