gauss.f

Go to the documentation of this file.

      SUBROUTINE lu(A,N,ldim,IR,IC)
C  IT IS THE FIRST SUBROUTINE TO COMPUTE THE MX. INV.
      dimension a(ldim,1),ir(1),ic(1)
      DO 10i=1,n
      ir(i)=i
      ic(i)=i
10      CONTINUE
      k=1
      l=k
      m=k
      xmax=abs(a(k,k))
      DO 100i=k,n
      DO 100j=k,n
      y=abs(a(i,j))
      IF(xmax.GE.y) GOTO 100
      xmax=y
      l=i
      m=j
100     CONTINUE
      DO 1000k=1,n
      irl=ir(l)
      ir(l)=ir(k)
      ir(k)=irl
      icm=ic(m)
      ic(m)=ic(k)
      ic(k)=icm
      IF(l.EQ.k) GOTO 300
      DO 200j=1,n
      b=a(k,j)
      a(k,j)=a(l,j)
      a(l,j)=b
200     CONTINUE
300     IF(m.EQ.k) GOTO 500
      DO 400i=1,n
      b=a(i,k)
      a(i,k)=a(i,m)
       a(i,m)=b
400    CONTINUE
500     c=1./a(k,k)
      a(k,k)=c
      IF(k.EQ.n) GOTO 1000
      k1=k+1
      xmax=abs(a(k1,k1))
      l=k1
      m=k1
      DO 600i=k1,n
       a(i,k)=c*a(i,k)
600     CONTINUE
      DO 800i=k1,n
      b=a(i,k)
      DO 800j=k1,n
      a(i,j)=a(i,j)-b*a(k,j)
      y=abs(a(i,j))
      IF(xmax.GE.y) GOTO 800
      xmax=y
      l=i
      m=j
800    CONTINUE
1000  CONTINUE
      RETURN
      END
      SUBROUTINE solve(F,A,K,N,ldim,IR,IC)
C   IT IS THE SECOND PART OF THE MATRIX INVERSION
      dimension a(ldim,1),f(ldim,1),ir(1),ic(1)
      COMMON /ctmpg/ g(2000)
C
C
      IF (n.GT.2000) THEN
         write(6,*) 'Abort IN Subrtouine SOLVE: N>2000, N=',n
         call exitt
      ENDIF
C
      n1=n+1
      DO 1000kk=1,k
      DO 100i=1,n
      iri=ir(i)
        g(i)=f(iri,kk)
100     CONTINUE
      DO 400i=2,n
      i1=i-1
      b=g(i)
      DO 300j=1,i1
        b=b-a(i,j)*g(j)
300     CONTINUE
        g(i)=b
400     CONTINUE
      DO 700it=1,n
      i=n1-it
      i1=i+1
      b=g(i)
      IF(i.EQ.n) GOTO 701
      DO 600j=i1,n
        b=b-a(i,j)*g(j)
600     CONTINUE
701     g(i)=b*a(i,i)
700     CONTINUE
      DO 900i=1,n
      ici=ic(i)
        f(ici,kk)=g(i)
900     CONTINUE
1000    CONTINUE
      RETURN
      END