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
. In hindsight I probably should have waited until I implemented kronecker products/ eigenvectors /eigenvalues beforer making this post.
.