Matrix Inversion

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

No. You will need a third party library for this.

Here is an example:

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

... and here a test program:

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

Quoted from stfark1 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.

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

!--------------------------------------------------------------------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