Fortran function to overload multiplication between derived types with allocatable components

Question

Foreword

In order to store banded matrices whose full counterparts can have both rows and columns indexed from indices other than 1, I defined a derived data type as

TYPE CDS
  REAL, DIMENSION(:,:), ALLOCATABLE :: matrix
  INTEGER, DIMENSION(2) :: lb, ub
  INTEGER :: ld, ud
END TYPE CDS

where CDS stands for compressed diagonal storage. Given the declaration TYPE(CDS) :: A,

The rank-2 component matrix is supposed to contain, as columns, the diagonals of the actual full matrix (like here, except that I store the diagonals as columns and not as rows).
The components ld and ud are supposed to contain the number of lower and upper diagonals respectively, that is -lbound(A%matrix,2) and +ubound(A%matrix,2).
The 2-elements components lb and ub are supposed to contain the lower bounds and upper bounds of the actual full matrix along the two dimensions. In particular lb(1) and ub(1) should be the same as lbound(A%matrix,1) and lbound(A%matrix,2).

As you can see in points 2. and 3., the derived type contains some redundant information, but I don't care, since they are just 3 couples of integers. Furthermore, in the code that I'm writing, the information about the bounds and the band of the actual full matrix is know before the matrix can be filled. So I first assign the values to the components ld, ud, lb and ub, and then I used these components to ALLOCATE the matrix component (then I can fill it properly).

The problem

I have to perform matrix multiplication between such sparse matrices, so I wrote a FUNCTION to perform such product and used it to overload the * operator.

At the moment the function is as follows,

FUNCTION CDS_mat_x_CDS_mat(A, B)
IMPLICIT NONE
TYPE(CDS), INTENT(IN) :: A, B
TYPE(CDS) :: cds_mat_x_cds_mat

! determine the lower and upper bounds and band of the result based on those of the operands
CDS_mat_x_CDS_mat%lb(1) = A%lb(1)
CDS_mat_x_CDS_mat%ub(1) = A%ub(1)
CDS_mat_x_CDS_mat%lb(2) = B%lb(2)
CDS_mat_x_CDS_mat%ub(2) = B%ub(2)
CDS_mat_x_CDS_mat%ld = A%ld + B%ld
CDS_mat_x_CDS_mat%ud = A%ud + B%ud

! allocate the matrix component
ALLOCATE(CDS_mat_x_CDS_mat%matrix(CDS_mat_x_CDS_mat%lb(1):CDS_mat_x_CDS_mat%ub(1),&
                               & -CDS_mat_x_CDS_mat%ld:+CDS_mat_x_CDS_mat%ud))

! perform the product
:
:

END FUNCTION

This means that, if I have to do the product multiple times the allocation is done multiple times inside the function. I think this is not good from a performance point of view.

I ask for suggestions on how to accomplish the task of banded sparse matrix times banded sparse matrix. I would like to use the type the I defined because I need it to be as general, in terms of bounds, as it is at the moment. But I could change the procedure to perform the product (from FUNCTION to SUBROUTINE, if necessary).

Ideas

I could rewrite the procedure as a SUBROUTINE in order to declare CDS_mat_x_CDS_mat with INTENT(INOUT) do the assignment of components other than matrix, as well as the allocation, outside the SUBROUTINE. The drawback would be that I could not overload the * operator.

I noticed the intrinsic function MATMUL can operate on any rank-2 operands, whichever the upper and lower bounds along the two dimensions. This means that the allocation is performed inside the function. I suppose it's efficient (since it is an intrinsic). The difference with respect to my function is that it accepts rank-2 arrays of any shape, wheres mine accept derived data types objects with a rank-2 array component of any shape.

Have you profiled to determine the (relative) cost of that allocation (which seems you have to do exactly once, whether in the function or before entering the subroutine)? But then, perhaps I don't understand: you suggest re-using the allocation when multiplying more than once. But why not just cache the result and have done? — francescalus
– francescalus, Commented Jul 6, 2016 at 17:42

IanH · Accepted Answer · 2016-07-06 20:52:59Z

2

The intrinsic function MATMUL has the equivalent of an automatic (F2008 5.2.2) result - the shape of the result is expressed in such a way that it becomes a characteristic of the function (F2008 12.3.3) - the shape of the function result is determined in the specification part of the function and (in terms of implementation) the compiler therefore knows how to calculate the shape of the function result ahead of actually executing the function proper.

As a consequence there are no Fortran language ALLOCATABLE variables associated with the equivalent of the intrinsic MATMUL function result. This is not the same thing as there being "no memory allocation" - the compiler may still need to allocate memory behind the scenes for its own purposes - for things like expressions temporaries and the like.

(I say "equivalent of" above, because intrinsic procedures are inherently special - but pretend for a moment that MATMUL was just a user function.)

The equivalent of that sort of automatic result for your case can be achieved by using length type parameters. This is Fortran 2003 feature - the same base language standard that introduced allocatable components - but it is not one that has yet been implemented by all actively maintained compilers.

MODULE xyz
  IMPLICIT NONE

  TYPE CDS(lb1, ub1, ld, ud)
    INTEGER, LEN :: lb1, ub1, ld, ud
    REAL :: matrix(lb1:ub1, ld:ud)
    INTEGER :: lb2, ub2
  END TYPE CDS

  INTERFACE OPERATOR(*)
    MODULE PROCEDURE CDS_mat_x_CDS_mat
  END INTERFACE OPERATOR(*)
CONTAINS
  FUNCTION CDS_mat_x_CDS_mat(A, B) RESULT(C)
    TYPE(CDS(*,*,*,*)), INTENT(IN) :: A, B
    TYPE(CDS(A%lb1, A%ub1, A%ld+B%ld, A%ud+B%ud)) :: C

    C%lb2 = B%lb2
    C%ub2 = B%ub2

    ! perform the product.
    ! :

  END FUNCTION CDS_mat_x_CDS_mat
END MODULE xyz

Theoretically this gives the compiler more opportunity for optimisation, because it has more insight ahead of the call to the function for the storage required for the function result. Whether this actually results in better real world performance depends on the implementation of the compiler and the nature of the function references.

answered Jul 6, 2016 at 20:52

IanH

21.6k2 gold badges41 silver badges61 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

Enlico Over a year ago

This seems just what I was looking for. Just two question about it: (1) I read a classic example of parametrized derived data types on Chapman's Fortran95/2003, but I see non LEN after the declaration of the parameteres. Is there a difference? (2) Can't I declare rank-1 2-elements lb and ub and use lb(1)use ub(1) in the declaration of the matrix component?

IanH Over a year ago

Kind parameters for user defined derived types are either KIND parameters - which have to be specified by compile time constants (and which can then be used in situations which require a compile time constant) - or LEN parameters, which can be specified at "run time". Type parameters are always scalar.

IanH Over a year ago

That should have been "Type parameters for user defined..."

IanH Over a year ago

If I understand what you are asking - then no - you cannot. The array spec for a non-allocatable, non-pointer component array cannot directly depend on the value of a variable - that is why the example in the answer uses a length type parameter. Type parameters are always scalar - they are never arrays.

Enlico Over a year ago

Oh, yes, I erroneously read your last words "Type [I've read Kind] parameters are always scalar". Ok! I'll try with your answer.

Collectives™ on Stack Overflow

Fortran function to overload multiplication between derived types with allocatable components

1 Answer 1

5 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related