      PROGRAM MAIN
      PARAMETER (NDIM=3)
      DOUBLE PRECISION NORM1, NORMUN
      DOUBLE PRECISION A(NDIM,NDIM), B(NDIM), HILF(0:NDIM),
     $ MA(NDIM,NDIM), MI(NDIM,NDIM), EI(NDIM)
      INTEGER N, I, J, INDEX(NDIM)

      N=3
      
      DO 4, I=1,N
         INDEX(I)=I
 4    CONTINUE

      A(1,1)=2.0d0
      A(1,2)=11.0d0
      A(1,3)=14.0d0
      A(2,1)=4.0d0
      A(2,2)=2.0d0
      A(2,3)=2.0d0
      A(3,1)=2.0d0
      A(3,2)=26.0d0
      A(3,3)=11.0d0

      B(1)=5.0D0
      B(2)=4.0D0
      B(3)=-13.0D0

      DO 1, J=1,N
         DO 2, I=1,N
            MA(I,J)=A(I,J)
 2       CONTINUE
 1    CONTINUE

      CALL GAUSEL(NDIM,N,A,B,INDEX,HILF)

****** IST MATRIX A SINGULŽR, BRECHE AB! ************
      IF(N .EQ. -1) THEN
         WRITE(*,*) 'SINGULŽRE MATRIX !    =>  ABBRUCH '
         GOTO 13
      ENDIF

****** BERECHNUNG DER KONDITION VON A
      CALL MATINV(NDIM,N,A,MI,EI,INDEX)

      WRITE(*,*)
      WRITE(*,*) 'KONDITION VON A NACH DER SPALTENSUMMENNORM: ',
     $            NORM1(NDIM,N,MA)*NORM1(NDIM,N,MI)
      WRITE(*,*)
      WRITE(*,*) 'KONDITION VON A NACH DER ZEILENSUMMENNORM:  ',
     $            NORMUN(NDIM,N,MA)*NORMUN(NDIM,N,MI)

****** AUSGABE                           ************
      WRITE(*,*)
      WRITE(*,*) 'L”sung:'
      WRITE(*,*)
      DO 17, I=1,N
         WRITE(*,500) I,B(I)
 17   CONTINUE
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) 'Indexvektor:'
      WRITE(*,*)
      DO 18, I=1,N
         WRITE(*,*) 'i =',I, '   index(i) =',INDEX(I)
 18   CONTINUE
      WRITE(*,*)
 500  FORMAT(1X,'i = ',I2,3x,'x(i) =',1E24.16E2)
 13   END
      


      SUBROUTINE GAUSEL(ND,N,A,B,INDEX,HILF)
      INTEGER PIVOT
      INTEGER ND, N, I, J, K, P, INDEX(ND), HI
      DOUBLE PRECISION A(ND,ND), B(ND), HILF(0:ND)
      DO 10, I=1,N-1
********* SUCHE PIVOTELEMENT     *****************
         P=PIVOT(I,N,A,ND)

********* IST MATRIX A SINGULŽR, BEREITE ABBRUCH VOR *********
         IF(P .EQ. 0) THEN
            N=-1
            GOTO 30
         ENDIF  

********* ZEILENTAUSCH ERFORDERLICH ?
         IF(P .EQ. I) THEN
            GOTO 12
         ELSE
********* ZEILENTAUSCH !!
            DO 11, K=1,N
               HILF(K)=A(I,K)
               A(I,K)=A(P,K)
               A(P,K)=HILF(K)
 11         CONTINUE
            HILF(0)=B(I)
            B(I)=B(P)
            B(P)=HILF(0)
            HI=INDEX(I)
            INDEX(I)=INDEX(P)
            INDEX(P)=HI
 12      ENDIF

********* GAUSSELIMINATION
         DO 15, J=I+1,N
            DO 20, K=I+1,N
               A(J,K)=A(J,K)-A(J,I)/A(I,I)*A(I,K)
 20         CONTINUE

********************************************************************
********* UMFORMUNG VON B AUSKOMMENTIEREN, FALLS LCB BENUTZT WIRD !!
            B(J)=B(J)-A(J,I)/A(I,I)*B(I)
********************************************************************

            A(J,I)=A(J,I)/A(I,I) 
 15      CONTINUE
 10   CONTINUE 

********* IST MATRIX A SINGULŽR, BEREITE ABBRUCH VOR *********  
      IF(A(N,N) .EQ. 0.0D0) THEN
         N=-1
         GOTO 30
      ENDIF

***********************************************************
********* LCB NUR AUFRUFEN, WENN B NICHT UMGEFORMT WURDE !!     
*      CALL LCB(ND,N,A,B)
***********************************************************

      CALL RXC(ND,N,A,B)

 30   END



****** L™ST Rx=c, WENN R IN A UND c IN b VORLIEGT; L™SUNG x GESPEICHERT IN b   
      SUBROUTINE RXC(ND,N,A,B)
      INTEGER ND, N, I, J
      DOUBLE PRECISION A(ND,ND), B(ND)
 
      DO 200, I=N,1,-1
         DO 210, J=N,I+1,-1
            B(I)=B(I)-A(I,J)*B(J)
 210     CONTINUE
         B(I)=B(I)/A(I,I)
 200  CONTINUE
      END


****** L™ST Lc=b, WENN L IN A VORLIEGT; L™SUNG c GESPEICHERT IN b 
      SUBROUTINE LCB(ND,N,A,B)
      INTEGER ND, N, I, J
      DOUBLE PRECISION A(ND,ND), B(ND)
 
      DO 100, I=1,N
         DO 110, J=1,I-1
            B(I)=B(I)-A(I,J)*B(J)
 110     CONTINUE
 100  CONTINUE
      END


****** SUCHT PIVOTELEMENT
      INTEGER FUNCTION PIVOT(I,N,A,ND)
      INTEGER I, N, ND, J, PIV
      DOUBLE PRECISION A(ND,ND)

      PIV=I
      DO 300, J=I+1,N
         IF (ABS(A(J,I)) .GT. ABS(A(PIV,I))) PIV=J
 300  CONTINUE
****** IST MATRIX A SINGULŽR, BEREITE ABBRUCH VOR
      IF (A(PIV,I) .EQ. 0.0D0) THEN
         PIVOT=0
      ELSE
         PIVOT=PIV
      ENDIF
      END


      SUBROUTINE MATINV(ND,N,MA,MI,EI,INDEX)
      INTEGER ND, N, I, J, INDEX(ND)
      DOUBLE PRECISION MA(ND,ND), MI(ND,ND), EI(ND)

      DO 310, J=1,N
         DO 320, I=1,N
            EI(I)=0.0D0
 320     CONTINUE
         EI(INDEX(INDEX(J)))=1.0D0
         CALL LCB(ND,N,MA,EI)
         CALL RXC(ND,N,MA,EI)
         DO 325, I=1,N
            MI(I,J)=EI(I)
 325     CONTINUE   
 310  CONTINUE  
 330  END


****** SPALTENSUMMENNORM
      DOUBLE PRECISION FUNCTION NORM1(ND,N,MA)
      INTEGER ND, N, I, J
      DOUBLE PRECISION MA(ND,ND), MAXIM, SUMME
 
      MAXIM=0.0D0
      DO 340, J=1,N
         SUMME=0.0D0
         DO 350, I=1,N
            SUMME=SUMME+ABS(MA(I,J))
 350     CONTINUE
         IF (SUMME .GT. MAXIM) MAXIM=SUMME
 340  CONTINUE
      NORM1=MAXIM
      END
     

 
****** ZEILENSUMMENNORM
      DOUBLE PRECISION FUNCTION NORMUN(ND,N,MA)
      INTEGER ND, N, I, J
      DOUBLE PRECISION MA(ND,ND), MAXIM, SUMME

      MAXIM=0.0D0
      DO 360, I=1,N
         SUMME=0.0D0
         DO 370, J=1,N
            SUMME=SUMME+ABS(MA(I,J))
 370     CONTINUE
         IF (SUMME .GT. MAXIM) MAXIM=SUMME
 360  CONTINUE
      NORMUN=MAXIM
      END
      
