Any intrinsic functions to compute leading dimensions?

[styc]

Hi,

I need to pass some assumed-shape 2D arrays to BLAS and LAPACK. Those subroutines require leading dimensions of arrays be passed along. To avoid triggering the compiler to emit code that generates contiguous copies of the arrays, I can use something like A(1, 1). But then I am stuck. I need the real leading dimension. SIZE(A, 1) does not help because A can be an array slice. After some twist, I figure out that I can compute it via (LOC(A(2, 1)) - LOC(A(1, 1)) / 8 assuming double precision. Are there any neater ways to do this?

Thanks,

[Steven_L_Intel1]

I'm a bit uncertain of what you are looking for - do LBOUND and UBOUND help? Note that if the argument is an array slice, LBOUND will always be 1.

[styc]

Quoting - Steve Lionel (Intel)

I'm a bit uncertain of what you are looking for - do LBOUND and UBOUND help? Note that if the argument is an array slice, LBOUND will always be 1.

LBOUND and UBOUND do not quite help either. Consider this example. Suppose I have an A(-1000:1000, -1000:1000). I pass A(-1000:0, 0:1000) to a subroutine inside which the slice is called B. So B has extents (1:1001, 1:1001). Then I want to feed B to BLAS and LAPACK. So I write CALL SOME_SUBROUTINE(..., B(1, 1), LDB, ...).

The problem here ishow to find outLDB = 2001. According to the 11.1.056 compiler, SIZE(B, 1) = UBOUND(B, 1) - LBOUND(B, 1) + 1 = 1001. So this is not what I want. The LOC trick works, but according to what I have read,the LOC intrinsic functionis not in the language standard. I am wondering if there is any standard way of computing LDB.

[Steven_L_Intel1]

SUM(SIZE(B)) ? SIZE(B,1) just gets you the extent of the first dimension. Isn't the correct value for LDB 2002 and not 2001?

[styc]

Quoting - Steve Lionel (Intel)

SUM(SIZE(B)) ? SIZE(B,1) just gets you the extent of the first dimension. Isn't the correct value for LDB 2002 and not 2001?

It does not make sense to write SUM(SIZE(B)).That does not even compile.

LDBshould beequal to SIZE(A, 1) (= 1000 - (-1000) + 1 = 2001)so that B can be declared as B(LDB, *). My question put in another way, how do I figure out SIZE(A, 1)from inside a subroutine that accepts B(:, :), knowing that B is a slice of A?

[Steven_L_Intel1]

Oh, I get it now.... Sort of. Please explain why 1001 is not the leading dimension. The array slice you have is dimension (1001,1001)

[styc]

Quoting - Steve Lionel (Intel)

Oh, I get it now.... Sort of. Please explain why 1001 is not the leading dimension. The array slice you have is dimension (1001,1001)

Maybe I should really not have used the phrase "leading dimension".

When I pass B(:, :) to BLAS and LAPACK, I want to simply tell them that B is a noncontiguous slice of a bigger matrix A of 2001 rows instead of making a contiguous copy of B. I need a portable way to compute the number 2001.

[Steven_L_Intel1 · Accepted solution]

Ok. There is no standard way of obtaining the stride of a non-contiguous array. I would suggest that you go ahead and use LOC - it may not be standard but it is widely implemented. If you wanted standard syntax, you could use TRANSFER(C_LOC(B(whatever)),1_C_INTPTR_T), though I think I'd rather use LOC here.

[Steven_L_Intel1]

It was pointed out to me that the implementations of BLAS and LAPACK we are familiar with either require a contiguous array and don't have the concept of a stride, or have an explicit Fortran interface with assumed-shape arrays, in which case the strides are handled automatically. So it is not clear to me under what circumstances this LDB would be useful.

[styc]

Quoting - Steve Lionel (Intel)

It was pointed out to me that the implementations of BLAS and LAPACK we are familiar with either require a contiguous array and don't have the concept of a stride, or have an explicit Fortran interface with assumed-shape arrays, in which case the strides are handled automatically. So it is not clear to me under what circumstances this LDB would be useful.

I know that MKL has Fortran 95 interfaces for BLAS and LAPACK. But if I understand the interface code correctly, the accompanying libraries (libmkl_{blas,lapack}95_something.a) actually have array copying code under the hood. There is no point in letting array copying happen as BLAS and LAPACK do indeed support noncontiguous matrices.

[Hirchert__Kurt_W]

Quoting - styc

Quoting - Steve Lionel (Intel)

It was pointed out to me that the implementations of BLAS and LAPACK we are familiar with either require a contiguous array and don't have the concept of a stride, or have an explicit Fortran interface with assumed-shape arrays, in which case the strides are handled automatically. So it is not clear to me under what circumstances this LDB would be useful.

I know that MKL has Fortran 95 interfaces for BLAS and LAPACK. But if I understand the interface code correctly, the accompanying libraries (libmkl_{blas,lapack}95_something.a) actually have array copying code under the hood. There is no point in letting array copying happen as BLAS and LAPACK do indeed support noncontiguous matrices.

I don't know about the MKL versions of LAPACK and the BLAS, but the original FORTRAN 77 versions support only a limited form of noncontiguous matrices. In particular, all strides have to be exact multiples of the element size. This is true of most subsections of arrays, but there are ways in Fortran 90/95/2003 to reference arrays where this property does not hold. Thus, there is no feature for converting an arbitrary assumed-shape array into the form expected by LAPACK.

If you are writing general procedures and have no control over the code that will call them, then it is impossible (in the general case) to avoid the copying.

If you do have control over the callers as well as the callees, then you might be better off emulating LAPACK by using an array form that is guaranteed to be contiguous and separately passing the information that indicates which part of that array is actually to be used.

If you are committed to using the general argument form while depending on the arguments being less general, appropriate use of LOC should be reasonably portable in handling those cases that are possible.

-Kurt