overload MATMUL(a,b) as a*b

[nooj]

This has got to be in the FAQ, but I didn't see it in any forum search or the googles.

I'd like to use the notation "C = A*B" to mean "C = MATMUL(A,B)". If I write a simple operator interface,

[bash]interface operator(*)
  intrinsic matmul
end interface operator [/bash]

I get "This statement is incorrectly positioned" because I don't know what I'm doing. If I write a more complicated interface that actually works, I lose much of the beauty and versatility of the builtin MATMUL interface. Reading reference manual entries for the "interface" statement left me confused.

Is there a good way to say simply "A*B" means "MATMUL(A,B)" (without knowing the size or type of A or B, for instance)? I'd also like "A*B*C" to mean "MATMUL(MATMUL(A,B),C)" as one might expect.

Overloading DOTPRPDUCT in a similar fashion (perhaps not using the same symbol) would also be nice.

- Nooj

[david_car · Accepted solution]

Nooj,

To be able to override operators you need to be define a derived type that you can override the * operator on. Heres an example:

[bash]module class_Matrix

  type Matrix
    real :: data(10, 10) = 1.
  end type Matrix

  interface operator(*)
    module procedure multMatrix
  end interface

contains

  function multMatrix( this, other ) result( res )
    type (Matrix), intent(in) :: this, other
    type (Matrix) :: res
    
    res % data = MATMUL( this % data, other % data )
  end function multMatrx

end module class_Matrix

program main
  
   use class_Matrix
   type (Matrix) :: a, b, c

   c = a * b
end program
[/bash]

This is a very simple version to get you thinking and started. The Matrix class stores only a 10x10. So you really need to make this allocatable data and provided a constructor to create a matrix type. This can all be placed in a module and you have yourself a nice matrix type. You may also want to look at reference counting your objects, but I won't get into any details here. You can check my project page below and download it to see examples of things like that. Have fun.

[david_car]

Sorry. The code syntaxhighlighter did some mangling to the code above. Not sure why.

[Steven_L_Intel1]

The problem with the syntaxhighlighter has been reported.

[nooj]

David -

> To be able to override operators you need to define a derived type that you can override the * operator on.

Thanks so much! This was the key I was looking for, and explains why I kept seeing examples like the one you gave.

I'll take a look at your modules and see if they work for me. I've found that for new code, I'm interested in readability and ease of use; but after things mature and become bulletproofed, I adjust them for speed. So I would like to be able to use compatible data structures for both.

- Nooj

[joseph-krahn]

My preference would be to use named operators, such as ".matmul." and ".dot.", rather than replace the scalar multipl ication operator.

David's example uses a matrix class object, rather than directly using Fortran matrices. A matrix class can provide more features, but makes your code specific to that specific class implementation.

&nb sp;

If you want to operate on standard Fortran matrices, you will need to provide wrapper procedures for each al lowed combination of arguments. Here is an example for the three forms of REAL(4) calls. The intrinsic MATMUL handles al l combinations of numeric types, as well as logical types. To support all valid combinations of 4 and 8 byte logical, in teger, real, and complex types, you will need 2*2+6*6=40 sets of three procedures shown here. (Logical and numeric types don't mix.) So, it is not hard to support one or two types, but is excessive to make a completely general interface.

The source-code formatter seems to be randomly chopping lines, so I will include this as an attachment as well.

[fortran]module matmul_op
  interface operator(.matmul.)
    module procedure matmul22_r4
    module procedure matmul12_r4
    module procedure matmul21_r4
  end interface operator(.matmul.)
contains
  pure &
  function matmul22_r4(m1,m2) result(prod)
    real(4), intent(in) :: m1(:,:), m2(:,:)
    real(4) :: prod(size(m1,1),size(m2,2))
    prod = matmul(m1,m2)
  end function matmul22_r4
  pure &
  function matmul21_r4(m1,m2) result(prod)
    real(4), intent(in) :: m1(:,:), m2(:)
    real(4) :: prod(size(m1,1))
    prod = matmul(m1,m2)
  end function matmul21_r4
  pure &
  function matmul12_r4(m1,m2) result(prod)
    real(4), intent(in) :: m1(:), m2(:,:)
    real(4) :: prod(size(m2,2))
    prod = matmul(m1,m2)
  end function matmul12_r4
end module matmul_op
program matmul_test
  use matmul_op
  real(4) :: v1(3) = (/ 2,3,4 /)
  real(4) :: m1(2,3) = reshape((/ 2,3,4, 3,4,5 /),(/2,3/))
  real(4) :: m2(3,2) = reshape((/ 2,3, 3,4, 4,5 /),(/3,2/))
  write(*,*) m1.matmul.m2
  write(*,*) matmul(m1,m2)
  write(*,*)
  write(*,*) v1.matmul.m2
  write(*,*) matmul(v1,m2)
  write(*,*)
  write(*,*) m1.matmul.v1
  write(*,*) matmul(m1,v1)
  write(*,*)
  write(*,*) m1*transpose(m2) ! just to show that * is still usable
end program matmul_test
[/fortran]

[david_car]

Hey Joe. I'm glad you brought this up. The PyF95++ link below is to my Fortran templating open source project. If you use it, you only have to maintain just one code base. It mimics the templating you'll find in C++ and also comes with the a template library with lists, queues, stacks, hashtable/dictionary, etc... I welcome input and contributions. With it, a templated version of the original Matrix class I posted would look like:

[bash]template 
module class_Matrix 
 
  AutoUse
  end AutoUse

  implicit none

  AutoPrivate
  end AutoPrivate

  type Matrix
    __OBJECT__ :: data(__ROW__, __COL__) = 1.
  end type Matrix

 AutoInterface
 end AutoInterface

contains
 
  function multMatrix( this, other ) result( res )  as operator(*)
     type (Matrix<__OBJECT__, __ROW__, __COL__>), intent(in) :: this, other
     type (Matrix<__OBJECT__,__ROW__,__COL__>) :: res

     res = MATMUL( this, other )
  end function multMatri

end module class_Matrix
end template 

program main
  AutoUse
  end AutoUse

  type (Matrix) :: a, b, c
  type (Matrix) :: ia, ib, ic

  ic = a * b
  ic = ia * ib

end program [/bash]

Now you can just stamp out any type of matrix you would like with only one code base. There is an upcoming ACM Fortran Forum article with more details. The "Auto*" blocks take care of generating the interfaces and importing the correct modules. All the best.

[david_car]

Obviously the syntaxhighlighter is not properly display the less than < and greater than > symbols. Sorry about that. Just click "view plain" and you should see it properly.

[Steven_L_Intel1]

I've reported the angle bracket problem. This is apparently a side-effect of a recent fix for the spurious blanks,