|
forums.silverfrost.com Welcome to the Silverfrost forums
|
View previous topic :: View next topic |
Author |
Message |
stfark1
Joined: 02 Sep 2008 Posts: 210
|
Posted: Tue Nov 18, 2014 5:26 pm Post subject: 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 |
|
Back to top |
|
|
PaulLaidler Site Admin
Joined: 21 Feb 2005 Posts: 7927 Location: Salford, UK
|
Posted: Tue Nov 18, 2014 6:02 pm Post subject: |
|
|
No. You will need a third party library for this. |
|
Back to top |
|
|
Wilfried Linder
Joined: 14 Nov 2007 Posts: 314 Location: Düsseldorf, Germany
|
Posted: Tue Nov 18, 2014 6:18 pm Post subject: |
|
|
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 |
|
|
Wilfried Linder
Joined: 14 Nov 2007 Posts: 314 Location: Düsseldorf, Germany
|
Posted: Tue Nov 18, 2014 6:19 pm Post subject: |
|
|
... 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 |
|
|
davidb
Joined: 17 Jul 2009 Posts: 560 Location: UK
|
Posted: Tue Nov 18, 2014 7:11 pm Post subject: Re: Matrix Inversion |
|
|
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 |
|
|
Kenneth_Smith
Joined: 18 May 2012 Posts: 697 Location: Hamilton, Lanarkshire, Scotland.
|
Posted: Wed Nov 19, 2014 11:03 pm Post subject: |
|
|
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 |
|
|
|
|
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
|