c     ex0502.f
      program Matrix_inv_11
      implicit none
      integer i, j, k, N, N2, Ns, Ns2
      parameter( N=11, N2=N*2 )
      real*8 A(N,N), E(N,N), D(N,N2), Ainv(N,N), B(N,N), C(N,N)
      real*8 Dkk, Dik
      data E/121*0d0/, A/121*0d0/

      open(1,file='ex0502.dat')
      read(1,*) Ns
      Ns2 = Ns * 2
      do i=1,Ns
       read(1,*) ( A(i,j), j=1, Ns )
      enddo
      close(1)
      write(*,*) 'A'
      call print_Matrix( A, Ns )

      do i=1,N
        E(i,i)=1d0
      enddo
      write(*,*) 'E'
      call print_Matrix( E, Ns )

       do i=1,Ns
         do j=1,Ns
            D(i,j   ) = A(i,j)
            D(i,j+Ns) = E(i,j)
         enddo
      enddo
      write(*,*) 'D'
      call print_Matrix_D( D, Ns )

      k   = 1
      write(*,*) 'k=', k
      Dkk = D(1,1)
      do j = 1, Ns2
         D(1,j) = D(1,j) / Dkk
      enddo
      call print_Matrix_D( D, Ns )
      do i = 2, Ns
         Dik = D(i,1)
         do j = 1, Ns2
            D(i,j) = D(i,j) - Dik * D(1,j)
         enddo
      enddo
      call print_Matrix_D( D, Ns )

      k   = 2
      write(*,*) 'k=', k
      Dkk = D(2,2)
      do j = 1, Ns2
         D(2,j) = D(2,j) / Dkk
      enddo
      call print_Matrix_D( D, Ns )
      do i = 1, Ns
         Dik = D(i,2)
         if( i.ne.2 ) then
            do j = 2, Ns2
               D(i,j) = D(i,j) - Dik * D(2,j)
            enddo
         endif
      enddo
      call print_Matrix_D( D, Ns )

      k   = 3
      write(*,*) 'k=', k
      Dkk = D(3,3)
      do j = 1, Ns2
         D(3,j) = D(3,j) / Dkk
      enddo
      call print_Matrix_D( D, Ns )
      do i = 1, Ns
         Dik = D(i,3)
         if( i.ne.3 ) then
            do j = 3, Ns2
               D(i,j) = D(i,j) - Dik * D(3,j)
            enddo
         endif
      enddo
      write(*,*) 'D'
      call print_Matrix_D( D, Ns )

      do i=1,Ns
         do j=1,Ns
            B(i,j)    = D(i,j   )
            Ainv(i,j) = D(i,j+Ns)
         enddo
      enddo
      write(*,*) 'B'
      call print_Matrix( B,    Ns )
      write(*,*) 'A'
      call print_Matrix( A,    Ns )
      write(*,*) 'Ainv'
      call print_Matrix( Ainv, Ns )

      call multiply_AB_C( A, Ainv, C, Ns )
      write(*,*) 'C = A * Ainv'
      call print_Matrix( C,    Ns )

      stop
      end

      subroutine print_Matrix( A, Ns )
      implicit none
      integer i, j, N, Ns
      parameter( N=11 )
      real*8 A(N,N)
      do i=1,Ns
       write(*,'(1p10e13.5)') ( A(i,j), j=1, Ns )
      enddo
      write(*,*)
      return
      end

      subroutine print_Matrix_D( A, Ns )
      implicit none
      integer i, j, N, N2, Ns, Ns2
      parameter( N=11, N2=N*2 )
      real*8 A(N,N2)
      Ns2 = Ns*2
      do i=1,Ns
c      write(*,'(   3f11.4,4x,3f12.4)') (A(i,j),j=1,Ns2)
       write(*,'(1p10e13.5)')           (A(i,j),j=1,Ns2)
      enddo
      write(*,*)
      return
      end

      subroutine multiply_AB_C( A, B, C, Ns )
      implicit none
      integer i, j, k, N, Ns
      parameter( N=11 )
      real*8 A(N,N), B(N,N), C(N,N)
      do i=1,Ns
       do j=1,Ns
        C(i,j) = 0d0
        do k=1,Ns
         C(i,j) = C(i,j) + A(i,k) * B(k,j)
        enddo
       enddo
      enddo
      return
      end