SUBROUTINE GNLRD(A,N,NDIM1,IDX)
*     LR decomposition of general N*N quadratic matrix in A(NDIM1,N)
*     with result overwriting matrix A. 
*     IDX(N) is an array, returned with information about the order,   
*     with IDX(1)=0 in case of singular matrix.
*     Limitation to N <= NDIMP! 
      INTEGER N,IDX(N),NDIM1,NDIMP,I,J,K,IMAX 
      PARAMETER (NDIMP=100) 
      REAL A(NDIM1,N),SCALE(NDIMP),ELMAX,TMP,SUM
      IF(N.GT.NDIMLD) STOP ' GNLRD: argument N too large '
      DO I=1,N
       ELMAX=0.0
       DO J=1,N 
        IF(ABS(A(I,J)).GT.ELMAX) ELMAX=ABS(A(I,I))
       END DO 
       IF(ELMAX.EQ.0.0) THEN
          IDX(1)=0 ! matrix singular
          RETURN   
       END IF 
       SCALE(I)=1.0/ELMAX 
      END DO 
      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)
        END DO 
        A(I,J)=SUM 
       END DO 
       ELMAX=0.0
       DO I=J,N
        SUM=A(I,J)
        DO K=1,J-1
         SUM=SUM-A(I,K)*A(K,J)
        END DO  
        A(I,J)=SUM
        TMP=SCALE(I)*ABS(SUM)
        IF(TMP.GE.ELMAX) THEN 
           IMAX=I
           ELMAX=TMP 
        END IF 
       END DO  
       IF(J.NE.IMAX) THEN 
          DO K=1,N
           TMP=A(IMAX,K)
           A(IMAX,K)=A(J,K)
           A(J,K)=TMP 
          END DO  
          SCALE(IMAX)=SCALE(J)  
       END IF
       IDX(J)=IMAX
       IF(J.NE.N) THEN 
          TMP=1.0/A(J,J)
          DO I=J+1,N
           A(I,J)=A(I,J)*TMP 
          END DO  
       END IF 
      END DO  
      END