netlib.org/lapack참조. Lapack(라팍)은 선형 시스템, 선형 시스템의 최소제곱법, 고유값문제, singular value문제 등을 해결하는 routine들을 제공한다.
Numpy, Matlab, Maple, R등에서도 이미 쓰이고, mac os x에도 쓰이고 fortran77로 쓰여있다.
lapack 사이트에서 보면 내가 원하는 subroutine을 찾을수 있다. 참고로 난 이 documentation이 매우 마음에 들지 않는다.
All driver/ computational routine들은 XYYZZZ라는 모양의 이름을 가지고있다. 6번째가 blank이기도 하다.
첫번째 X는 data type을 알려준다:
S : Real
D : Double Precision
C : Complex
Z : Complex*16 or Double Complex
그다음의 YY는 type of matrix를 알려준다. 참고로 여긴 몇가지만 적겠지만 총 28가지의 matrix종류가 있다.
BD : Bidiagonal
DI : Diagonal
GB : General Band
GE : General (i.e. unsymmetric, in some cases rectangular)
GG : General matrices, generalized problem(i.e. a pair of general matrices)
GT : General Tridiagonal
HB : (complex) Hermitian band
HE : (complex) Hermitian
그래서 예를들어, DGESV이면 , Double Precision routine for SolVing(SV) systems with general matrices라는 뜻.
내가 이번 과제에서 사용했던 DSYEV 는 Double precision routine for computing Eigenvalues and optionally eigenvectors of a real symmetric matrix A. 이다.
그냥, 실대칭행렬의 고유값과 고유벡터를 계산해주는 routine.
근데 너무 힘들었던게 ㅋㅋㅋㅋㅋ
matlab에선 그냥 eig(M) 이러면 고유값 나오지않나여..?
얘는 documentation이 너무 복잡해서 힘들었음.
일단 보면 알겠지만 매우 복잡함. (나만 그런가? 처음써봐서?) 설명도 별로 친절하지가않음.
subroutine dsyev를 사용하려면 (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO) 를 적어줘야 한다.
작명최악…–;
하나씩 풀어서 설명하고 예제를 보여줄 것이다.
JOBZ : 는 인풋이다. character*1이다. ‘N’또는 ‘V’를 적어주면 된다. 각각 고윳값만 구할것인지, 고윳값과 고유벡터 둘다 구할것인지 먼저 셋팅해주는 것이다.
UPLO : 는 인풋이다. character*1이다. ‘U’ 를 적어주면 A의 상삼각행렬이 저장되고, ‘L’을 적어주면 A의 하삼각행렬이 저장된다.
N : 는 integer 인풋이다. A의 order를 적어주면 되고, 0보다 같거나 커야 한다.
A : 인풋이자 아웃풋이다. 우리가 구하고자하는 행렬, double precision array로 인풋을 주면 된다. dimension은 (LDA,N)이다. 고유벡터 아웃풋을 이 행렬에 다시 담아서 돌려주는 것 같다.
LDA : 인풋이고 integer이다. Leading dimension of A이다.
W : double precision array이고 아웃풋이다. dimension (N)이다. Info값이 0이면 고윳값이 ascending order로 여기에 담겨서 리턴된다.
WORK : double precision array이고 아웃풋이다. dimension은 (Max(1,LWORK))이다. subroutine이 끝날 때, INFO 값이 0이면, WORK(1)은 optimal LWORK를 리턴해준다.
LWORK : Integer이고 input이다. WORK어레이의 length이다. LWORK는 >= max(1,3*N-1)이라는 조건을 만족시켜주어야한다.
Optimal efficiency를 위해서는 LWORK >= (NB + 2)*N 이면 된다. NB는 blocksize for DSYTRD by ILAENV라는데 모르겠다 뭔말인지….
하지만 중요한건! LWORK = -1이면 workspace query를 가정하고 돌린다. 그러면 routine은 optimal size of WORK array만 계산해준다. 그리고 WORK(1)값에 그 값이 담긴다. 에러메세지도 생기지 않는다. 이게 무슨말인지 잠시후 예제에서 다룰 것이다.
INFO : Integer 이고 아웃풋인데, 쉽게말해선 이게 제대로 돌았는지 아닌지를 보여준다. 0이나오면 제대로된것이고, 음수가 나오면 , -i가 나오면 ith argument가 illegal value를 가졌다는 뜻이다. 양수로 i가 나오면, algorithm이 수렴하지 못했다는 것이다. I off-diagonal elements of an intermediate tridiagonal form 이 0으로 수렴하지 않았다는 뜻이다.
일단나는 이 documentation보다더 도움이됬던건 인텔 문서였다.
https://software.intel.com/en-us/node/469176
여기에 보면 dsyev에대해 다시 한번 설명을 해준다.
그리고 인텔에서 제공하는 예제코드가 있는데 난 이게 되게 많이 도움이 됐다.
예제 코드를 보면 dsyev를 두번 call하는걸 볼 수 있다.
처음에 lwork의 optimal값을 모르기때문에 -1로 놓고 돌리고,
거기서 리턴된 work(1)값을 lwork에 다시 담아주고, work를 deallocate해준 후에 새로운 lwork값으로 dsyev를 돌려주면 우리가 원하는 고유값을 W에서 구할수 있다 !
