Understanding DGETRF/SGETRF interface

I was working on the inverse matrix the other day, it can be done through DGETRI of LAPACK. The DGETRI requires DGETRF to do LU decomposition first in order to calculate the inverse matrix. Then I found DGETRF interesting to start with since it has been awhile for me to touch linear algebra stuff. I wrote a "hello world" code for it. I was stuck on the issue to construct permutation matrix from IPIV, basically IPIV is not like vecP in Mathematica. This is how I understand it, I am not confident about the P matrix, this requires some LAPACK friends' help. On Mac OS X, just do g95 -framework Accelerate hello.f90 -o hellolapack and run the program by ./hellolapack , sorry I am too lazy to install a software lapack. :p

program hello
   integer INFO
   integer IPIV(3)
   real*8 A(3,3),B(3,3),C(3,3),D(3,3)
   real*8 P(3,3),L(3,3),U(3,3)
   real*8 myVEC(3)
   real*8 WORK(64*3)
   A(1,1)=1d0
   A(1,2)=2d0
   A(1,3)=2d0
   A(2,1)=2d0
   A(2,2)=7d0
   A(2,3)=7d0
   A(3,1)=2d0
   A(3,2)=7d0
   A(3,3)=9d0
   P=0d0
   P(1,1)=1d0
   P(2,2)=1d0
   P(3,3)=1d0
   ! dgetrf(M,N,A,LDA,IPIV,INFO)
   ! input
   ! M          The number of rows of the matrix A.  M >= 0. Integer
   ! N          The number of columns of the matrix A.  N >= 0. Integer
   ! A          The DP array
   ! output
   ! A          The DP array
   ! LDA        The leading dimension of the array A.  LDA >= max(1,M).
   ! IPIV       INTEGER array, dimension (min(M,N))
   !            The pivot indices; for 1 <= i <= min(M,N), row i of the
   !            matrix was interchanged with row IPIV(i).
   print *,"A"
   write (*,'(F8.3,F8.3,F8.3)') A(1,1:3)
   write (*,'(F8.3,F8.3,F8.3)') A(2,1:3)
   write (*,'(F8.3,F8.3,F8.3)') A(3,1:3)
   B=A
   call dgetrf(3,3,A,3,IPIV,INFO)
   if (INFO.ne.0) stop
   !IPIV(1)= 3 implies that the first and third rows were interchanged
   !when factoring the first column; IPIV(2)= 3 implies that the second
   !and third rows were interchanged when factoring the second column.
   !In this case, IPIV(3) must be 3 because there are only three rows.
   if (IPIV(1).eq.3) then
      myVEC=P(1,1:3)
      P(1,1:3)=P(3,1:3)
      P(3,1:3)=myVEC
   elseif(IPIV(1).eq.2) then
      myVEC=P(1,1:3)
      P(1,1:3)=P(2,1:3)
      P(2,1:3)=myVEC
   endif
   if (IPIV(2).eq.3) then
      myVEC=P(2,1:3)
      P(2,1:3)=P(3,1:3)
      P(3,1:3)=myVEC
   elseif(IPIV(2).eq.1) then
      myVEC=P(2,1:3)
      P(2,1:3)=P(1,1:3)
      P(1,1:3)=myVEC
   endif
   print *,"P"
   write (*,'(F8.3,F8.3,F8.3)') P(1,1:3)
   write (*,'(F8.3,F8.3,F8.3)') P(2,1:3)
   write (*,'(F8.3,F8.3,F8.3)') P(3,1:3)
   print *,"L"
   L(1,1)=1d0
   L(1,2:3)=0d0
   L(2,1)=A(2,1)
   L(2,2)=1d0
   L(2,3)=0d0
   L(3,1:2)=A(3,1:2)
   L(3,3)=1d0
   write (*,'(F8.3,F8.3,F8.3)') L(1,1:3)
   write (*,'(F8.3,F8.3,F8.3)') L(2,1:3)
   write (*,'(F8.3,F8.3,F8.3)') L(3,1:3)
   print *,"U"
   U(1,1:3)=A(1,1:3)
   U(2,1)=0d0
   U(2,2:3)=A(2,2:3)
   U(3,1:2)=0d0
   U(3,3)=A(3,3)
   write (*,'(F8.3,F8.3,F8.3)') U(1,1:3)
   write (*,'(F8.3,F8.3,F8.3)') U(2,1:3)
   write (*,'(F8.3,F8.3,F8.3)') U(3,1:3)
   call dgemm("N","N",3,3,3,1d0,L,3,U,3,0d0,C,3)
   print *,"L.U"
   write (*,'(F8.3,F8.3,F8.3)') C(1,1:3)
   write (*,'(F8.3,F8.3,F8.3)') C(2,1:3)
   write (*,'(F8.3,F8.3,F8.3)') C(3,1:3)
   call dgemm("N","N",3,3,3,1d0,P,3,L,3,0d0,C,3)
   call dgemm("N","N",3,3,3,1d0,C,3,U,3,0d0,D,3)
   print *,"P.L.U"
   write (*,'(F8.3,F8.3,F8.3)') D(1,1:3)
   write (*,'(F8.3,F8.3,F8.3)') D(2,1:3)
   write (*,'(F8.3,F8.3,F8.3)') D(3,1:3)

   call dgetri(3,A,3,IPIV,WORK,64*3,INFO)
   if (INFO.ne.0) stop
   call dgemm("N","N",3,3,3,1d0,B,3,A,3,0d0,C,3)
   print *,"A.A^-1"
   write (*,'(F8.3,F8.3,F8.3)') C(1,1:3)
   write (*,'(F8.3,F8.3,F8.3)') C(2,1:3)
   write (*,'(F8.3,F8.3,F8.3)') C(3,1:3)
end program hello