Hello everyone,

I have a physics background and have been mainly coding mathematica, python/cython , and julia so far (Though not efficiently, just to get things done). Lately, I realized that better code makes your life so much easier so I spent a fair amount of time improving some code I had. Due to the many python/cython/numba/scipy incompatible bits that I have now I decided to redoit, and take that time to learn more about proper SWE, since it’s not that much code anyway, I decided I wanted to try low level languages so I started doing the same code in fortran and cpp and though I’m doing it blindly in the sense I don’t relly know the programming languages properly I’m just trying to implement with lots of help from old stack overflow answers. Here’s what I got so far

module Qobjmod
implicit none
private
public :: Qobj, operator (==), operator (*),dag

type :: Qobj
complex(8), allocatable :: data(:,:)
contains
procedure :: init
procedure :: print_matrix
end type Qobj

interface operator (==)
module procedure is_equal
end interface

interface operator (*)
module procedure multMatrix
module procedure multScalar1
module procedure multScalar2
end interface

contains

! Initialize Qobj with a matrix
subroutine init(this, mat)
class(Qobj), intent(inout) :: this
complex(8), intent(in) :: mat(:,:)
allocate(this%data(size(mat, 1), size(mat, 2)))
this%data = mat
end subroutine init

! Overload == operator for Qobj
logical function is_equal(this, other)
class(Qobj), intent(in) :: this, other
real(8), parameter :: tol = 1.0e-10
is_equal = all(abs(this%data - other%data) < tol)
end function is_equal

! Overload * operator for Qobj
function multMatrix(this, other) result(res)
type(Qobj), intent(in) :: this, other
type(Qobj) :: res
! Allocate result
allocate(res%data(size(this%data, 1), size(other%data, 2)))
! Perform matrix multiplication
res%data = matmul(this%data, other%data)
end function multMatrix

! Overload * operator for scalar multiplication (scalar * Qobj)
function multScalar1(scalar, mat) result(res)
complex(8), intent(in) :: scalar
type(Qobj), intent(in) :: mat
type(Qobj) :: res

! Allocate result
allocate(res%data(size(mat%data, 1), size(mat%data, 2)))
! Perform scalar multiplication
res%data = scalar * mat%data
end function multScalar1

! Overload * operator for scalar multiplication (Qobj * scalar)
function multScalar2(mat, scalar) result(res)
type(Qobj), intent(in) :: mat
complex(8), intent(in) :: scalar
type(Qobj) :: res

! Allocate result
allocate(res%data(size(mat%data, 1), size(mat%data, 2)))
! Perform scalar multiplication
res%data = mat%data * scalar
end function multScalar2

! dagger operation
function dag(this) result(res)
type(Qobj), intent(in) :: this
type(Qobj) :: res
! Allocate result
allocate(res%data(size(this%data, 1), size(this%data, 2)))
! Perform matrix multiplication
res%data = transpose(conjg(this%data))
end function dag

! Procedure to print the matrix
subroutine print_matrix(this)
class(Qobj), intent(in) :: this
print *, "Matrix:"
print *, this%data
end subroutine print_matrix

end module Qobjmod


My question here is whether I’m doing something that is not a best practice in fortran? or something that’s just wrong/discouraged?

Aditionally I have some questions about testing and formatting. I’ve been using fpm for tests and so far so good, but couldn’t find any documentation on how to get code coverage from those test?

Is there any formatter for fortran code? similar to black or autopep for python?

Some other place to interact with fortran users, like a discord channel?

and lastly some general advice for you high performant coders about the problem I will try to tackle at the end. I basically will have some large number of matrices of N dimension (the number is about N^{2} in the worst case scenario) and for all pair combinations I need to compute something similar to

\sum_{i,j} \alpha(i,j,t) (A_{i}A_{j}^{T} \otimes 1)-( A_{i} \otimes A_{j}^{T} )

Where the A_{i}s are the matrices I mentioned and the alphas are some time dependent coefficients. For now I precompute the \alpha get an array (one element for each t) and multiply it by the matrix for each term of the sum. I need to use generators in python due to the memory cost of storing all the individual sum terms. So I would like general advice on how to implement it as memory efficient as possible (of course also as fast but memory is a bigger priority for now). I’m sorry if this post is all over the place, lots of questions and don’t really know anyone I could ask

EDIT: I just realized getting eigenvectors and eigenvalues is not so straight-forward. I did see however that functions were recently added to the stdlib. are they as performant as calling lapack? I also couldn’t find documentation for them

I’ll give you some micro feedback on your first question. This (8) is not best practice because the kind numbers are processor dependent. There’s some debate around here about whether you should use the iso_fortran_env module constants real64, real128, etc which define a bit storage size, or use the selected_real_kind() intrinsic to specify a level of decimal precision and let the compiler choose the storage size.

On assignment to allocatable arrays, you generally don’t need to pre-allocate, it will do so automatically. The only downside I can think of to explicitly allocating is that it always performs the allocation step, while auto-allocation is supposedly only done if there’s a dimension mismatch with the already-allocated left hand side, potentially saving time.

I’m on my phone and not really able to analyze your design. But what do you gain by using a derived type with a single field, versus using intrinsic complex arrays? Is this just a demo code for something more complicated?

1 Like

In is_equal and mult_matrix you probably should add a dimension compatibility check before performing operations that could suffer a runtime error.

1 Like

The allocated dimensions of the transpose result are incorrect, they need to be swapped.

1 Like

Thanks a lot for your comments, wasn’t aware that was processor dependent. I think I’ll go with the iso_fortran_env.

I also wasn’t aware that preallocation wasn’t necessary, I kind of took that from a forum and then assumed it needed to be used everywhere. will something like deallocation also happen if I assign another value to the same variable? I had some funny issues with python jax reserving too much memory for 16x16 matrices when I started using it

It’s a demo for something more complicated.I’m just taking it one step at a time and thought that I had enough code to seek for advice, since I’ve only been doing fortran for a couple days.

On the checks for dimensionality that is true, didn’t think of it.

Didn’t notice the issue with the transpose, well spotted, I only tested with square matrices, but also need matrix vector multiplication, so definitely need to fix it.

Yes, my understanding is that deallocation and reallocation happen automatically if the dimensions don’t match. Someone else might correct me if that’s wrong.

1 Like

You might have already learned this the hard way, but this will not print something that looks like a matrix. It will print a list of numbers left-to-right with arbitrary line breaks at some regular interval, which are the elements of the matrix in column order. In other words, the matrix that you might write in Matlab as:

[ 1+0i  2+0i  3+0i  4+0i;
5+0i  6+0i  7+0i  8+0i;
9+0i 10+0i 11+0i 12+0i]


with list-directed output print *,data might print in a processor-dependent way as

 (1.000000,0.0000000E+00) (5.000000,0.0000000E+00) (9.000000,0.0000000E+00)
(2.000000,0.0000000E+00) (6.000000,0.0000000E+00) (10.00000,0.0000000E+00)
(3.000000,0.0000000E+00) (7.000000,0.0000000E+00) (11.00000,0.0000000E+00)
(4.000000,0.0000000E+00) (8.000000,0.0000000E+00) (12.00000,0.0000000E+00)


While this looks like a clean transpose, the number of values per line is processor-dependent.

You’ll want to develop a format string to display your matrix in a pretty way. It will likely involve implied do-loops. For example:

Program with complex matrix format examples
implicit none

complex, dimension(3,4) :: a = reshape(cmplx([1,2,3,4,5,6,7,8,9,10,11,12]),[3,4],order=[2,1])
integer :: i

a = a + (0.0,1.0)

print *, a
print *
print '(4("(",g0,",",g0,")",1x))', (a(i,:), i=1,3)
print *
print '(("[",1x,4(S,G0,SP,G0,"i",1x),"]"))', (a(i,:), i=1,3)
print *

a = a - (0.0,2.0)

print *, a
print *
print '(4("(",g0,",",g0,")",1x))', (a(i,:), i=1,3)
print *
print '(("[",1x,4(S,G0,SP,G0,"i",1x),"]"))', (a(i,:), i=1,3)
print *

end

Output of program using ifx 2024.1 on WSL
 (1.000000,1.000000) (5.000000,1.000000) (9.000000,1.000000) (2.000000,1.000000)
(6.000000,1.000000) (10.00000,1.000000) (3.000000,1.000000) (7.000000,1.000000)
(11.00000,1.000000) (4.000000,1.000000) (8.000000,1.000000) (12.00000,1.000000)

(1.000000,1.000000) (2.000000,1.000000) (3.000000,1.000000) (4.000000,1.000000)
(5.000000,1.000000) (6.000000,1.000000) (7.000000,1.000000) (8.000000,1.000000)
(9.000000,1.000000) (10.00000,1.000000) (11.00000,1.000000) (12.00000,1.000000)

[ 1.000000+1.000000i 2.000000+1.000000i 3.000000+1.000000i 4.000000+1.000000i ]
[ 5.000000+1.000000i 6.000000+1.000000i 7.000000+1.000000i 8.000000+1.000000i ]
[ 9.000000+1.000000i 10.00000+1.000000i 11.00000+1.000000i 12.00000+1.000000i ]

(1.000000,-1.000000) (5.000000,-1.000000) (9.000000,-1.000000)
(2.000000,-1.000000) (6.000000,-1.000000) (10.00000,-1.000000)
(3.000000,-1.000000) (7.000000,-1.000000) (11.00000,-1.000000)
(4.000000,-1.000000) (8.000000,-1.000000) (12.00000,-1.000000)

(1.000000,-1.000000) (2.000000,-1.000000) (3.000000,-1.000000) (4.000000,-1.000000)
(5.000000,-1.000000) (6.000000,-1.000000) (7.000000,-1.000000) (8.000000,-1.000000)
(9.000000,-1.000000) (10.00000,-1.000000) (11.00000,-1.000000) (12.00000,-1.000000)

[ 1.000000-1.000000i 2.000000-1.000000i 3.000000-1.000000i 4.000000-1.000000i ]
[ 5.000000-1.000000i 6.000000-1.000000i 7.000000-1.000000i 8.000000-1.000000i ]
[ 9.000000-1.000000i 10.00000-1.000000i 11.00000-1.000000i 12.00000-1.000000i ]


To handle arbitrary-size matrices you’ll probably need to build your format strings programmatically, generating repeat counts and implied do-loop limits based on the matrix dimensions.

That’s correct. And it also means you can actually omit the allocate statements for the code you posted.

A few things I noticed:

• If the data component of the Qobj derived type is not intended to be private, then the init procedure is unnecessary. Otherwise, you could make it a custom constructor with interface Qobj; module procedure ...; end interface.
• Your code is not checking whether the data component is allocated or not, before performing operations on it.
• If the compiler(s) you’re using support the generic statement for the type-bound-procedure-part of derived types, then you can do the operator overloading there instead.
• the print_matrix subroutine will print one column after another (not one row after another).
• With the same generic statement, you could overload write(formatted) as well, which will also cover the print statement.
1 Like

When we started writing our electronic structure code (about 20 years ago), we also started by defining objects (like wavefunctions) and using operator overloading to manipulate them.

After about three weeks we went back to flat arrays and procedural code, and we’ve never looked back. There is only one derived type in the entire code, used for a large ragged array. The reason is somewhat subtle, but we were losing the flexibility offered by procedural code and having to define increasingly complicated object operations. Linus Torvalds puts it better but less politely for C++.

If your problem is similar, i.e. flat arrays, basic real/complex/integer types, straight-forward mathematical algorithms and long-term development, then you may be better off by avoiding OOP code altogether.

You may also like to use BLAS or LAPACK code where possible. There are highly optimised libraries BLAS/LAPACK libraries such as Intel’s MKL, flame/blis and OpenBLAS which can dramatically increase the speed of your code.

3 Likes

Welcome to the Fortran Discourse, @mcditoos

Concerning Discord channels, you can try:

2 Likes

From C++ core guidelines,

Overload only for operations that are roughly equivalent.
Core Guideline C.163

It is a bad idea to equate an “almost equal” with “equal”. Other programmers will expect you to adhere to the expected meaning of operators like == (strict equality). Moreover, a fixed absolute tolerance won’t work in all scenarios. Instead, I’d recommend defining a function like numpy.isclose. This could be an elemental function, which would automatically cover also comparisons between matrices and scalars. If you really want to, you could define a custom operator,

type(Qobj) :: A, B
logical :: predicate

predicate = is_close(A, B)
predicate = (A .is_close_to. B)


even if this is just a form of syntactic sugar.

I’ll also note that the previous C++ core guideline is,

Overload operations that are roughly equivalent
Core guideline C.162

There is a tension between these two guidelines, and it’s up to you to decide if introducing a function into an overload set is the right choice. Keep in mind that overloaded functions may carry expectations on their behaviour. If these expectations are broken, things will not work well.

2 Likes

An fpm package for displaying small matrices is at

GitHub - urbanjost/M_display: An fpm(1) package for displaying small matrices based on dispmodule(3f)

which wraps dispmodule in fpm.

2 Likes

I don’t know if it’s only me but whenever I click on those links I get the message “invitation link has expired”

Thanks, Indeed the function didn’t print anything nice, just had it there to check by eye that things were computed correctly . Thanks for all your useful comments !

There’s a point to it, I could do it all with procedural code. At some point I want certain matrices to behave differently than this ones and that’s why I do OOP since it’s the way I did it in other programming languages.

Directly calling BLAS or LAPACK was one of the main reasons to learn a compiled language but it seems a lot more daunting than I thought. I saw the stdlib implementation makes it easier but haven’t seen any benchmarks and haven’t been able to set it up with fpm correctly . In hindsight I probably should have waited until I implemented kronecker products/ eigenvectors /eigenvalues beforer making this post.

Thanks. I’m afraid I don’t get the type-bound-procedure part. Any reference or tutorial on it?

That is true I tend to think of == as np.isclose for matrices so something like that was the intention but since I’m only learning fortran for now decided to use the fixed tolerance, due to some code implementations I’ve been working with and === as strict equality.

I wasn’t familiar with those core guidelines, I guess my main reason to overload most operations is to have something similar to what I have been working with in python so that in case the Fortran implementation does good. It would work exactly the same way (mainly so that my collaborators would use it, there’s always certain inertia to use new notation) though there’s probably better ways to do it, will need to read up .

Thanks a lot for the guidelines, wasn’t familiar with them but it’s nice to see a pep equivalent. Should get familiar with it

Operator overloading was a Fortran 90 thing, and the OOP facilities were added in Fortran 2003, so you have two ways to do the same thing, e.g.:

....
type, public :: qobj
complex(real64), allocatable :: data(:,:)
contains
procedure, private ::  mult_matrix, mult_scalar1, mult_scalar2
generic :: operator(*) => mult_matrix, mult_scalar1, mult_scalar2
procedure, private :: write_qobj
generic :: write(formatted) => write_qobj
end type
...
contains
....
subroutine write_qobj(this, unit, iotype, v_list, iostat, iomsg)
class(qobj), intent(in) :: this
integer, intent(in) :: unit
character(*), intent(in) :: iotype
integer, intent(in) :: v_list(:)
integer, intent(out) :: iostat
character(*), intent(inout) :: iomsg
....
end subroutine


In Fortran, type-bound procedure is the fancy way of saying method (i.e., passed object or receiver marked with class, dynamic dispatch, etc.).

This post lists some interesting stuff related to OOP in Fortran.

1 Like

Indeed, those posts are two years old. Either the invitation link is not valid anymore or the discord channel has disappeared.

Another place to discuss Fortran is:

https://fortran-lang.discourse.group/t/comp-lang-fortran-37-years-of-discourse/

1 Like

@jkd2022 100% agree with your whole comment. I have exactly the same experience.

I would add, that it is possible to program in OO way with Go-like interfaces and compose your project that way and have things reasonably decoupled. Rust also allows this kind of workflow.

LFortran is written in C++ and we use a few classes, mainly visitors to visit the intermediate representation and change it. We had to write our own vector class and custom allocators and such, for performance. In fact, for performance, the visitor class is just a zero-cost abstraction over a C style struct and casting, which I found the fastest way (after a lot of benchmarking) to represent and visit a tree in memory. So we could have written LFortran in C only. However, it would be more verbose, and that is probably the biggest advantage of C++ that it can keep the verbosity down by having some carefully designed abstractions and simpler syntax.

For array-based numerical code, the main abstraction is a multi-dimensional array, and Fortran has it, so the need for OO is greatly reduced. For my Fortran codes, I typically also have no classes and usually just one (or very few) derived types for the whole project, if needed. (But I also know other Fortran programmers that like to have a lot of classes.)

1 Like