Fortran program for calculating matrices inverses with jacobi algorithm

[jomiro87]

I am trying to set a program in fortran 95 for calculating inverses of matrices by using the jacobi algorithm. I have compiled it successfully, but I execute it, it is only able to calculate correctly the first column of the new matrix. The code is:

Program algo
implicit none
real(8), dimension,, allocatable :: A, ID, AINVI, AS, D, DINV, PM, SM, AINVF
real(8) :: sub, e, tru, trd, FindDet
integer :: i, j, n, iter
print *, 'Dimensió matriu'
read *, n
allocate (A(n,n))
allocate (ID(n,n))
allocate (AINVI(n,n))
allocate (AS(n,n))
allocate (D(n,n))
allocate (DINV(n,n))
allocate (PM(n,n))
allocate (SM(n,n))
allocate (AINVF(n,n))
print *, 'matriu'
do i=1, n
print *, 'fila', i
read *, (A(i,j), j=1, n)
end do
ID=0.0d0
do i=1, n
ID(i,i)=1.0d0
end do
AINVI=0.0d0
do i=1, n
do j=1, n
if (i == j) then
AS(i,j)=0.0d0
else
AS(i,j)=A(i,j)
end if
end do
end do
do i=1, n
do j=1, n
if (i==j) then
D(i,j)=A(i,j)
else
D(i,j)=0.0d0
end if
end do
end do
do i=1, n
do j=1, n
if (i==j) then
DINV(i,j)=1/D(i,j)
else
DINV(i,j)=0.0d0
end if
end do
end do
do i=1, n
print *, (AS(i,j), j=1, n)
end do
do i=1, n
print *, (D(i,j), j=1, n)
end do
do i=1, n
print *, (DINV(i,j), j=1, n)
end do
PM=matmul(DINV, ID)
SM=matmul(DINV, AS)
print *, 'pm'
do i=1, n
print *, (PM(i,j), j=1, n)
end do
print *, 'sm'
do i=1, n
print *, (SM(i,j), j=1, n)
end do
iter=0
e=0.0001
sub=1.0d0
tru=0.0d0
trd=0.0d0
do while (sub>=e)
iter=iter+1
AINVF=PM-matmul(SM, AINVI)
do i=1, n
tru=tru+AINVI(i,i)
trd=trd+AINVF(i,i)
end do
print *, tru, trd, AINVI(1,1), AINVF(1,1), FindDet(AINVI,n), FindDet(AINVF,n)
sub=abs(tru-trd)
AINVI=AINVF
end do
print *, 'inversa'
do i=1, n
print *, (AINVF(i,j), j=1, n)
end do


print *, 'iteracions', iter
print *, sub
print *, e
end program

REAL(8) FUNCTION FindDet(matrix, n)
IMPLICIT NONE
REAL(8), DIMENSION(n,n) :: matrix
INTEGER, INTENT(IN) :: n
REAL(8) :: m, temp
INTEGER :: i, j, k, l
LOGICAL :: DetExists = .TRUE.
l = 1
!Convert to upper triangular form
DO k = 1, n-1
IF (matrix(k,k) == 0) THEN
DetExists = .FALSE.
DO i = k+1, n
IF (matrix(i,k) /= 0) THEN
DO j = 1, n
temp = matrix(i,j)
matrix(i,j)= matrix(k,j)
matrix(k,j) = temp
END DO
DetExists = .TRUE.
l=-l
EXIT
ENDIF
END DO
IF (DetExists .EQV. .FALSE.) THEN
FindDet = 0
return
END IF
ENDIF
DO j = k+1, n
m = matrix(j,k)/matrix(k,k)
DO i = k+1, n
matrix(j,i) = matrix(j,i) - m*matrix(k,i)
END DO
END DO
END DO

!Calculate determinant by finding product of diagonal elements
FindDet = l
DO i = 1, n
FindDet = FindDet * matrix(i,i)
END DO

END FUNCTION FindDet

For example, for matrix

2 1
3 2

it shows

1.9999999999 -Infinity
-2.9999999999 NaN

When actually, it should give the following:

2 -1
-3 2

I would be grateful if you could provide me any suggestion for improving my code.

Cheers,

Joan

 

[gummibaer]

jomiro

(1) If you post code, then include it in code-tags so the text format of your code is maintained and your code is readable. Put the codetags first and then copy and paste your code between them.

(2) Why reinvent the wheel? It will not get rounder from it. Inverting matrices is a standard procedure since computer-time begun and there are a number of utilities available to do the job. Google for 'fortran libraries' and you will find a lot of math-libraries to pick the solution you like, optimzed for memory usage or speed. This will be a much faster approach than falling into all the pits a lot of programmers have fallen in before.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[FJacq]

Sorry but I don't think that the Jacobi algorithm is used to inverse matrices.

As far as I remember, this is an iterative algorithm used to solve a linear system when the matrix exhibits a dominant diagonal.

But you perhaps confused with the Jordan algorithm ? Here is the simplest example, without looking for the best pivot and without managing possible division by zero :

PROGRAM test
  IMPLICIT NONE
  INTEGER,PARAMETER:: n=5
  DOUBLE PRECISION :: a(n,n),ainv(n,n),z
  INTEGER :: i,j,k
  DO i=1,n
    DO j=1,n
      CALL random_number(a(i,j))
    ENDDO
  ENDDO
  ainv=a
  CALL invers(ainv,n)
  DO i=1,n
    DO j=1,n
      z=0
      DO k=1,n
        z=z+a(i,k)*ainv(k,j)
      ENDDO
      WRITE(*,*) i,j,z
    ENDDO
  ENDDO
END PROGRAM

SUBROUTINE invers(a,n)

  IMPLICIT NONE
  
  INTEGER         ,INTENT(in)    :: n
  DOUBLE PRECISION,INTENT(inout) :: a(n,n)
  
  INTEGER :: k,i,j
  DOUBLE PRECISION :: p,sp
  
  DO k=1,n
    p=1/a(k,k)
    DO j=1,N
      a(k,j)=a(k,j)*P
    ENDDO
    a(k,k)=p
    DO i=1,n
      IF(i == k) CYCLE
      sp=-a(i,k)
      DO j=1,n
        a(i,j)=a(i,j)+sp*a(k,j)
      ENDDO
      a(i,k)=sp*p
    ENDDO
  ENDDO
  
END SUBROUTINE

François Jacq

 

[jomiro87]

I knew that. In fact, I was trying to use Jacobi algorithm by dividing the identity matrix in subarrays composed by its columns like that:


program inverse
real(8), dimension,, allocatable :: a, ID
real(8), dimension), allocatable :: b, x, y
real(8) :: s1, su, sm
integer :: i, j, k, n
print *, 'dimensió'
read *, n
allocate (a(n,n))
allocate (ID(n,n))
allocate (b)
allocate (x)
allocate (y)
print *, 'Introduïu matriu per columnes'
do j=1, n
read *, (A(i,j), i=1, n)
end do
ID=0.0d0
do i=1, n
do j=1, n
if (i == j) ID(i,j)=1.0d0
end do
end do


invc = 0.0d0
sub = 1.0d0
sm =1.0d0
x =0.0d0
print *, '==========='
do while (sm.ge.0.001)
s1 = 0.0
do j=1, n
b=ID,j)

do i = 1,n
su = b(i)


do k = 1,i-1
su = su-a(i,k)*x(k)
end do
do k = i+1,n
su = su-a(i,k)*x(k)
end do
y(i) = su/a(i,i)
s1 = amax1(s1,abs(y(i)))
end do
sm = 0.0
do i = 1,n
sm = amax1(abs(x(i)-y(i))/s1,sm)
x(i) = y(i)
end do
print *,x
end do

But it doesn't work properly.

PS: how can I put codetags?

 

[jomiro87]

I have just done it. Thank you.

rogram inverse
real(8), dimension,, allocatable :: a, ID
real(8), dimension), allocatable :: b, x, y
real(8) :: s1, su, sm
integer :: i, j, k, n
print *, 'dimensió'
read *, n
allocate (a(n,n))
allocate (ID(n,n))
allocate (b)
allocate (x)
allocate (y)
print *, 'Introduïu matriu per columnes'
do j=1, n
read *, (A(i,j), i=1, n)
end do
ID=0.0d0
do i=1, n
do j=1, n
if (i == j) ID(i,j)=1.0d0
end do
end do



x =0.0d0
print *, '==========='
do j=1, n
sm=1.0d0
b=ID,j)
print *, '?????????????'
do while (sm.ge.0.001)
s1 = 0.0
do i = 1,n
su = b(i)
do k = 1,i-1
su = su-a(i,k)*x(k)
end do
do k = i+1,n
su = su-a(i,k)*x(k)
end do
y(i) = su/a(i,i)
s1 = amax1(s1,abs(y(i)))
end do
sm = 0.0
do i = 1,n
sm = amax1(abs(x(i)-y(i))/s1,sm)
x(i) = y(i)
end do
print *,x
end do
end do


end program