Hello. I have been using BLAS/LAPACK with Fortran for some time now, and recently I have used it with C which lead me to a confusion on how it is called with Fortran.
In Fortran, if you want to call a subroutine such as dgemm for matrix-matrix multiply and you want to use a submatrix of say a matrix A in the multiplication (for example the submatrix A(2:, 2:)) as opposed to the entire matrix, the usual process I would follow would be to just specify the top-left entry A(2, 2) of the submatrix and the proper submatrix size and its leading dimension. An example is given in the Fortran code below, where the multiplication C = A(2:, 2:) * B is performed.
I just took this as a way of life in Fortran, but recently when I used BLAS/LAPACK in C it seemed much more clear why you specify the first entry of the submatrix in the subroutine call. Below the Fortran code is an equivalent C code (I am storing C matrices as 1d arrays in column-major order). When calling dgemm in C, you specify a pointer to the top-left entry of the submatrix (in the example code this is &A[1 + lda*1]) and this makes sense to me since an array in C is secretly a pointer to the first element of the array (so in calling a submatrix you specify a pointer to the first element of the array and the code knows how far in memory to jump to adjacent entries).
Thinking back to Fortan, I don’t really know why you must also specify the first entry of the matrix. Are arrays also stored basically as pointers in Fortran? I tried looking through Modern Fortran Explained, but I haven’t been able to find an answer (it could be there but I just haven’t been able to find it). How would I write my own code that looks like BLAS/LAPACK calls, since to me it looks like an argument such as A(2, 2) is just a double, so I would get compiler errors trying to index a double. If there is a reference for how arrays in Fortran are stored/used, that would be appreciated also.
Thanks for the time.
program main
implicit none
double precision :: A(3, 3), B(2, 2), C(2, 2)
A = reshape([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0], shape=[3, 3])
B = reshape([1.5, 2.5, 3.5, 4.5], shape=[2, 2])
! A = 1.0 4.0 7.0 A(2:, 2:) = 5.0 8.0
! 2.0 5.0 8.0 6.0 9.0
! 3.0 6.0 9.0
!
! B = 1.5 3.5
! 2.5 4.5
! C = 1.0 * A(2:, 2:) * B
call dgemm('n', 'n', 2, 2, 2, 1.0d0, A(2, 2), 3, B, 2, 0.0d0, C, 2)
! C = 27.5 53.5
! 31.5 61.5
endprogram main
#include <stdlib.h>
void dgemm_(char* transa, char* transb, int* m, int* n, int* k,
double* alpha, double* A, int* lda, double* B, int* ldb,
double* beta, double* C, int* ldc);
int main()
{
int m, n, k;
int lda, ldb, ldc;
m = 2; n = 2; k = 2;
lda = 3; ldb = 2; ldc = 2;
char transa, transb;
transa = 'n'; transb = 'n';
double alpha, beta;
alpha = 1.0; beta = 0.0;
double *A, *B, *C;
A = (double*) malloc(lda * lda * sizeof(double));
B = (double*) malloc(ldb * ldb * sizeof(double));
C = (double*) malloc(ldc * ldc * sizeof(double));
A[0 + lda*0] = 1.0; A[0 + lda*1] = 4.0; A[0 + lda*2] = 7.0;
A[1 + lda*0] = 2.0; A[1 + lda*1] = 5.0; A[1 + lda*2] = 8.0;
A[2 + lda*0] = 3.0; A[2 + lda*1] = 6.0; A[2 + lda*2] = 9.0;
B[0 + ldb*0] = 1.5; B[0 + ldb*1] = 3.5;
B[1 + ldb*0] = 2.5; B[1 + ldb*1] = 4.5;
dgemm_(&transa, &transb, &m, &n, &k,
&alpha, &A[1 + lda*1], &lda, B, &ldb,
&beta, C, &ldc);
free(A); free(B); free(C);
return 0;
}