gauss.F90

Go to the documentation of this file.

SUBROUTINE lu(A,N,NDIM,IR,IC)
  use kinds, only : dp
  implicit none
  integer :: ndim, n, ir(*), ic(*)
  real(DP) :: A(ndim,*)

  integer :: i, j, k, l, m, irl, icm, k1
  real(DP) :: xmax, y, b, c

  DO i=1,n
      ir(i)=i
      ic(i)=i
  END DO
  k=1
  l=k
  m=k
  xmax=abs(a(k,k))
  DO i=k,n
      DO j=k,n
          y=abs(a(i,j))
          IF(xmax >= y) cycle
          xmax=y
          l=i
          m=j
      END DO
  enddo 
  DO k=1,n
      irl=ir(l)
      ir(l)=ir(k)
      ir(k)=irl
      icm=ic(m)
      ic(m)=ic(k)
      ic(k)=icm
      IF(l /= k) then 
         DO j=1,n
             b=a(k,j)
             a(k,j)=a(l,j)
             a(l,j)=b
         END DO
      endif
      IF(m /= k) then 
         DO i=1,n
             b=a(i,k)
             a(i,k)=a(i,m)
             a(i,m)=b
         END DO
      endif
      c=1./a(k,k)
      a(k,k)=c
      IF(k == n) cycle
      k1=k+1
      xmax=abs(a(k1,k1))
      l=k1
      m=k1
      DO i=k1,n
          a(i,k)=c*a(i,k)
      END DO
      DO i=k1,n
          b=a(i,k)
          DO j=k1,n
              a(i,j)=a(i,j)-b*a(k,j)
              y=abs(a(i,j))
              IF(xmax >= y) cycle
              xmax=y
              l=i
              m=j
          END DO
       END DO 
  END DO 
  RETURN
END SUBROUTINE lu

SUBROUTINE solve(F,A,K,N,NDIM,IR,IC)
  use kinds, only : dp
  implicit none

  integer :: n, k, ndim
  real(DP) :: A(ndim,*),F(ndim,*)
  integer :: IR(*),IC(*) 

  real(DP) :: G(2000)

  integer :: n1, kk, i, iri, i1, j, it, ici
  real(DP) :: B

  IF (n > 2000) THEN
      write(6,*) 'Abort IN Subrtouine SOLVE: N>2000, N=',n
      call exitt
  ENDIF

  n1=n+1
  DO 1000 kk=1,k
      DO 100 i=1,n
          iri=ir(i)
          g(i)=f(iri,kk)
      100 END DO
      DO 400 i=2,n
          i1=i-1
          b=g(i)
          DO 300 j=1,i1
              b=b-a(i,j)*g(j)
          300 END DO
          g(i)=b
      400 END DO
      DO 700 it=1,n
          i=n1-it
          i1=i+1
          b=g(i)
          IF(i == n) goto 701
          DO 600 j=i1,n
              b=b-a(i,j)*g(j)
          600 END DO
          701 g(i)=b*a(i,i)
      700 END DO
      DO 900 i=1,n
          ici=ic(i)
          f(ici,kk)=g(i)
      900 END DO
  1000 END DO
  RETURN
END SUBROUTINE solve