I was having a look at the C interoperability after reading some post of @RonShepard about a proposed shallow_transpose.
Now the real point about my post in another thread is not strictly related to have just a way to calculate the transpose but to have a transpose that is definable, so it can be pointed by a pointer or passed as an argument to a procedure with a intent(inout) dummy argument.
Before asking for adding a new feature to Fortran I wanted to see what is the present situation and what we can get with present compilers.
First I was curious if I could use the C interoperability to look at the effect of the transpose function.
That said I wrote a small c function that will show the descriptor of an array:
void display_info(CFI_cdesc_t *src)
{
printf("Address: %p\n", src->base_addr);
for (int i = 0; i < src->rank; i++)
{
printf(" %d", i);
printf(" lower_bound: %ld", src->dim[i].lower_bound);
printf(" extent: %ld", src->dim[i].extent);
printf(" stride: %ld\n", src->dim[i].sm);
}
printf("elem_len %ld\n\n", src->elem_len);
}
And a small Fortran program that will call this C function:
...
interface
subroutine display_info(src) bind(C)
use iso_c_binding
integer(C_INT), intent(in) :: src(..)
end subroutine
end interface
...
integer :: a(3,4)
integer :: i
a = reshape([(i, i=1, size(a))], shape(a))
call display_info(a)
call display_info(transpose(a))
Using fpm to compile with gfortran version 11.3
I get the following output:
Address: 0x7ffd228190c0
0 lower_bound: 0 extent: 3 stride: 4
1 lower_bound: 0 extent: 4 stride: 12
elem_len 4
Address: 0x7ffd228190c0
0 lower_bound: 0 extent: 4 stride: 12
1 lower_bound: 0 extent: 3 stride: 4
elem_len 4
Wow @RonShepard you’ve got your shallow transpose (at least with gfortran).
To check that this is linked to the transpose function I wrote another routine that I called in a non conformant way:
subroutine change(a, val)
integer :: a(:,:), val
a(1,2) = val
end subroutine
An then given an a matrix like:
1 4 7 10
2 5 8 11
3 6 9 12
calling change with: call change(transpose(a), 200)
I got:
1 4 7 10
200 5 8 11
3 6 9 12
Of course it is non conformant and a compiler can do whatever is best to do.
By the way it seems that ifort instead make a copy or, at least, pass a copy of the transpose to the c routine, this desn’t tell anything on the internal behaviour of the ifort compiler, of course.
Then I wrote a shallow_transpose using CFI_Fortran_binding:
where a shallow_transpose procedure is implemented with the following interface:
subroutine shallow_transpose(src, dst, status)
f_type, intent(in), pointer :: src(:,:)
f_type, intent(out), pointer :: dst(:,:)
integer, intent(out) :: status
end subroutine
where f_type is integer(C_INT) and real(C_FLOAT).
The first argument is a pointer with intent in in order to force the calling program to pass at least an array with the target attribute, while the second argument is the output pointer that will contain the transposed matrix.
It should be said that due to the restriction on the modification of the descriptor the routine shallow_transpose it’s still not conformant to the Fortran Standard. I have not found any CFI function able to reorder the index and so I did it on reorder_dim c function that you can see in the repository. The shallow_transpose function is the following one:
void shallow_transpose(CFI_cdesc_t *src, CFI_cdesc_t *dst, int *status) {
int order[2] = {1,0};
*status = CFI_setpointer(dst, src, 0);
if (*status != CFI_SUCCESS) return;
*status = reorder_dim(dst, order, 2);
}
Where first of all I’m using a standard CFI function to create a pointer. But then I reorder some of the dimensions.
At a certain point inside reorder_dim I had to write:
...
// update the descriptor
// This part is not conforming
for(int i=0; i<n; i++) {
v->dim[i].extent = dim[i].extent;
v->dim[i].lower_bound = dim[i].lower_bound;
v->dim[i].sm = dim[i].sm;
}
...
Now, I could have a look at the gfortran source of CFI_setpointer and I’m not expecting that (except for bugs) my reorder dim will break anything.
I checked the program with the compilers gfortran and ifort and both gived correct results, and both output matrix were correctly recognized as not contiguous.
That said I think it should be added to the Standard another CFI function able to reorder the dimensions. That could be handy when coming to interface Fortran with other languages.
My other suggestion would be to have a definable designator to get the transpose of a matrix, mainly when one needs to pass a matrix to an inout dummy argument. Of course if the dummy argument has only the intent(in) attribute the current transpose function will be fine.
That said I have never needed that, so it is just academic.
By the way I have found that the designators:%im and %re are definable for gfortran but not for ifort and I think that ifort is wrong. But tell me if I’m wrong.
More important on the C interoperability will be to have custom deallocator for allocatable arrays, but this is another story.