forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Matrix Inversion

 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Plato
View previous topic :: View next topic  
Author Message
stfark1



Joined: 02 Sep 2008
Posts: 210

PostPosted: Tue Nov 18, 2014 5:26 pm    Post subject: Matrix Inversion Reply with quote

Wondering if there is a routine in the Silverfrost Fortran library that will inverse a 3d matrix, please let me know and how to access, Sid Kraft
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 7916
Location: Salford, UK

PostPosted: Tue Nov 18, 2014 6:02 pm    Post subject: Reply with quote

No. You will need a third party library for this.
Back to top
View user's profile Send private message AIM Address
Wilfried Linder



Joined: 14 Nov 2007
Posts: 314
Location: Düsseldorf, Germany

PostPosted: Tue Nov 18, 2014 6:18 pm    Post subject: Reply with quote

Here is an example:

Code:
subroutine invert(mdim,A,rtcode)

IMPLICIT NONE

integer*4  mdim,zdim,i,j,k,rtcode
real*8     A(1:50,1:50),B(1:50,1:100)
real*8     r,fastnix

rtcode  = -1
if (mdim > 50) return
fastnix = 10.D0*tiny(1.D0)
zdim    = 2*mdim

! B = input matrix plus identity matrix

do k = 1,mdim
  do j = 1,mdim
    B(k,j) = A(k,j)
    if (k == j) then
      B(k,j+mdim) = 1.D0
    else
      B(k,j+mdim) = 0.D0
    end if
  end do
end do

do k = 1,mdim

! if diagonal element equal 0, change rows

  if (B(k,k) < fastnix) then
    do i = k+1,mdim
      if (B(i,k) > fastnix) then
        do j = 1,zdim
          r      = B(k,j)
          B(k,j) = B(i,j)
          B(i,j) = r
        end do
        exit
      end if
    end do
  end if

! transform to upper triangle matrix

  do i = k,mdim
    if (abs(B(i,k)) > fastnix) then
      r = B(i,k)
      do j = 1,zdim
        B(i,j) = B(i,j)/r
      end do
      if (i > k) then
        do j = 1,zdim
          B(i,j) = B(i,j)-B(k,j)
        end do
      end if
    end if
  end do
end do

! diagonal element equal 0 --> return

do k = 1,mdim
  if (abs(B(k,k)) < fastnix) return
end do

! transform to diagonal matrix

do k = mdim-1,1,-1
  do i = k,1,-1
    r = B(i,k+1)
    do j = 1,zdim
      B(i,j) = B(i,j)-r*B(k+1,j)
    end do
  end do
end do

! copy result to matrix A

do k = 1,mdim
  do j = 1,mdim
    A(k,j) = B(k,j+mdim)
  end do
end do

rtcode = 0
return
end

Back to top
View user's profile Send private message
Wilfried Linder



Joined: 14 Nov 2007
Posts: 314
Location: Düsseldorf, Germany

PostPosted: Tue Nov 18, 2014 6:19 pm    Post subject: Reply with quote

... and here a test program:

Code:
WINAPP
PROGRAM MAT_INV

IMPLICIT NONE

! ========================================================= !
!
! Matrix inversion, method from Gauss
!
! ========================================================= !

integer*4      mdim,ctrl,i,j,k,rtcode
real*8         A(1:50,1:50),B(1:50,1:50),P(1:50,1:50)
real*8         random,dr
character*500  string

mdim = 9
dr = 5.D0

j = winio@('%ca[Invert matrix]%sy[3d_thin]&')
j = winio@('Dimension  %dd%il%3rd%ff&',1,1,50,mdim)
j = winio@('Entries from intervall +/- %2rf%ff&',dr)
j = winio@('%nl%cn%`7bt[OK]')

j = winio@('%lw%sy[3d_thin]%fn[Arial]%ts&',ctrl,1.D0)
j = winio@('%120.40cw[vscroll,hscroll]',6)

dr = max(dr,1.D0)

! set co-efficients with random generator:

call date_time_seed@()

write(6,'(/,A,/)')'Start matrix:'
do i = 1,mdim
  string = ' '
  do j = 1,mdim
    A(i,j) = dr*random()-dr
    k = 10*(j-1)+1
    write(string(k:k+9),'(F10.3)')A(i,j)
  end do
  write(6,'(A)')string
end do
B = A

! invert matrix

call invert(mdim,A,rtcode)

if (rtcode < 0) then
  write(6,'(/,A)')'*** Rank defect ***'
else

! result:

  write(6,'(/,A,/)')'Invers matrix:'
  do i = 1,mdim
    string = ' '
    do j = 1,mdim
      k = 10*(j-1)+1
      write(string(k:k+9),'(F10.3)')A(i,j)
    end do
    write(6,'(A)')string
  end do

! test:

  write(6,'(/,A,/)')'Test:'
  do i = 1,mdim
    string = ' '
    do j = 1,mdim
      P(i,j) = 0.D0
      do k = 1,mdim
        P(i,j) = P(i,j)+B(i,k)*A(k,j)
      end do
      k = 10*(j-1)+1
      write(string(k:k+9),'(F10.3)')P(i,j)
    end do
    write(6,'(A)')string
  end do

end if
end

Wilfried
Back to top
View user's profile Send private message
davidb



Joined: 17 Jul 2009
Posts: 560
Location: UK

PostPosted: Tue Nov 18, 2014 7:11 pm    Post subject: Re: Matrix Inversion Reply with quote

stfark1 wrote:
Wondering if there is a routine in the Silverfrost Fortran library that will inverse a 3d matrix, please let me know and how to access, Sid Kraft


What is a 3d matrix?

If just a matrix then try Wilfried's code above.
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl
Back to top
View user's profile Send private message
Kenneth_Smith



Joined: 18 May 2012
Posts: 697
Location: Hamilton, Lanarkshire, Scotland.

PostPosted: Wed Nov 19, 2014 11:03 pm    Post subject: Reply with quote

An alternative for a 2d matrix is that by Shipley and Coleman. It's not well know, but I've found it very useful for large complex matrices. Here's a very simple example.

Cheers Ken


Code:
     
!--------------------------------------------------------------------72
!
!     SUBROUTINE INV_SHIPLEY_COLEMAN              KSS 5th August 2014
!
!--------------------------------------------------------------------72     
!     Subroutine  INV_SHIPLEY_COLEMAN takes the input matrix Y(N,N)
!     and returns the inverse matrix in Y using the Shipley-Coleman
!     Inversion method.
!
!     Reference
!     =========
!     R. B. Shipley and D. Coleman  "A New Direct Matrix Inversion
!     Method",  AIEE Trans. (Communication and Electronics),  vol. 78,
!     pp.568 -572 1959
!--------------------------------------------------------------------72
      SUBROUTINE INV_SHIPLEY_COLEMAN (Y,N)
      IMPLICIT NONE
      REAL*8, INTENT(INOUT) :: Y(N,N)
      INTEGER,   INTENT(IN)  :: N
      INTEGER J, K, M
      REAL*8, PARAMETER:: ONE = 1.d0
      REAL*8, PARAMETER:: ZERO = 0.d0

      DO M = 1, N, 1
       
          IF ( ABS( Y(M,M) ) .EQ. ZERO) THEN
              WRITE(6,*)
              WRITE(6,*) 'ERROR in subroutine INV_SHIPLEY_COLEMAN'
              WRITE(6,*) 'One or more main diagonal elements is zero'
              STOP
          END IF
         
      END DO
     
!     Do for each pivot m = 1,2,3, ... ,n
     DO M  = 1, N, 1
     
!         Kron reduce all rows and columns, except the pivot row and
!         pivot column
          DO J  = 1, N, 1
             IF (J .NE. M) THEN
                  DO K  = 1, N, 1
                     IF (K .NE. M) THEN
                          Y(J,K) = Y(J,K) - Y(J,M) * Y(M,K) / Y(M,M)
                      END IF
                  END DO
              END IF
         END DO
         
C         Replace the pivot with its negative inverse
         Y(M,M)   = -ONE / Y(M,M)
         
!         Work across the pivot row and down the pivot column, multiplying
!         by the new pivot value
         DO K  = 1, N, 1
             IF (K .NE. M) THEN
                 Y(M,K) = Y(M,K) * Y(M,M)
                 Y(K,M) = Y(K,M) * Y(M,M)
              END IF
          END DO
         
      END DO
     
C     Negate the results
     Y = -Y
     RETURN
      END SUBROUTINE INV_SHIPLEY_COLEMAN
 
Back to top
View user's profile Send private message Visit poster's website
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Plato All times are GMT + 1 Hour
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group