      SUBROUTINE LUDCMP(A,N,INDX,D,KHAD)
!
!  Dieses Programm fuehrt eine LU-Aufspaltung einer reellen quadratischen
!  Matrix durch; gleichzeitig werden die Determinante und die
!  Hadamardsche Konditionszahl berechnet.
!
!  Last update: 2-10-96
!
!  Quelle: NUMREC 1986  S. 38 f  (mit Aenderungen)
!
!  DOUBLE PRECISION VERSION

      implicit none

      integer n,indx(n),i,j,k,imax
      real*8  a(n,n),d,khad,sum,aamax,dum,tiny
      PARAMETER (TINY=1.0D-20)

      khad=1.d0
      do i=1,n
        sum=0.d0
        do j=1,n
          sum=sum+a(i,j)**2
        enddo
        khad=khad*sqrt(sum)
      enddo

      D=1.d0
      DO J=1,N
        DO I=1,J-1
          SUM=A(I,J)
          DO K=1,I-1
            SUM=SUM-A(I,K)*A(K,J)
          enddo
          A(I,J)=SUM
        enddo
        AAMAX=0.d0
        DO I=J,N
          SUM=A(I,J)
          DO K=1,J-1
            SUM=SUM-A(I,K)*A(K,J)
          enddo
          A(I,J)=SUM
          DUM=ABS(SUM)
          IF (DUM.GE.AAMAX) THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
        enddo
        IF (J.NE.IMAX)THEN
          DO K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(J,K)
            A(J,K)=DUM
          enddo
          D=-D
        ENDIF
        INDX(J)=IMAX
        IF(A(J,J).EQ.0.d0)THEN
	  A(J,J)=TINY
	ENDIF
        DUM=1.d0/A(J,J)
        DO I=J+1,N
          A(I,J)=A(I,J)*DUM
        enddo
        d=d*a(j,j)
      enddo

      khad=abs(d)/khad

      RETURN
      END




      SUBROUTINE LUBKSB(A,N,INDX,B,X)
!
!  Dieses Programm berechnet den Loesungsvektor X des Systems
!
!                   LU X = B
!
!  unter Verwendung der Gleichungen (2.10-2.13) im Skriptum.
!
!  Last update: 2-10-96
!
!  Quelle: NUMREC 1986  S. 39  (mit Aenderungen)
!
!  DOUBLE PRECISION VERSION

      implicit none

      integer n,indx(n),i,j,ll
!
      real*8  a(n,n),b(n),x(n),bb(1000),sum,y(1000)
!
      do i=1,n
        bb(i)=b(i)
      enddo
      DO I=1,N
        LL=INDX(I)
        SUM=BB(LL)
        BB(LL)=BB(I)
        DO J=1,I-1
          SUM=SUM-A(I,J)*Y(J)
        enddo
        Y(I)=SUM
      enddo
      DO I=N,1,-1
        SUM=Y(I)
        DO J=I+1,N
          SUM=SUM-A(I,J)*X(J)
        enddo
        X(I)=SUM/A(I,I)
      enddo
      RETURN
      END
