Thanks. Now I can define some C functions with matrix arguments and use them from Fortran.
The C code
sum_mat.c
#include "sum_mat.h"
#include <stdio.h>
void print_mat(const int nr, const int nc, const float x[nr][nc])
{
for (int ir=0; ir<nr; ir++) {
for (int ic=0; ic<nc; ic++) {
printf(" %f",x[ir][ic]);
}
printf("\n");
}
}
void print_mat_transpose(const int nr, const int nc, const float x[nr][nc])
{
for (int ic=0; ic<nc; ic++) {
for (int ir=0; ir<nr; ir++) {
printf(" %f",x[ir][ic]);
}
printf("\n");
}
}
float sum_mat(const int nr, const int nc, const float x[nr][nc])
{
float xsum = 0.0;
for (int ir=0; ir<nr; ir++)
for (int ic=0; ic<nc; ic++)
xsum += x[ir][ic];
return xsum;
}
void sum_rows(const int nr, const int nc, const float x[nr][nc], float sums[nr])
{
for (int ir=0; ir<nr; ir++) {
sums[ir] = 0.0;
for (int ic=0; ic<nc; ic++) {
sums[ir] += x[ir][ic];
}
}
}
void scale_mat(float xscale, const int nr, const int nc, float x[nr][nc])
{
for (int ir=0; ir<nr; ir++) {
for (int ic=0; ic<nc; ic++) {
x[ir][ic] *= xscale;
}
}
}
with header file
void print_mat(const int nr, const int nc, const float x[nr][nc]);
void print_mat_transpose(const int nr, const int nc, const float x[nr][nc]);
float sum_mat(const int nr, const int nc, const float x[nr][nc]);
void sum_rows(const int nr, const int nc, const float x[nr][nc], float sums[nr]);
void scale_mat(float xscale, const int nr, const int nc, float x[nr][nc]);
can be called from a Fortran code xxmat.f90
program xxmat
use iso_c_binding, only: c_int, c_float
implicit none
integer, parameter :: nr = 2, nc = 3
real(kind=c_float) :: x(nr,nc),row_sums(nr),col_sums(nc)
integer :: ir,ic
interface
subroutine print_mat(nrow_c, ncol_c, x) bind(c,name="print_mat_transpose")
import c_int, c_float
integer(kind=c_int), intent(in), value :: nrow_c, ncol_c
real(kind=c_float) , intent(in) :: x(ncol_c,nrow_c)
end subroutine print_mat
!
function sum_mat(nrow_c,ncol_c,x) result(xsum) bind(c)
import c_int, c_float
integer(kind=c_int), intent(in), value :: nrow_c, ncol_c
real(kind=c_float) , intent(in) :: x(ncol_c,nrow_c)
end function sum_mat
!
subroutine sum_cols(nrow_c,ncol_c,x,col_sums) bind(c,name="sum_rows")
import c_int, c_float
integer(kind=c_int), intent(in), value :: nrow_c, ncol_c
real(kind=c_float) , intent(in) :: x(ncol_c,nrow_c)
real(kind=c_float) , intent(out) :: col_sums(nrow_c)
end subroutine sum_cols
!
subroutine scale_mat(xscale, nrow_c, ncol_c, x) bind(c,name="scale_mat")
import c_int, c_float
real(kind=c_float) , intent(in), value :: xscale
integer(kind=c_int), intent(in), value :: nrow_c, ncol_c
real(kind=c_float) , intent(in) :: x(ncol_c,nrow_c)
end subroutine scale_mat
!
end interface
print*,"matrix"
do ir=1,nr
do ic=1,nc
x(ir,ic) = real(10*ir + ic, kind = c_float)
end do
print "(*(1x,f5.1))",x(ir,:)
end do
print "(/,a)", "calling print_mat"
call print_mat(nrow_c=nc, ncol_c=nr, x=x)
print*
print*,"sum_mat(nr,nc,x),sum(x)=",sum_mat(nr,nc,x),sum(x)
call sum_cols(nrow_c=nc, ncol_c=nr, x=x, col_sums=col_sums)
print*,"col_sums=",col_sums
print*,"sum(x,dim=1) =",sum(x,dim=1)
call scale_mat(xscale=10.0, nrow_c=nc, ncol_c=nr, x=x)
print "(/,a)","matrix after scaling:"
call print_mat(nrow_c=nc, ncol_c=nr, x=x)
end program xxmat
Compiling with
gcc -c sum_mat.c
gfortran sum_mat.o xxmat.f90
gives output
matrix
11.0 12.0 13.0
21.0 22.0 23.0
calling print_mat
11.000000 12.000000 13.000000
21.000000 22.000000 23.000000
sum_mat(nr,nc,x),sum(x)= 102.000000 102.000000
col_sums= 32.0000000 34.0000000 36.0000000
sum(x,dim=1) = 32.0000000 34.0000000 36.0000000
matrix after scaling:
110.000000 120.000000 130.000000
210.000000 220.000000 230.000000
The only tricky part was that real :: x(n1,n2) corresponds to float x[n2][n1], which I often stumbled over. A C function to compute row sums corresponds to a Fortran function to compute column sums, as the interface for sum_cols illustrates, and the argument names in that interface remind me what is being done. By defining a wrapper with assumed shape argument in a module, such as
module mat_mod
use iso_c_binding, only: c_int, c_float
implicit none
interface
subroutine sum_cols_base(nrow_c,ncol_c,x,col_sums) bind(c,name="sum_rows")
import c_int, c_float
integer(kind=c_int), intent(in), value :: nrow_c, ncol_c
real(kind=c_float) , intent(in) :: x(ncol_c,nrow_c)
real(kind=c_float) , intent(out) :: col_sums(nrow_c)
end subroutine sum_cols_base
end interface
contains
!
function sum_cols(x) result(col_sums)
real(kind=c_float), intent(in) :: x(:,:)
real(kind=c_float) :: col_sums(size(x,2))
call sum_cols_base(size(x,2),size(x,1),x,col_sums)
end function sum_cols
end module mat_mod
you can define a procedure that can be accessed without having to remember how Fortran and C matrix ordering differs, for example
program xxmat
use iso_c_binding, only: c_float
use mat_mod, only: sum_cols
implicit none
integer, parameter :: nr = 2, nc = 3
real(kind=c_float) :: x(nr,nc)
integer :: ir,ic
print*,"matrix"
do ir=1,nr
do ic=1,nc
x(ir,ic) = real(10*ir + ic, kind = c_float)
end do
print "(*(1x,f5.1))",x(ir,:)
end do
print "(a,*(1x,f5.1))","column sums =",sum_cols(x)
end program xxmat
! output:
! matrix
! 11.0 12.0 13.0
! 21.0 22.0 23.0
! column sums = 32.0 34.0 36.0
The C code shown cannot be compiled with g++. I see how to pass matrices from Fortran to C, but what about C++?