Fortran+Lapackで連立一次方程式Ax=bを解く(dgesv())

ここを参考に。
www.nag-j.co.jp

引数は以下

subroutine dgesv	(	integer 	N,
integer 	NRHS,
double precision, dimension( lda, * ) 	A,
integer 	LDA,
integer, dimension( * ) 	IPIV,
double precision, dimension( ldb, * ) 	B,
integer 	LDB,
integer 	INFO 
)

N→連立方程式の数（＝Aの次数）
NRHS→B行列の列数（ベクトルの場合は１）
A→係数行列（転置行列で入力する）
LDA→AのLeading dimension。とりあえずNと同じ。
IPIV→ピボットのインデックスを格納するベクトル。Bと同じ大きさ。
B→bベクトル
LDB→BのLeading dimension。Nを入れておけば良いらしい。
INFO→結果を格納する整数。

ハマったのは、dgesv()にA行列を転置で入力しないといけない点。
なぜそういう仕様になっているのか。
lapacktest.f95

program call_lapack
  implicit none
  integer :: n =2
  integer info

  ! 静的に読み込む場合
  ! double precision a(2,2), b(2), ipiv(2)
  ! data a /1,2, 3,4/
  ! data b /5,6/

  ! 動的に読み込む場合
  double precision, allocatable :: a(:,:), b(:), ipiv(:)
  allocate(a(n,n))
  allocate(b(n))
  allocate(ipiv(n))
  a = reshape([1,2,3,4],[2,2])
  b = [5,6]

  print *,'A='
  print 10,a
  print *,'B='
  print 20,b
  10 format(2f10.4)
  20 format(f10.4)

  ! ライブラリルーチンの呼び出し
  ! なぜかAを転置行列で入力しないといけない
  call dgesv(n,1,transpose(a),n,ipiv,b,n,info)

  ! 結果を出力
  print *, "Solution"
  print "((f9.4))", b

end program call_lapack

コンパイル

$ gfortran -o lapacktest lapacktest.f95 -lblas -llapack

実行結果

./lapacktest

 A=
    1.0000    2.0000
    3.0000    4.0000
 B=
    5.0000
    6.0000
 Solution
  -4.0000
   4.5000