forums.silverfrost.com
Welcome to the Silverfrost forums

Author Message
stfark1

Joined: 02 Sep 2008
Posts: 115

 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

Joined: 21 Feb 2005
Posts: 5353
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: 553
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: 175
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