diff --git a/ext/math/Makefile.in b/ext/math/Makefile.in
index a09febbcc..84197212e 100755
--- a/ext/math/Makefile.in
+++ b/ext/math/Makefile.in
@@ -15,7 +15,7 @@ LIB = @buildlib@/libctmath.a
 all: $(LIB)
 
 
-F_FLAGS = @FFLAGS@ @F77FLAGS@
+F_FLAGS = @FFLAGS@ 
 
 OBJS = \
 mach.o \
diff --git a/ext/math/dgbfa.f b/ext/math/dgbfa.f
new file mode 100644
index 000000000..c26e6f579
--- /dev/null
+++ b/ext/math/dgbfa.f
@@ -0,0 +1,174 @@
+      subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
+      integer lda,n,ml,mu,ipvt(1),info
+      double precision abd(lda,1)
+c
+c     dgbfa factors a double precision band matrix by elimination.
+c
+c     dgbfa is usually called by dgbco, but it can be called
+c     directly with a saving in time if  rcond  is not needed.
+c
+c     on entry
+c
+c        abd     double precision(lda, n)
+c                contains the matrix in band storage.  the columns
+c                of the matrix are stored in the columns of  abd  and
+c                the diagonals of the matrix are stored in rows
+c                ml+1 through 2*ml+mu+1 of  abd .
+c                see the comments below for details.
+c
+c        lda     integer
+c                the leading dimension of the array  abd .
+c                lda must be .ge. 2*ml + mu + 1 .
+c
+c        n       integer
+c                the order of the original matrix.
+c
+c        ml      integer
+c                number of diagonals below the main diagonal.
+c                0 .le. ml .lt. n .
+c
+c        mu      integer
+c                number of diagonals above the main diagonal.
+c                0 .le. mu .lt. n .
+c                more efficient if  ml .le. mu .
+c     on return
+c
+c        abd     an upper triangular matrix in band storage and
+c                the multipliers which were used to obtain it.
+c                the factorization can be written  a = l*u  where
+c                l  is a product of permutation and unit lower
+c                triangular matrices and  u  is upper triangular.
+c
+c        ipvt    integer(n)
+c                an integer vector of pivot indices.
+c
+c        info    integer
+c                = 0  normal value.
+c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
+c                     condition for this subroutine, but it does
+c                     indicate that dgbsl will divide by zero if
+c                     called.  use  rcond  in dgbco for a reliable
+c                     indication of singularity.
+c
+c     band storage
+c
+c           if  a  is a band matrix, the following program segment
+c           will set up the input.
+c
+c                   ml = (band width below the diagonal)
+c                   mu = (band width above the diagonal)
+c                   m = ml + mu + 1
+c                   do 20 j = 1, n
+c                      i1 = max0(1, j-mu)
+c                      i2 = min0(n, j+ml)
+c                      do 10 i = i1, i2
+c                         k = i - j + m
+c                         abd(k,j) = a(i,j)
+c                10    continue
+c                20 continue
+c
+c           this uses rows  ml+1  through  2*ml+mu+1  of  abd .
+c           in addition, the first  ml  rows in  abd  are used for
+c           elements generated during the triangularization.
+c           the total number of rows needed in  abd  is  2*ml+mu+1 .
+c           the  ml+mu by ml+mu  upper left triangle and the
+c           ml by ml  lower right triangle are not referenced.
+c
+c     linpack. this version dated 08/14/78 .
+c     cleve moler, university of new mexico, argonne national lab.
+c
+c     subroutines and functions
+c
+c     blas daxpy,dscal,idamax
+c     fortran max0,min0
+c
+c     internal variables
+c
+      double precision t
+      integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
+c
+c
+      m = ml + mu + 1
+      info = 0
+c
+c     zero initial fill-in columns
+c
+      j0 = mu + 2
+      j1 = min0(n,m) - 1
+      if (j1 .lt. j0) go to 30
+      do 20 jz = j0, j1
+         i0 = m + 1 - jz
+         do 10 i = i0, ml
+            abd(i,jz) = 0.0d0
+   10    continue
+   20 continue
+   30 continue
+      jz = j1
+      ju = 0
+c
+c     gaussian elimination with partial pivoting
+c
+      nm1 = n - 1
+      if (nm1 .lt. 1) go to 130
+      do 120 k = 1, nm1
+         kp1 = k + 1
+c
+c        zero next fill-in column
+c
+         jz = jz + 1
+         if (jz .gt. n) go to 50
+         if (ml .lt. 1) go to 50
+            do 40 i = 1, ml
+               abd(i,jz) = 0.0d0
+   40       continue
+   50    continue
+c
+c        find l = pivot index
+c
+         lm = min0(ml,n-k)
+         l = idamax(lm+1,abd(m,k),1) + m - 1
+         ipvt(k) = l + k - m
+c
+c        zero pivot implies this column already triangularized
+c
+         if (abd(l,k) .eq. 0.0d0) go to 100
+c
+c           interchange if necessary
+c
+            if (l .eq. m) go to 60
+               t = abd(l,k)
+               abd(l,k) = abd(m,k)
+               abd(m,k) = t
+   60       continue
+c
+c           compute multipliers
+c
+            t = -1.0d0/abd(m,k)
+            call dscal(lm,t,abd(m+1,k),1)
+c
+c           row elimination with column indexing
+c
+            ju = min0(max0(ju,mu+ipvt(k)),n)
+            mm = m
+            if (ju .lt. kp1) go to 90
+            do 80 j = kp1, ju
+               l = l - 1
+               mm = mm - 1
+               t = abd(l,j)
+               if (l .eq. mm) go to 70
+                  abd(l,j) = abd(mm,j)
+                  abd(mm,j) = t
+   70          continue
+               call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
+   80       continue
+   90       continue
+         go to 110
+  100    continue
+            info = k
+  110    continue
+  120 continue
+  130 continue
+      ipvt(n) = n
+      if (abd(m,n) .eq. 0.0d0) info = n
+      return
+      end