0

I am somewhat new to programming in fortran and I have been doing small activities, one of them is a subroutine to solve a system of equations with the Gauss Jordan Method:

     PROGRAM Program2
     REAL(8), DIMENSION(:,:), ALLOCATABLE ::  a,b
     INTEGER                              ::  i, j, p, n
     REAL(8), DIMENSION(:), ALLOCATABLE   ::  temp,x
     REAL(8)                              ::   suma,m
     PRINT "(A)", "Enter number of equations"
     READ *, n

     ALLOCATE (a(n,n+1), temp(n+1),x(n),b(n,n+1))
     PRINT*, ' Write the rows (are of size ',n+1,')'
     DO i = 1, n
     READ (*,*) (a(i,j), j=1,N+1)
     END DO
     b=a   
     DO i = 1, n
     PRINT *, (b(i,j), j=1,n+1)
     
     END DO

     call Gauss_Jordan(n,a,b,temp,x)
    
     END PROGRAM Program2
      !***************** Subroutine that performs Gauss-Jordan elimination*****************
     subroutine Gauss_Jordan (n,a,b,temp,x)  
     implicit none
     INTEGER                              ::  i, j, p, n
     REAL(8), DIMENSION(n,n+1)            ::  a,b 
     REAL(8), DIMENSION(n+1)              ::  temp 
     REAL(8), DIMENSION(n)                ::  x
      REAL(8)                              ::  sum,m  
      DO i = 1, n-1                                        
      P=0   
      DO j = i, n                                       
      IF (a(j,i) /= .0 ) THEN
     p=j
        EXIT
       ENDIF
    ENDDO
     IF ( P .eq. 0) print*, 'The solution is not unique'
     IF ( P .eq. 0) stop                             
     IF ( P  /=  i) THEN                                  
     ! Swap row p with row i.
     temp = a(p,:);  a(p,:) = a(i,:);  a(i,:) = temp
     ENDIF
     DO j = i+1, n                                    
     m = a(j,i)/a(i,i)                             
     a(j,:) = a(j,:) - m * a(i,:)                 
     ENDDO
     ENDDO
     IF ( a(n,n) .eq. 0) THEN
     PRINT*,' There is no unique solution '
     STOP     
     ENDIF
     !******************Start substitution backwards******************** 
     x(n)=a(n,n+1)/a(n,n)
     DO i = n-1,1,-1
     sum=.0
     DO j = i+1,n
     sum = sum + a(i,j)*x(j)
     ENDDO
     x(i)=(a(i,n+1)-sum)/a(i,i)
     ENDDO
     Print*, x
     print*, 'residuo =',b(:,n+1)- MATMUL(b(:,1:n),x)
    !print*,B(:,N+1),MATMUL(B(:,1:N),X)
    end subroutine Gauss_Jordan

My question is if it is possible to find the inverse matrix, using this same subroutine, I think the system would be Ax=b where A is the n x n square matrix to find the inverse and b will be the vectors (1,0...,0) , (0,1,0,...,0), (0, 0,1,...,0) ...(0,...,1) (transposed) of size n, but I'm really lost on how to reflect this in the program

Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
Ixe Sánchez Ibarra
  • 3
  • 1
  • 3
    As you are learning Fortran can I point out you are picking up some bad habits - You should 1) Always, *Always*, **Always** use ` Implicit None` - I won't help people who don't use it in their questions 2) `Real(8)` is bad practice and not portable - see https://stackoverflow.com/questions/838310/fortran-90-kind-parameter 3) External subroutines as you have used here are not the best way, try to always use contained or, best, module procedures – Ian Bush May 25 '22 at 04:42
  • 2
    Also, use some consistent indentation. Whitespace at the start of lines that help visualize the structure of the program. Right now it is really hard to see which END DO belongs to which DO and which lines are inside some constructs. – Vladimir F Героям слава May 25 '22 at 05:49

0 Answers0