Merge pull request #363 from easifem-fortran/axpy

Updates AXPY in AbstractNodeField_Class
easifem · Jan 2, 2024 · 5074180 · 5074180
2 parents 904dd78 + bc1b49c
commit 5074180
Show file tree

Hide file tree

Showing 2 changed files with 189 additions and 9 deletions.
diff --git a/src/modules/AbstractNodeField/src/AbstractNodeField_Class.F90 b/src/modules/AbstractNodeField/src/AbstractNodeField_Class.F90
@@ -157,8 +157,11 @@ MODULE AbstractNodeField_Class
 
   ! GET:
   ! @BlasMethods
-  PROCEDURE, PUBLIC, PASS(obj) :: AXPY => obj_AXPY
-  !! Y = Y + scale * X
+  PROCEDURE, PASS(obj) :: AXPY1 => obj_AXPY1
+  PROCEDURE, PASS(obj) :: AXPY2 => obj_AXPY2
+  PROCEDURE, PASS(obj) :: AXPY3 => obj_AXPY3
+  GENERIC, PUBLIC :: AXPY => AXPY1, AXPY2, AXPY3
+  !! Y = Y + scale * X ...
   PROCEDURE, PUBLIC, PASS(obj) :: SCAL => obj_SCAL
   !! X = scale * X
   PROCEDURE, PUBLIC, PASS(obj) :: COPY => obj_Copy
@@ -651,20 +654,58 @@ END SUBROUTINE obj_ApplyDirichletBC2
 ! summary:  y = y + s * x
 
 INTERFACE
-  MODULE SUBROUTINE obj_AXPY(obj, x, scale)
+  MODULE SUBROUTINE obj_AXPY1(obj, x, scale)
     CLASS(AbstractNodeField_), INTENT(INOUT) :: obj
     CLASS(AbstractNodeField_), INTENT(INOUT) :: x
     REAL(DFP), INTENT(IN) :: scale
-  END SUBROUTINE obj_AXPY
+  END SUBROUTINE obj_AXPY1
 END INTERFACE
 
 !----------------------------------------------------------------------------
 !                                                          AXPY@BlasMethods
 !----------------------------------------------------------------------------
 
+!> author: Shion Shimizu
+! date:   2023-12-29
+! summary:  y = y + a1 * x1 + a2 * x2
+
+INTERFACE
+  MODULE SUBROUTINE obj_AXPY2(obj, x1, x2, a1, a2)
+    CLASS(AbstractNodeField_), INTENT(INOUT) :: obj
+    CLASS(AbstractNodeField_), INTENT(INOUT) :: x1
+    CLASS(AbstractNodeField_), INTENT(INOUT) :: x2
+    REAL(DFP), INTENT(IN) :: a1
+    REAL(DFP), INTENT(IN) :: a2
+  END SUBROUTINE obj_AXPY2
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                          AXPY@BlasMethods
+!----------------------------------------------------------------------------
+
+!> author: Shion Shimizu
+! date:   2023-12-29
+! summary:  y = y + a1 * x1 + a2 * x2 + a3 * x3
+
+INTERFACE
+  MODULE SUBROUTINE obj_AXPY3(obj, x1, x2, x3, a1, a2, a3)
+    CLASS(AbstractNodeField_), INTENT(INOUT) :: obj
+    CLASS(AbstractNodeField_), INTENT(INOUT) :: x1
+    CLASS(AbstractNodeField_), INTENT(INOUT) :: x2
+    CLASS(AbstractNodeField_), INTENT(INOUT) :: x3
+    REAL(DFP), INTENT(IN) :: a1
+    REAL(DFP), INTENT(IN) :: a2
+    REAL(DFP), INTENT(IN) :: a3
+  END SUBROUTINE obj_AXPY3
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                          ACAL@BlasMethods
+!----------------------------------------------------------------------------
+
 !> author: Vikas Sharma, Ph. D.
 ! date:  2023-12-17
-! summary:  y = y + s * x
+! summary: scaling y = s * y
 
 INTERFACE
   MODULE SUBROUTINE obj_SCAL(obj, scale)

diff --git a/src/submodules/AbstractNodeField/src/[email protected] b/src/submodules/AbstractNodeField/src/[email protected]
@@ -24,11 +24,11 @@
 CONTAINS
 
 !----------------------------------------------------------------------------
-!                                                                     AXPY
+!                                                                     AXPY1
 !----------------------------------------------------------------------------
 
-MODULE PROCEDURE obj_AXPY
-CHARACTER(*), PARAMETER :: myName = "obj_AXPY()"
+MODULE PROCEDURE obj_AXPY1
+CHARACTER(*), PARAMETER :: myName = "obj_AXPY1()"
 LOGICAL(LGT) :: problem
 INTEGER(I4B) :: ierr
 
@@ -83,7 +83,146 @@
 CALL e%RaiseInformation(modName//'::'//myName//' - '// &
   & '[END]')
 #endif
-END PROCEDURE obj_AXPY
+END PROCEDURE obj_AXPY1
+
+!----------------------------------------------------------------------------
+!                                                                   AXPY2
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE obj_AXPY2
+CHARACTER(*), PARAMETER :: myName = "obj_AXPY2()"
+LOGICAL(LGT) :: problem
+INTEGER(I4B) :: ierr
+
+#ifdef DEBUG_VER
+CALL e%RaiseInformation(modName//'::'//myName//' - '// &
+  & '[START]')
+#endif
+
+SELECT CASE (obj%engine%chars())
+CASE (TypeEngineName%native_serial)
+
+#ifdef DEBUG_VER
+  problem = x1%engine .NE. "NATIVE_SERIAL" .OR.  &
+          & x2%engine .NE. "NATIVE_SERIAL"
+  IF (problem) THEN
+    CALL e%RaiseError(modName//'::'//myName//' - '// &
+      & '[INTERNAL ERROR] :: engine of x1 and x2'// &
+      & 'should be NATIVE_SERIAL.')
+    RETURN
+  END IF
+#endif
+
+  CALL AXPY(X=x1%realvec, Y=obj%realvec, A=a1)
+  CALL AXPY(X=x2%realvec, Y=obj%realvec, A=a2)
+
+#ifdef USE_LIS
+CASE (TypeEngineName%lis_omp)
+
+#ifdef DEBUG_VER
+  CALL lis_vector_is_null(obj%lis_ptr, ierr)
+  CALL CHKERR(ierr)
+  CALL lis_vector_is_null(x1%lis_ptr, ierr)
+  CALL CHKERR(ierr)
+  CALL lis_vector_is_null(x2%lis_ptr, ierr)
+  CALL CHKERR(ierr)
+#endif
+
+  ! alpha, x, y, ierr
+  CALL lis_vector_axpy(a1, x1%lis_ptr, obj%lis_ptr, ierr)
+  CALL lis_vector_axpy(a2, x2%lis_ptr, obj%lis_ptr, ierr)
+
+#ifdef DEBUG_VER
+  CALL CHKERR(ierr)
+#endif
+
+#endif
+
+CASE DEFAULT
+
+  CALL e%RaiseError(modName//'::'//myName//' - '// &
+    & '[INTERNAL ERROR] :: No case found given engine = '//  &
+    & obj%engine%chars())
+  RETURN
+END SELECT
+
+#ifdef DEBUG_VER
+CALL e%RaiseInformation(modName//'::'//myName//' - '// &
+  & '[END]')
+#endif
+END PROCEDURE obj_AXPY2
+
+!----------------------------------------------------------------------------
+!                                                                   AXPY3
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE obj_AXPY3
+CHARACTER(*), PARAMETER :: myName = "obj_AXPY3()"
+LOGICAL(LGT) :: problem
+INTEGER(I4B) :: ierr
+
+#ifdef DEBUG_VER
+CALL e%RaiseInformation(modName//'::'//myName//' - '// &
+  & '[START]')
+#endif
+
+SELECT CASE (obj%engine%chars())
+CASE (TypeEngineName%native_serial)
+
+#ifdef DEBUG_VER
+  problem = x1%engine .NE. "NATIVE_SERIAL" .OR.  &
+          & x2%engine .NE. "NATIVE_SERIAL" .OR. &
+          & x3%engine .NE. "NATIVE_SERIAL"
+  IF (problem) THEN
+    CALL e%RaiseError(modName//'::'//myName//' - '// &
+      & '[INTERNAL ERROR] :: engine of x1, x2 and x3'// &
+      & 'should be NATIVE_SERIAL.')
+    RETURN
+  END IF
+#endif
+
+  CALL AXPY(X=x1%realvec, Y=obj%realvec, A=a1)
+  CALL AXPY(X=x2%realvec, Y=obj%realvec, A=a2)
+  CALL AXPY(X=x3%realvec, Y=obj%realvec, A=a3)
+
+#ifdef USE_LIS
+CASE (TypeEngineName%lis_omp)
+
+#ifdef DEBUG_VER
+  CALL lis_vector_is_null(obj%lis_ptr, ierr)
+  CALL CHKERR(ierr)
+  CALL lis_vector_is_null(x1%lis_ptr, ierr)
+  CALL CHKERR(ierr)
+  CALL lis_vector_is_null(x2%lis_ptr, ierr)
+  CALL CHKERR(ierr)
+  CALL lis_vector_is_null(x3%lis_ptr, ierr)
+  CALL CHKERR(ierr)
+#endif
+
+  ! alpha, x, y, ierr
+  CALL lis_vector_axpy(a1, x1%lis_ptr, obj%lis_ptr, ierr)
+  CALL lis_vector_axpy(a2, x2%lis_ptr, obj%lis_ptr, ierr)
+  CALL lis_vector_axpy(a3, x3%lis_ptr, obj%lis_ptr, ierr)
+
+#ifdef DEBUG_VER
+  CALL CHKERR(ierr)
+#endif
+
+#endif
+
+CASE DEFAULT
+
+  CALL e%RaiseError(modName//'::'//myName//' - '// &
+    & '[INTERNAL ERROR] :: No case found given engine = '//  &
+    & obj%engine%chars())
+  RETURN
+END SELECT
+
+#ifdef DEBUG_VER
+CALL e%RaiseInformation(modName//'::'//myName//' - '// &
+  & '[END]')
+#endif
+END PROCEDURE obj_AXPY3
 
 !----------------------------------------------------------------------------
 !                                                                     SCAL