      PROGRAM MAIN
      PARAMETER (NDIM=3)
      DOUBLE PRECISION VNORM2
      DOUBLE PRECISION A(NDIM,NDIM), Y(NDIM), HILF(0:NDIM), B(NDIM),
     $ AT(NDIM,NDIM), MA(NDIM,NDIM), X(NDIM)
      INTEGER N, M, I, J, INDEX(NDIM)

      N=2
      M=3

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

      A(1,1)=1.0d0/DEXP(1.0D0)
      A(2,1)=1.0d0
      A(3,1)=DEXP(1.0D0)
      A(1,2)=-1.0D0
      A(2,2)=0.0D0
      A(3,2)=1.0D0

****** BERECHNE DIE TRANSPONIERTE VON A
      CALL TRANSP(NDIM,M,N,A,AT)

      Y(1)=-1.0D0
      Y(2)=A(3,1)
      Y(3)=2.0D0+A(3,1)**2

****** BERECHNE AT MAL Y      
      CALL MATVEK(NDIM,N,M,AT,Y,X)
****** BERECHNE AT MAL A
      CALL MATMUL(NDIM,N,M,N,AT,A,MA)

****** L™SE DIE NORMALENGLEICHUNGEN
      CALL GAUSEL(NDIM,N,MA,X,INDEX,HILF)

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

****** AUSGABE                           ************
      WRITE(*,*)
      WRITE(*,*) 'L”sung:'
      WRITE(*,*)
      DO 17, I=1,N
         WRITE(*,500) I,X(I)
 17   CONTINUE
      WRITE(*,*)


****** BERECHNE DAS RESIDUUM
      CALL MATVEK(NDIM,M,N,A,X,B)
      DO 5, J=1,M
         B(J)=B(J)-Y(J)
 5    CONTINUE
      WRITE(*,*) 'RESIDUUM: ', VNORM2(NDIM,M,B)
      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 MATMUL(ND,N,M,L,B,C,A)
      INTEGER ND, N, M, L, I, J, K
      DOUBLE PRECISION B(ND,ND), C(ND,ND), A(ND,ND)

      DO 400, I=1,N
         DO 410, J=1,L
            A(I,J)=0.0D0
            DO 420, K=1,M
               A(I,J)=A(I,J)+B(I,K)*C(K,J)
 420        CONTINUE
 410     CONTINUE
 400  CONTINUE
      END


      SUBROUTINE MATVEK(ND,N,M,B,C,X)
      INTEGER ND, N, M, I, J
      DOUBLE PRECISION B(ND,ND), C(ND), X(ND)

      DO 450, I=1,N
         X(I)=0.0D0
         DO 460, J=1,M
            X(I)=X(I)+B(I,J)*C(J)
 460     CONTINUE
 450  CONTINUE
      END


      SUBROUTINE TRANSP(ND,N,M,A,AT)
      INTEGER ND, N, M, I, J
      DOUBLE PRECISION A(ND,ND), AT(ND,ND)

      DO 480, I=1,N
         DO 490, J=1,M
            AT(J,I)=A(I,J) 
 490     CONTINUE
 480  CONTINUE
      END


      DOUBLE PRECISION FUNCTION VNORM2(ND,N,X)
      INTEGER ND, N, I
      DOUBLE PRECISION X(ND), SUM
      
      SUM=0.0D0
      DO 600, I=1,N
         SUM=SUM+X(I)**2
 600  CONTINUE
      VNORM2=DSQRT(SUM)
      END
