lapack_gen.F

#if defined LINUX || defined APPLE
!----------------------------------------------------------------------
! From Numerical Recipes
! added two Numerical Recipes routines for matrix inversion
! LUBKSB - solves a system of linear equations, follows call to LUDCMP
! LUDCMP - replaces an NxN matrix a with the LU decomposition
! 
!  2012-11-05 E.Mirvis separated this generic LU from the splat.F
!----------------------------------------------------------------------
      SUBROUTINE lubksb(A,N,NP,INDX,B)
      REAL A(NP,NP),B(N)
      INTEGER INDX(N)
      ii=0  
      DO 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        IF (ii.NE.0)THEN
          DO 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
   11     CONTINUE
        ELSE IF (sum.NE.0.) THEN
          ii=i
        ENDIF
        b(i)=sum
   12 CONTINUE
      DO 14 i=n,1,-1  
        sum=b(i)
        IF(i.LT.n)THEN
          DO 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
   13 CONTINUE
        ENDIF
        b(i)=sum/a(i,i)
   14 CONTINUE
      RETURN
      END
 
      SUBROUTINE ludcmp(A,N,NP,INDX)
C      PARAMETER (NMAX=400,TINY=1.0E-20)
      parameter(tiny=1.0e-20)
C==EM==^^^
C
      REAL A(NP,NP),VV(N),D
C     REAL A(NP,NP),VV(NMAX),D
C==EM==^^^
      INTEGER INDX(N)
      d=1.  
      DO 12 i=1,n
        aamax=0.
        DO 11 j=1,n
          IF (abs(a(i,j)).GT.aamax) aamax=abs(a(i,j))
   11   CONTINUE
        IF (aamax.EQ.0.) print *, 'SINGULAR MATRIX.'
        vv(i)=1./aamax
   12 CONTINUE
      DO 19 j=1,n
        IF (j.GT.1) THEN
          DO 14 i=1,j-1
            sum=a(i,j)
            IF (i.GT.1)THEN
              DO 13 k=1,i-1
                sum=sum-a(i,k)*a(k,j)
   13         CONTINUE
              a(i,j)=sum
            ENDIF
   14     CONTINUE
        ENDIF
        aamax=0.
        DO 16 i=j,n
          sum=a(i,j)
          IF (j.GT.1)THEN
            DO 15 k=1,j-1
              sum=sum-a(i,k)*a(k,j)
   15       CONTINUE  
            a(i,j)=sum
          ENDIF
          dum=vv(i)*abs(sum)
          IF (dum.GE.aamax) THEN
            imax=i
            aamax=dum
          ENDIF
   16   CONTINUE
        IF (j.NE.imax)THEN
          DO 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
   17     CONTINUE
          d=-d
          vv(imax)=vv(j)
        ENDIF
        indx(j)=imax
        IF(j.NE.n)THEN
          IF(a(j,j).EQ.0.)a(j,j)=tiny
          dum=1./a(j,j)
          DO 18 i=j+1,n
            a(i,j)=a(i,j)*dum
   18     CONTINUE
        ENDIF
   19 CONTINUE
      IF(a(n,n).EQ.0.)a(n,n)=tiny
      RETURN
      END
#endif