leading dimension

[serlancelot]

I understood how leading dimension works, it's simply in order to let the computer find the element in the consecutive column with the same row index of the element I'm starting from (at least I hope it is right); but what I don't understood is, why should I put an lda different from the simple number of rows aof the matrix? Which difference can it make, if with a 2 x 6 matrix I put 8 as lda instead of 5 for exemple? Which advantage could give me choosing one instead of the other? How do the people usually choose the lda's?
hope you can help me

[Julien Langou]

Hello,

the leading dimension argument appears for all 2D arrays in the BLAS and
LAPACK interfaces.

If your matrix is simply a 2-by-6 matrix, that starts from A(1,1), then sure
you can work with a 2-by-6 array:

Code: Select all

N=2, M=6, A(1,1), LDA=2.


And LDA = N = number of rows of the matrix.

Now imagine that you have a 20-by-20 array, and you want to work on the
2-by-6 block starting in position A(4,10), in Matlab notation this is
A(4:5,10:15) then you can simply work on this submatrix with

Code: Select all

N=2, M=6, A(4,10), LDA=20.


If you look at any LAPACK subroutines (for example DGETRF), you will
soon understand how useful the leading dimension can be useful.

Another advantage, slightly more complex, is that some leading dimension
are bad. In general, if you take N=LDA=2^9=512 or another power
of 2 or an integer made of high power of 2 (e.g. 2^9+2^8=768), you'll
end up with lots of cache line conflicts and performance can be very poor.
So if you need to work with a matrix of size 512, a good idea is to set
LDA=513 to avoid this problem. This effect and its cure were mentionned
several times in the past. Recently, Sven Hammarling (NAG) and some
colleagues from the University of Manchester came up with nice examples.

Julien

[serlancelot]

Yes, ok, it makes sense, thanks; however I still have some doubts:
how do you communicate to the computers then these codes?

N=2, M=6, A(1,1), LDA=2

N=2, M=6, A(4,10), LDA=20

How do you communicate to the computer that you wanna start from A(4,10) in a 20x20 array?
I guess you do it in the subroutine calling, in the main program, but how?
Could you give me the examples?

What does it mean when they write "On entry, LDA specifies the first dimension of A as declared in the calling (sub) program"? Do they mean the subroutine calling line in the main program, or the first line of the arguments in the code of the subroutine? It's always the same doubt, where do you communicate that you wanna start from A(4,10) in the 20x20 A array, and you wanna work only with a part of it, a 2x6 matrix, exactly starting from A(4,10)?
thanks again and again

[Julien Langou]

Hello,

How do you communicate to the computer that you want to start from A(4,10) in a 20x20 array?

You simply give the pointer to A(4,10). If you are in Fortran, since Fortran passes the
arguments by pointer, you siomply give A(4,10).

For example DGETRF. Here is the interface:

Code: Select all

      SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      DOUBLE PRECISION   A( LDA, * )


In your application code, you write:

Code: Select all

     DOUBLE PRECISION A( 20, 20 )
     INTEGER IPIV(2), INFO


Then you affect values in the array A.
Then you can call DGETRF with:

Code: Select all

      CALL DGETRF( 2, 6, A(4, 10), 20, IPIV, INFO)


This will preform the LU factorization of A( 4:5, 10:15).

Julien.

[serlancelot]

Hey, I tried it but it doesn't work... :cry:
should it work also with Fortran77 ? The compiler gives me "(..)argument arrayness mismatch(..) blablabla (..)dummy arg is a whole array(..)actual arg is array element(..)"...
Maybe it's because Fortran77 doesn't use pointers... I don't know
I tried with sgemm and with dgetrf; in the last case, exactly the way you said me. With sgemm, I built a 4x5 array, I extracted 2 matrixes, 2x5 each, and multiplyied the first with the second transposed, giving a(1,1) as indication for A, and a(3,1) as indication for B in the sub calling.
I give you the line:
call sgemm('n','t',2,2,5,1.0,a(1,1),4,a(3,1),4,1.0,c,2)
But it gave me the same problem...
I can't understand...
And by the way, when you use a subroutine like dgetrf do you copy too, all the codes of the subroutines called by this subroutine, and the subroutines called by the subroutines called by this subroutine, and so on, to add them after the main program? To try "dgetrf" I had to do it, but it's really a hard work, It's never over, I had 13 lines of the main program, to generate the matrix and to call the dgetrf, and more than 2000 lines of subroutine codes pasted together. Is it normal? Do all the people do so?
You know, I'm a very beginner, and I think I ignore something.
Anyway thank you alot, I'm very grateful for the help

[sven]

You might find the example programs at:

http://www.nag.co.uk/lapack-ex/lapack-ex.html

helpful. They show straightforward examples of calling the LAPACK drivers. I don't think that they show examples of passing sections that do not start at (1,1), but at least they should help you get started.

Best wishes,

Sven Hammarling.

[Julien Langou]

Your call to SGEMM makes some sense, SGEMM should perform something like

Code: Select all

   c(1:2,1:2) = 1.0e+00 * c(1:2,1:2) +   1.0e+00 * a(1:2,1:5) * transpose(a(3:4,1:5))


There might have two problems with your code:

1. you need to declare the array a and the array c as follows:
   

   Code: Select all

         real a(4,5), c(2,2)


2. I am not sure exactly how Fortran will interpret the '1.0', you want them to be
   of type REAL (single precision), and a better way to go is to type 1.0E+00. This
   will force the type to be REAL. On my machine/compiler, 1.0 is interpreted as
   the single precision, i.e 1.0E+00. So your call is correct. I am not sure this
   is a standard thing. (Anybody knows?)

What about trying:

Code: Select all

      real a(4,5), c(2,2)
      call sgemm('n','t', 2, 2, 5, 1.0e+00, a(1,1), 4, a(3,1), 4,
     &              1.0e+00, c, 2)


If you are still struggling, copy-paste a small self content code, rather than a line.

Julien.

[Julien Langou]

Hello,
Sven gave me an answer for the point (b) that I have raised. If you type 1.00 in
Fortran this implies the REAL type, so it's the same as 1.00E+00. So your code was
correct on this point (and my compiler as well:). Thanks Sven, for the insight in the
subtility of the Fortran language. And so well, are you sure you did not messed up when
you have declared your arrays?
Julien.

[serlancelot]

eh...I'm still whimpering...
What about, if I sent you my code? so you can take a look, if you want...
maybe I messed in the previous part, I began to write programs in Fortran77 three months ago, maybe I did some mistakes...
you don't care of the italian comments


About the compiler I'm using is the Force 2.0.8
What do you think about? Maybe it is the compiler not very reliable
I cross my fingers