Inverseion of a Matrix

[mshokrey]

Dear friend,

Can any one of you pass me a subroutine for inversion of a matrix which support F90.

I know there are many libraries but I need a simple subroutine.

Thanks

 

[xwb]

Can provide you with a matrix inversion routine but what are you using it for?

If it is for solving simultaneous equations, Choleski LU decomposition would be a better technique.

 

[mshokrey]

No it is not for solving simultaneous equation. I just need a routine which estimate the inverse of a matrix for normal matrix operations as multiplication ...etc.

 

[xwb]

This isn't the most efficient method in the world and it assumes you do not have any ill conditioning.

module ModMatrix
   type Matrix
      character*16:: name
      integer:: xmax
      integer:: ymax
      real, allocatable, dimension(:,:):: d
   end type Matrix
contains
   ! Create a square matrix
   subroutine MatCreateSquare (mat, cols, tag)
      type (Matrix), intent(inout):: mat
      integer, intent(in):: cols
      character*(*), intent(in):: tag

      mat%name = tag
      mat%xmax = cols
      mat%ymax = cols
      allocate (mat%d(mat%xmax, mat%ymax))
   end subroutine MatCreateSquare
   
   ! Create an identity matrix
   subroutine MatIdent (mat)
      type (Matrix), intent(inout):: mat
      integer:: xix, yix

      do yix = 1, mat%ymax
         do xix = 1, mat%xmax
            mat%d(xix,yix) = 0.0
         enddo
      enddo

      do xix = 1, mat%xmax
         mat%d(xix,xix) = 1.0
      enddo
   end subroutine
   
   ! Copy cpy to mat
   subroutine MatCopy (mat, cpy)
      type (Matrix), intent(inout):: mat
      type (Matrix), intent(in):: cpy
      integer:: xix, yix
   
      do yix = 1, cpy%ymax
         do xix = 1, cpy%xmax
            mat%d(xix,yix) = cpy%d(xix,yix)
         enddo
      enddo
   end subroutine
   
   ! Invert a matrix
   subroutine MatInv (mat, inv)
      type (Matrix), intent(in):: mat
      type (Matrix), intent(inout):: inv
      type (Matrix):: wkg
      integer:: pix, xix, yix, ix
      real:: mul

      call MatIdent (inv)
      call MatCreateSquare (wkg, mat%xmax, 'wkg')
      call MatCopy (wkg, mat)

      do pix = 1, inv%ymax
         ! Convert row to unit
         mul = wkg%d(pix,pix)
         if (mul .eq. 0) then
             write (*, *) 'Unable to invert matrix'
             return
         end if
         mul = 1.0 / mul
         
         ! Create a pivot row
         do xix = 1, wkg%xmax
            wkg%d(xix,pix) = wkg%d(xix,pix) * mul
            inv%d(xix,pix) = inv%d(xix,pix) * mul
         enddo

         ! wkg will eventually become identity matrix
         ! inv will become the inverse
         do yix = 1, inv%ymax
            if (yix .ne. pix) then
               mul = wkg%d(pix,yix)
               do xix = 1, inv%xmax
                  wkg%d(xix,yix) = wkg%d(xix,yix) - wkg%d(xix,pix) * mul
                  inv%d(xix,yix) = inv%d(xix,yix) - inv%d(xix,pix) * mul
               enddo
               ! call MatDump (wkg)
               ! call MatDump (inv)
            end if
         enddo
      enddo
   end subroutine
   
   ! Dump the matrix to the screen
   subroutine MatDump (mat)
      type (Matrix), intent(in):: mat
      integer:: xix, yix
      character*16:: fmt

      write (*, *) mat%name
      ! this will work for matrices from 1 to 9
      write (fmt, 10) mat%xmax
10    format ('(',I1, 'F7.2)')
      do yix = 1, mat%ymax
         write (*, fmt) (mat%d(xix,yix), xix = 1, mat%xmax)
      enddo
   end subroutine

   ! Fill a matrix from a sequential array
   subroutine MatFill (mat, info)
      type (Matrix), intent(inout):: mat
      real, dimension(:), intent (in):: info
      integer:: ix, xix, yix
      
      ix = 0
      do yix = 1, mat%ymax
         do xix = 1, mat%xmax
            ix = ix + 1
            mat%d(xix,yix) = info(ix)
         enddo
      enddo
   end subroutine MatFill
end module

program main
   use ModMatrix
   type (Matrix):: lft, rgt, sng, inv
   real, dimension(9):: lhs, rhs
   data lhs / 1., 3., 3., 1., 4., 3., 1., 3., 4. /
   data rhs / 7., -3., -3., -1., 1., 0., -1., 0., 1. /

   call MatCreateSquare (lft, 3, 'lft')
   call MatFill (lft, lhs)
   call MatDump (lft)

   call MatCreateSquare (inv, 3, 'inv')
   call MatInv (lft, inv)
   call MatDump (inv)
   stop
end program

 

[mshokrey]

Thanks for help,

I really don't know much of Matrices, but I faced a formula in my code that need matrix inversion. I'll try to use this code, If I have any forward problems, I'll contact you later.

Again Thanks for your help