0

As per title, what's the best algorithm to numerically solve a linear system in Fortran, if this system has 0s along the main diagonal?

Up to now, I had been fine using simple Gaussian elimination:

SUBROUTINE solve_lin_sys(A, c, x, n)
! =====================================================
! Uses gauss elimination and backwards substitution
! to reduce a linear system and solve it.
! Problem (0-div) can arise if there are 0s on main diag.
! =====================================================
  IMPLICIT NONE
  INTEGER:: i, j, k
  REAL*8::fakt, summ
  INTEGER, INTENT(in):: n
  REAL*8, INTENT(inout):: A(n,n)
  REAL*8, INTENT(inout):: c(n)
  REAL*8, INTENT(out):: x(n)

  DO i = 1, n-1                   ! pick the var to eliminate
    DO j = i+1, n                 ! pick the row where to eliminate
      fakt = A(j,i) / A(i,i)      ! elimination factor
      DO k = 1, n                 ! eliminate
        A(j,k) = A(j,k) - A(i,k)*fakt
      END DO
      c(j)=c(j)-c(i)*fakt         ! iterate on known terms
    END DO
  END DO

   ! Actual solving:
    x(n) = c(n) / A(n,n)          ! last variable being solved
    DO i = n-1, 1, -1
      summ = 0.d0
      DO j = i+1, n
        summ = summ + A(i,j)*x(j)
      END DO
      x(i) = (c(i) - summ) / A(i,i)
    END DO

END SUBROUTINE solve_lin_sys

As you can see I'm dividing by A(i,i) in the calculations. Same problem arises using Gauss-Jordan transformation or Gauss-Seidel elimination.

What's the best solution? I know i'm probably missing some really basic step but I'm a beginner programmer and apparently my linear algebra is getting rusty.

Ian Bush
  • 6,996
  • 1
  • 21
  • 27
Marco Leonardi
  • 149
  • 7
  • 2
    Instead of straight Gaussian elimination, you need to using pivoting where the equations are re-ordered to give a nonzero pivot value. I'll suggest that you should also look into LU decomposition. – steve Aug 31 '22 at 18:06
  • 4
    Use LAPACK. Nobody should be writing linear solvers for dense systems nowadays, and they shouldn't have been for the whole of this millennium. – Ian Bush Aug 31 '22 at 20:32
  • Also note real*8 is not part of standard fortran, has never been part of standard fortran, might not be supported by your next compiler, and might not do what you expect. The standard way is via the kind mechanism - see https://stackoverflow.com/questions/838310/fortran-90-kind-parameter – Ian Bush Sep 01 '22 at 09:02
  • @steve: pivoting and LU decomposition are two orthogonal concepts. What the OP is after is (full or partial) pivoting. –  Sep 01 '22 at 13:18
  • @YvesDaoust, Yeah, I know pivoting and LU decomposition are different concepts. If you read what I wrote. I tell OP he needs pivoting for his Gaussian elimination method (that's the 1st sentence). I then **suggest** (in the 2nd sentence) that he may want to look at LU decomposition. It seemed obvious that LU decomposition is another method to solve a linear system. – steve Sep 01 '22 at 17:54
  • @steve: this is only clear in your mind, I guess. By the way, LU decomposition *is* Gaussian elimination. Just an efficient storage scheme is specified. Gauss-Jordan elimination would be a truly different method. –  Sep 01 '22 at 18:37

0 Answers0