간단히 역행렬 구하기(포트란에서 슬라텍이용), 또는 다른 라이

소스가 오픈되어 있는 수치해석용 라이브리 중에서 슬라텍으로 간단히 테스트를 해본 것입니다. 그런데 답이 영 엉뚱하게 나옵니다. 아마도 제가 잘 못 알고 있는 것 같은데, 해결 방법이 없을까요?

program SGEDI_test
  implicit none
  real :: a(3,3)
  integer :: ipvt(3)
  integer :: i, j
  real, dimension(3) :: work
  real, dimension(2) :: det
  
  open(1, file='a',  status='old')
  open(2, file='aa', status='old')
  open(3, file='b',  status='old')
  open(4, file='bb', status='old')
  open(5, file='inva', status='unknown')
  
  do i =1,3
     read(1,*) a(i,1:3)        
  enddo

  write(*,*) 'Matrix a'
  do i=1,3
     write(*,*) a(i,1:3)

  enddo

  call SGEDI(a, 3, 3, ipvt, det, work, 11)
  
  write(*,*) 'Invert Matrix a'
  do i=1,3
     write(*,*) a(i,1:3)
     write(5,*) a(i,1:3)
  enddo
  
  write(*,*)  det
  close(1); close(2); close(3); close(4); close(5)
end program SGEDI_test

slatec에 있는 SGEDI의 설명문 입니다. 아래 부분의 소스는 생략을 했습니다.

지금 가장 이해가 안되는 부분은 N이라는 변수의 역할입니다. 그냥 보기에는 A행렬의 열을 넣으면 되는 것 같은데 그게 아닌 모양입니다.

subroutine SGEDI (A, LDA, N, IPVT, DET, WORK, JOB)
!
!! SGEDI computes the determinant and inverse of a matrix using the ...
!            factors computed by SGECO or SGEFA.
!
!***LIBRARY   SLATEC (LINPACK)
!***CATEGORY  D2A1, D3A1
!***TYPE      SINGLE PRECISION (SGEDI-S, DGEDI-D, CGEDI-C)
!***KEYWORDS  DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX
!***AUTHOR  Moler, C. B., (U. of New Mexico)
!***DESCRIPTION
!
!     SGEDI computes the determinant and inverse of a matrix
!     using the factors computed by SGECO or SGEFA.
!
!     On Entry
!
!        A       REAL(LDA, N)
!                the output from SGECO or SGEFA.
!
!        LDA     INTEGER
!                the leading dimension of the array  A .
!
!        N       INTEGER
!                the order of the matrix  A .
!
!        IPVT    INTEGER(N)
!                the pivot vector from SGECO or SGEFA.
!
!        WORK    REAL(N)
!                work vector.  Contents destroyed.
!
!        JOB     INTEGER
!                = 11   both determinant and inverse.
!                = 01   inverse only.
!                = 10   determinant only.
!
!     On Return
!
!        A       inverse of original matrix if requested.
!                Otherwise unchanged.
!
!        DET     REAL(2)
!                determinant of original matrix if requested.
!                Otherwise not referenced.
!                Determinant = DET(1) * 10.0**DET(2)
!                with  1.0  <=  ABS(DET(1))  <  10.0
!                or  DET(1)  ==  0.0 .
!
!     Error Condition
!
!        A division by zero will occur if the input factor contains
!        a zero on the diagonal and the inverse is requested.
!        It will not occur if the subroutines are called correctly
!        and if SGECO has set RCOND  >  0.0 or SGEFA has set
!        INFO  ==  0 .
!
!***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
!                 Stewart, LINPACK Users' Guide, SIAM, 1979.
!***ROUTINES CALLED  SAXPY, SSCAL, SSWAP
!***REVISION HISTORY  (YYMMDD)
!   780814  DATE WRITTEN
!   890831  Modified array declarations.  (WRB)
!   890831  REVISION DATE from Version 3.2
!   891214  Prologue converted to Version 4.0 format.  (BAB)
!   900326  Removed duplicate information from DESCRIPTION section.
!           (WRB)
!   920501  Reformatted the REFERENCES section.  (WRB)
!***END PROLOGUE  SGEDI
  INTEGER LDA,N,IPVT(*),JOB
  REAL A(LDA,*),DET(2),WORK(*)
!
  REAL T
  REAL TEN
  INTEGER I,J,K,KB,KP1,L,NM1
!***FIRST EXECUTABLE STATEMENT  SGEDI
!
!     COMPUTE DETERMINANT

간단히 옥타브로 계산을 해 보면 얻고자 하는 답은 아래와 같습니다.

octave:1> load a
octave:2> a
a =

   2   0   1
  -2   3   4
  -5   5   6

octave:3> inv (a)
ans =

   -2.0000    5.0000   -3.0000
   -8.0000   17.0000  -10.0000
    5.0000  -10.0000    6.0000

octave:4> a* inv(a)
ans =

   1.00000   0.00000   0.00000
  -0.00000   1.00000  -0.00000
  -0.00000   0.00000   1.00000

어떻게 하면 될까요? 지금까지는 뉴메리컬레서피에 나오는 가우스 조던을 그냥 써 왔는데 저도 다른 방향으로 한번 가볼까 하는데 벌써 이틀동안 진도가 전혀 안나가네요. :oops:

혹시 뉴메릭 라이브러리를 이용하시는 분이 계시면 조언을 바랍니다.