forums.silverfrost.com
Welcome to the Silverfrost forums

Author Message
stfark1

Joined: 02 Sep 2008
Posts: 86

 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
PaulLaidler
Site Admin

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

 Posted: Tue Nov 18, 2014 6:02 pm    Post subject: No. You will need a third party library for this.
Wilfried Linder

Joined: 14 Nov 2007
Posts: 303
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
Wilfried Linder

Joined: 14 Nov 2007
Posts: 303
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
davidb

Joined: 17 Jul 2009
Posts: 522
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
Kenneth_Smith

Joined: 18 May 2012
Posts: 82
Location: Glasgow, 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
 Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT + 1 Hour Page 1 of 1

 Jump to: Select a forum Admin----------------Announcements FTN95----------------GeneralKBaseSupportSuggestionsClearWin+Plato64-bit FTN77----------------Support
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