Leading dimension

如果你用LAPACK解过矩阵本征值问题，你一定会接触到这样一个名词，“leading dimension”，比如在函数zheev中。我想绝大部分人在第一次接触这个词的时候都不明白它到底是什么意思。以前我也不明白，今天索性搜了一把，在下面找到了答案。
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?p=661&sid=67c66465dedfcbb6e0612cca7647698f

Suppose that you have a matrix A of size 100x100 which is stored in an array 100x100. In this case LDA is the same as N. Now suppose that you want to work only on the submatrix A(91:100 , 1:100); in this case the number of rows is 10 but LDA=100. Assuming the fortran column-major ordering (which is the case in LAPACK), the LDA is used to define the distance in memory between elements of two consecutive columns which have the same row index. If you call B = A(91:100 , 1:100) then B(1,1) and B(1,2) are 100 memory locations far from each other. 

其实之所以设LDA（leading dimension）这个参数主要是考虑到fortran是“列优先”存储数组的原因。这里要解本征值的矩阵是NxN大小的，但是存储这个矩阵的数组A却并不一定非得是NxN大小，可以是M1xM2大小，其中 M1≥N，M2≥N，NxN矩阵要存放在M1xM2数组的左上角，即A(1:N,1:N)部分。这样，当把数组A传给zheev时，zheev通过N来知道要解的矩阵是多大的，通过LDA来知道同一行中相邻两列的元素在内存中相距多远，可见，LDA=M2，其实就是数组A的一列的元素个数，也就是“实际存储时的第一维”的大小。正是由于fortran列优先存储数组，才使得概念上的第一维（行）与实际存储时的第一维（列）不一样。

注1：如果参数A的位置就用数组片段来调用的话则令当别论。例如同样A的大小为M1xM2，那么若调用zheev时，A参数位置处的实参是A(1:N,1:N)的话，则LDA位置处的实参应该是N，而不是M2！

注2：zheev的输出本征矢时，A的每一列代表一个本征矢。

附：调用zheev解厄米矩阵本征值
program main
implicit none
  CHARACTER  JOBZ, UPLO
  INTEGER::  INFO=0
  INTEGER,parameter::N=3, LDA=N, LWORK=2*N-1
  REAL*8     RWORK(3*N-2), W(N)
  COMPLEX*16 A(LDA,N), WORK(LWORK),B(5,3)
  data A/1,2,3,2,4,5,3,5,6/
  B(2:4,:)=A

 !输出中A的每一列是一个本征矢。
 !call ZHEEV('V', 'U', N, A, LDA, W, WORK, LWORK, RWORK, INFO )
  call myZHEEV(N,W,A)
  print*, W
  print*, '--------------------   --------------------    ------------------'
  print*, real(A(1,:))
  print*, real(A(2,:))
  print*, real(A(3,:))
  print*
  call ZHEEV('V', 'U', N, B(2:5,:), 4, W, WORK, LWORK, RWORK, INFO )
  print*, W
  print*, '--------------------   --------------------    ------------------'
  print '(3F25.15)', real(B(1,:))
  print '(3F25.15)', real(B(2,:))
  print '(3F25.15)', real(B(3,:))
  print '(3F25.15)', real(B(4,:))
  print '(3F25.15)', real(B(5,:))

  contains
  subroutine myzheev(N,W,A)
    implicit none
    integer::N,INFO=0
    real*8::W(N),RWORK(3*N-2)
    complex*16::A(N,N),WORK(2*N-1)
    call ZHEEV('V', 'U', N, A, N, W, WORK, 2*N-1, RWORK, INFO )
  endsubroutine

end

输出如下：

 -0.515729471589256       0.170915188827180        11.3448142827621
 --------------------   --------------------    ------------------
  0.736976229099579       0.591009048506103       0.327985277605682
  0.327985277605681      -0.736976229099578       0.591009048506103
 -0.591009048506103       0.327985277605682       0.736976229099578

 -0.515729471589256       0.170915188827180        11.3448142827621
 --------------------   --------------------    ------------------
        0.000000000000000        0.000000000000000        0.000000000000000
        0.736976229099579        0.591009048506103        0.327985277605682
        0.327985277605681       -0.736976229099578        0.591009048506103
       -0.591009048506103        0.327985277605682        0.736976229099578
        0.000000000000000        0.000000000000000        0.000000000000000