Reverse Communication

[holysword]

Hi fellows,

I am working with an (parallel) iterative solver implementation right now, and I was told to implement reverse communication mechanism on it, to "unplug" the code from the data structures.

I checked some old codes containing this strategy and they are, in majority, FORTRAN77 codes. I have to admit that there were no better way to do that in F77, but I wonder if with current intel compiler it would be possible to use some object orientation to obtain the very same effect.

I would like to have advices from people who are already familiarized with this stragegy, and would like to know if there is a better approach with the current features of intel fortran compiler (our code is FORTRAN90/2003 nevertheless).

Thank you in advance

Holysword

[mecej4]

> I have to admit that there were no better way to do that in F77

From what basis do you draw that conclusion? The idea of reverse communication applies to any procedural programming language that has functions that can be called with arguments and return to the caller. There are advantages and disadvantages associated with the reverse communication mechanism, and you can find many articles about those on the Web.

Fortran 90 and later allows data to be shared by using modules and by host association, so you have more choices. Still, reverse communication remains a contender.

The current MKL library has a few solvers that exploit the reverse communication mechanism.

[holysword]

Hi mecej4, thank you for answering.

>From what basis do you draw that conclusion?

I was refering to the abstraction of the data structure. As far as I know, F77 doesn't have many features allowing that.

>There are advantages and disadvantages associated with the reverse communication mechanism, and you can find many articles about those on the Web.

Unfortunatelly I found only discussions related to F77 (including ARPACK package), and as I mentioned before, Iam just wondering if with F90/2003 features already implemented in intel fotran compiler we can abstract the data structure from the code without using Reverse Communication, and what would be the drawbacks of each strategy.

[holysword]

I have implemented a "dummy algorithm" using Reverse Communication. It doesn't do anything significative, just manipualtion of some vectors, and the matrices. I implemented two simple examples, one reverse communication interface which works with diagonal matrices and the other which works with common full matrices. Follows the code:

[fxfortran]MODULE rc_mod

 IMPLICIT NONE
 
 INTEGER, PARAMETER :: ST_INIT=0, ST_LOOP=1, ST_DONE=999, ST_AY_OUT=3, ST_MX_OUT=4, &
                       ST_AY_IN=-3, ST_MX_IN=-4
 
 CONTAINS

!--------------------------------------------------------------------
SUBROUTINE RC_Inner(x,y,i,max_i, status)
	
	!-------- PARAMETERS --------
	REAL, INTENT(INOUT), DIMENSION(:) :: x, y
	INTEGER, INTENT(INOUT) :: i, status
	INTEGER, INTENT(IN)    :: max_i
	
	!-------- LOCAL --------
	REAL :: alpha
	
	!-------- INIT --------
	IF (status.EQ.ST_INIT) THEN
		i = 1
		x = 1
		y = 0
		status = ST_LOOP
	ENDIF
		
!--------------------- MAIN LOOP ----------------------
	DO WHILE (i.LE.max_i)
		SELECT CASE (status)
			CASE (ST_LOOP)
				x = 2*x + y
				
				alpha = DOT_PRODUCT(x,x)
				alpha = SQRT(alpha)
				y = x / alpha

		! 		y = A*y
				status = ST_AY_OUT
				RETURN

			!-------------------
			CASE (ST_AY_IN)
				x = x - y
				alpha = DOT_PRODUCT(x,x)
				alpha = SQRT(alpha)
				x = x / alpha
				
		! 		y = M*x
				status = ST_MX_OUT
				RETURN
			!-------------------	
			CASE (ST_MX_IN)
				i = i + 1
				status = ST_LOOP
		END SELECT
	END DO

! Will reach here only after finishing the loop	
	status = ST_DONE
	RETURN

END SUBROUTINE RC_Inner
!--------------------------------------------------------------------


END MODULE rc_mod


MODULE rc_interface
 
 USE rc_mod
 
 IMPLICIT NONE
 
 CONTAINS

!--------------------------------------------------------------------
SUBROUTINE RC_Interface_Common(A, M, max_i)
	!-------- PARAMETER --------
	REAL, INTENT(IN), DIMENSION(:,:) :: A,M
	INTEGER, INTENT(IN) :: max_i
	
	!-------- LOCAL --------
	REAL, DIMENSION(:), ALLOCATABLE :: x, y
	INTEGER :: i, size_A, status
	
	
	size_A = SIZE(A,1) ! Since its, just a test, the matrix must be square
	
	ALLOCATE(x(size_A))
	ALLOCATE(y(size_A))
	
	status=ST_INIT
	CALL RC_Inner(x,y,i,max_i, status)
	
	DO
		SELECT CASE (status)
			!---------------------------
			CASE (ST_DONE)
				EXIT
				
			!---------------------------
			CASE (ST_AY_OUT)
				y = MATMUL(A,y)
				status = ST_AY_IN
				CALL RC_Inner(x,y,i,max_i, status)
				
			!---------------------------
			CASE (ST_MX_OUT)
				y = MATMUL(M,x)
				status = ST_MX_IN
				CALL RC_Inner(x,y,i,max_i, status)
				
		END SELECT
	END DO
	
	WRITE(*,*) '[RC_Interface_Common]:', x
	
END SUBROUTINE RC_Interface_Common


!--------------------------------------------------------------------
SUBROUTINE RC_Interface_Diag(A, M, max_i)
	!-------- PARAMETER --------
	REAL, INTENT(IN), DIMENSION(:) :: A,M
	INTEGER, INTENT(IN) :: max_i
	
	!-------- LOCAL --------
	REAL, DIMENSION(:), ALLOCATABLE :: x, y
	INTEGER :: i, j,size_A, status
	
	
	size_A = SIZE(A)
	
	ALLOCATE(x(size_A))
	ALLOCATE(y(size_A))
	
	status=ST_INIT
	CALL RC_Inner(x,y,i,max_i, status)
	
	DO
		SELECT CASE (status)
			!---------------------------
			CASE (ST_DONE)
				EXIT
				
			!---------------------------
			CASE (ST_AY_OUT)
				DO j=1,size_A
					y(j) = y(j)*A(j)
				END DO
				status = ST_AY_IN
				CALL RC_Inner(x,y,i,max_i, status)
				
			!---------------------------
			CASE (ST_MX_OUT)
				DO j=1,size_A
					y(j) = x(j)*M(j)
				END DO
				status = ST_MX_IN
				CALL RC_Inner(x,y,i,max_i, status)
				
		END SELECT
	END DO
	
	WRITE(*,*) '[RC_Interface_Diag]:', x
	
END SUBROUTINE RC_Interface_Diag
!--------------------------------------------------------------------


END MODULE rc_interface[/fxfortran]

Since its just an example I didn't put error manipulation and etc. I also had in mind only square matrices, and something would go wrong otherwise. I hope I didn't mistake anything in this piece of code.

[holysword]

... and then I implemented the very same algorithm, using object orientation ideas. Follows the code:

[fxfortran]MODULE oo_mtx

 IMPLICIT NONE

PUBLIC 

! ----------------------------------------------------------
! Declaration of the ABSTRACT CLASS
! ----------------------------------------------------------
TYPE, ABSTRACT ::  mtx
	! Just some cosmetics parameters
	INTEGER :: nnz, col, row
	CONTAINS
		PROCEDURE(matvec_abstract), DEFERRED :: matvec
END TYPE mtx 

INTERFACE
	SUBROUTINE matvec_abstract(self,x,y)
		IMPORT :: mtx
		CLASS(mtx),         INTENT(IN)  :: self
		REAL, DIMENSION(:), INTENT(IN)  :: x
		REAL, DIMENSION(:), INTENT(OUT) :: y
	END SUBROUTINE matvec_abstract
END INTERFACE

END MODULE oo_mtx



MODULE oo_mtx_diag

 USE oo_mtx

 IMPLICIT NONE
 
! ----------------------------------------------------------
! Declaration of the type mtx_diag
! ----------------------------------------------------------
TYPE, EXTENDS(mtx) :: mtx_diag
	REAL, DIMENSION(:),ALLOCATABLE :: A
	CONTAINS
	PROCEDURE, PASS(self) :: matvec
END TYPE mtx_diag
 
PUBLIC

PRIVATE :: matvec
 
 CONTAINS
 
 
! mtx_diag constructor
! ----------------------------------------------------------
FUNCTION mtx_diag_cnstr(col, row) RESULT(self)
	INTEGER, INTENT(IN) :: col, row
	TYPE(mtx_diag) :: self
	
	self%col  = col
	self%row  = row
	self%nnz = min(col,row)
	ALLOCATE( self%A( self%nnz ) )
	self%A = sqrt(-1.0)

END FUNCTION mtx_diag_cnstr
 

! mtx_diag matvec
! ----------------------------------------------------------
SUBROUTINE matvec(self,x,y)
	CLASS(mtx_diag),    INTENT(IN)  :: self
	REAL, DIMENSION(:), INTENT(IN)  :: x
	REAL, DIMENSION(:), INTENT(OUT) :: y
	
	INTEGER :: j
	
	DO j=1,self%nnz
		y(j) = x(j)*self%A(j)
	END DO
END SUBROUTINE matvec
 
 
END MODULE oo_mtx_diag



MODULE oo_mtx_common

 USE oo_mtx

 IMPLICIT NONE
 
! ----------------------------------------------------------
! Declaration of the type mtx_common
! ----------------------------------------------------------
TYPE, EXTENDS(mtx) :: mtx_common
	REAL, DIMENSION(:,:), ALLOCATABLE :: A
	CONTAINS
		PROCEDURE, PASS(self) :: matvec
END TYPE mtx_common

PUBLIC

PRIVATE :: matvec

 CONTAINS


! mtx_common constructor
! ----------------------------------------------------------
FUNCTION mtx_common_cnstr(col, row, nnz) RESULT(self)
	INTEGER, INTENT(IN) :: col, row, nnz
	TYPE(mtx_common) :: self
	
	self%col = col
	self%row = row
	self%nnz = nnz
	ALLOCATE( self%A( row, col ) )
	self%A = sqrt(-1.0)

END FUNCTION mtx_common_cnstr


! mtx_diag matvec
! ----------------------------------------------------------
SUBROUTINE matvec(self,x,y)
	CLASS(mtx_common),  INTENT(IN)  :: self
	REAL, DIMENSION(:), INTENT(IN)  :: x
	REAL, DIMENSION(:), INTENT(OUT) :: y
	
	y = MATMUL(self%A,x)
END SUBROUTINE matvec


END MODULE oo_mtx_common


MODULE oo_mod
 
 USE oo_mtx
 
 IMPLICIT NONE

 CONTAINS
 
!--------------------------------------------------------------------
! MAIN SUBROUTINE
!--------------------------------------------------------------------
SUBROUTINE OO_Interface(A, M, max_i)
	CLASS(mtx), INTENT(IN)  :: A, M
	INTEGER,  INTENT(IN)    :: max_i
	
	REAL, DIMENSION(A%col) :: x, y
	REAL :: alpha
	INTEGER :: i
	
	x = 1
	y = 0
		
!--------------------- MAIN LOOP ----------------------
	DO i=1,max_i
		x = 2*x + y
		alpha = DOT_PRODUCT(x,x)
		alpha = SQRT(alpha)
		y = x / alpha

		CALL A%matvec(y,y)

		x = x - y
		alpha = DOT_PRODUCT(x,x)
		alpha = SQRT(alpha)
		x = x / alpha
				
		CALL M%matvec(x,y)
		
	END DO
	
	WRITE(*,*) '[OO_Interface]:', x
	
END SUBROUTINE OO_Interface
!--------------------------------------------------------------------

END MODULE oo_mod
[/fxfortran]

[holysword]

So, answering my initial question: is it possible to implement abstraction of data type using other strategy than Reverse Communication? The answer is "yes", as I have just shown. I wrote a code for testing this implementation, which follows now:

[fxfortran]PROGRAM main

USE rc_interface
USE oo_mtx_diag
USE oo_mtx_common
USE oo_mod


IMPLICIT NONE


REAL, DIMENSION(5,5) :: A1, M1
REAL, DIMENSION(5)   :: A2, M2
TYPE(MTX_COMMON)     :: A3, M3
TYPE(MTX_DIAG)       :: A4, M4
INTEGER :: max_i

! Standard
!------------------------------------------------------
A1(1,:) = (/ 1.0, 0.0, 0.0, 0.0, 0.0/)
A1(2,:) = (/ 0.0,-2.3, 0.0, 0.0, 0.0/)
A1(3,:) = (/ 0.0, 0.0, 1.0, 0.0, 0.0/)
A1(4,:) = (/ 0.0, 0.0, 0.0,-1.4, 0.0/)
A1(5,:) = (/ 0.0, 0.0, 0.0, 0.0, 5.0/)


M1(1,:) = (/-1.0, 0.0, 0.0, 0.0, 0.0/)
M1(2,:) = (/ 0.0, 2.0, 0.0, 0.0, 0.0/)
M1(3,:) = (/ 0.0, 0.0, 1.2, 0.0, 0.0/)
M1(4,:) = (/ 0.0, 0.0, 0.0, 1.0, 0.0/)
M1(5,:) = (/ 0.0, 0.0, 0.0, 0.0, 2.9/)

! Diagonal
!------------------------------------------------------
A2 = (/ 1.0,-2.3, 1.0,-1.4, 5.0/)
M2 = (/-1.0, 2.0, 1.2, 1.0, 2.9/)


! Class Standard
!------------------------------------------------------
A3 = mtx_common_cnstr(5,5,5)
M3 = mtx_common_cnstr(5,5,5)
A3%A = A1
M3%A = M1

! Class Diagonal
!------------------------------------------------------
A4 = mtx_diag_cnstr(5,5)
M4 = mtx_diag_cnstr(5,5)
A4%A = A2
M4%A = M2

max_i = 1

CALL RC_Interface_Common(A1, M1, max_i)
CALL RC_Interface_Diag(A2, M2, max_i)

CALL OO_Interface(A3, M4, max_i)
CALL OO_Interface(A4, M3, max_i)

END PROGRAM[/fxfortran]

Once again, its just a test. I find interesting also is that OO_Interface can receive BOTH mtx_diag and mtx_common (that is, they can be mixed), whereas for Reverse Communication I would have to create a new RC_Interface_CommonDiag for doing that (and perhaps anotherRC_Interface_DiagCommon). I would dare to say that the object orientation strategy is more general and flexible. Another important point (which is a matter of personal opinion): I find it much easier to read, write and manipulate the object oriented version.

But the most important question would still remain. What are the disvantages regarding the object oriented version? Would it have any (considerable) performance problems? The fact that matvec is a DEFERRED procedure would make it impossible to inline the function at compiling time (right?), but can the Reverse Communication version be inlined without problems? If the MATMULT instriscic function is what takes most of the computational time (let's say, 98%), would the context exchange of calling matvec function decrease significatively the performance of the whole algorithm?