new file mode 100644
--- /dev/null
+++ b/libcruft/balgen/balgen.f
@@ -0,0 +1,102 @@
+      subroutine balgen (n,ma,a,mb,b,low,igh,cscale,cperm,wk)
+c
+c     *****parameters:
+      integer igh,low,ma,mb,n
+      double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6)
+c
+c     *****local variables:
+c     none
+c
+c     *****functions:
+c     none
+c
+c     *****subroutines called:
+c     reduce, scaleg, gradeq
+c
+c     ---------------------------------------------------------------
+c
+c     *****purpose:
+c     this subroutine balances the matrices a and b to improve the
+c     accuracy of computing the eigensystem of the generalized
+c     eigenproblem a*x = (lambda)*b*x.  the algorithm is specifically
+c     designed to precede qz type algorithms, but improved performance
+c     is expected from most eigensystem solvers.
+c     ref.:  ward, r. c., balancing the generalized eigenvalue
+c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
+c     141-152.
+c
+c     *****parameter description:
+c
+c     on input:
+c
+c       ma,mb   integer
+c               row dimensions of the arrays containing matrices
+c               a and b respectively, as declared in the main calling
+c               program dimension statement;
+c
+c       n       integer
+c               order of the matrices a and b;
+c
+c       a       real(ma,n)
+c               contains the a matrix of the generalized eigenproblem
+c               defined above;
+c
+c       b       real(mb,n)
+c               contains the b matrix of the generalized eigenproblem
+c               defined above;
+c
+c       wk      real(n,6)
+c               work array that must contain at least 6*n storage
+c               locations.  wk is altered by this subroutine.
+c
+c     on output:
+c
+c       a,b     contain the balanced a and b matrices;
+c
+c       low     integer
+c               beginning -1 of the submatrices of a and b
+c               containing the non-isolated eigenvalues;
+c
+c       igh     integer
+c               ending -1 of the submatrices of a and b
+c               containing the non-isolated eigenvalues.  if
+c               igh = 1 (low = 1 also), the a and b matrices have
+c               been permuted into upper triangular form and have
+c               not been balanced;
+c
+c       cscale  real(n)
+c               contains the exponents of the column scaling factors
+c               in its low through igh locations and the reducing
+c               column permutations in its first low-1 and its
+c               igh+1 through n locations;
+c
+c       cperm   real(n)
+c               contains the column permutations applied in grading
+c               the a and b submatrices in its low through igh
+c               locations;
+c
+c       wk      contains the exponents of the row scaling factors
+c               in its low through igh locations, the reducing row
+c               permutations in its first low-1 and its igh+1
+c               through n locations, and the row permutations
+c               applied in grading the a and b submatrices in its
+c               n+low through n+igh locations.
+c
+c     *****algorithm notes:
+c     none
+c
+c     *****history:
+c     written by r. c. ward.......
+c
+c     ---------------------------------------------------------------
+c
+      call reduce (n,ma,a,mb,b,low,igh,cscale,wk)
+      if (low .eq. igh) go to 10
+      call scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk)
+      call gradeq (n,ma,a,mb,b,low,igh,cperm,wk(1,2))
+   10 continue
+      return
+c
+c     last line of balgen
+c
+      end

new file mode 100644
--- /dev/null
+++ b/libcruft/balgen/gradeq.f
@@ -0,0 +1,177 @@
+      subroutine gradeq (n,ma,a,mb,b,low,igh,cperm,wk)
+c
+c     *****parameters:
+      integer igh,low,ma,mb,n
+      double precision a(ma,n),b(mb,n),cperm(n),wk(n,2)
+c
+c     *****local variables:
+      integer i,ighm1,im,ip1,j,jm,jp1,k
+      double precision cmax,rmax,suma,sumb,temp
+c
+c     *****fortran functions:
+      double precision dabs
+c
+c     *****subroutines called:
+c     none
+c
+c     ---------------------------------------------------------------
+c
+c     *****purpose:
+c     this subroutine grades the submatrices of a and b given by
+c     starting -1 low and ending -1 igh in the generalized
+c     eigenvalue problem a*x = (lambda)*b*x by permuting rows and
+c     columns such that the norm of the i-th row (column) of the
+c     a submatrix divided by the norm of the i-th row (column) of
+c     the b submatrix becomes smaller as i increases.
+c     ref.:  ward, r. c., balancing the generalized eigenvalue
+c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
+c     141-152.
+c
+c     *****parameter description:
+c
+c     on input:
+c
+c       ma,mb   integer
+c               row dimensions of the arrays containing matrices
+c               a and b respectively, as declared in the main calling
+c               program dimension statement;
+c
+c       n       integer
+c               order of the matrices a and b;
+c
+c       a       real(ma,n)
+c               contains the a matrix of the generalized eigenproblem
+c               defined above;
+c
+c       b       real(mb,n)
+c               contains the b matrix of the generalized eigenproblem
+c               defined above;
+c
+c       low     integer
+c               specifies the beginning -1 for the rows and
+c               columns of a and b to be graded;
+c
+c       igh     integer
+c               specifies the ending -1 for the rows and columns
+c               of a and b to be graded;
+c
+c       wk      real(n,2)
+c               work array that must contain at least 2*n locations.
+c               only locations low through igh and n+low through
+c               n+igh are referenced by this subroutine.
+c
+c     on output:
+c
+c       a,b     contain the permuted and graded a and b matrices;
+c
+c       cperm   real(n)
+c               contains in its low through igh locations the
+c               column permutations applied in grading the
+c               submatrices.  the other locations are not referenced
+c               by this subroutine;
+c
+c       wk      contains in its low through igh locations the row
+c               permutations applied in grading the submatrices.
+c
+c     *****algorithm notes:
+c     none.
+c
+c     *****history:
+c     written by r. c. ward.......
+c
+c     ---------------------------------------------------------------
+c
+      if (low .eq. igh) go to 510
+      ighm1 = igh-1
+c
+c     compute column norms of a / those of b
+c
+      do 420 j = low,igh
+         suma = 0.0d0
+         sumb = 0.0d0
+         do 410 i = low,igh
+            suma = suma + dabs(a(i,j))
+            sumb = sumb + dabs(b(i,j))
+  410    continue
+         if (sumb .eq. 0.0d0) go to 415
+         wk(j,2) = suma / sumb
+         go to 420
+  415    continue
+         wk(j,2) = 1.0d38
+  420 continue
+c
+c     permute columns to order them by decreasing quotients
+c
+      do 450 j = low,ighm1
+         cmax = wk(j,2)
+         jm = j
+         jp1 = j+1
+         do 430 k = jp1,igh
+            if (cmax .ge. wk(k,2)) go to 430
+            jm = k
+            cmax = wk(k,2)
+  430    continue
+         cperm(j) = jm
+         if (jm .eq. j) go to 450
+         temp = wk(j,2)
+         wk(j,2) = wk(jm,2)
+         wk(jm,2) = temp
+         do 440 i = 1,igh
+            temp = b(i,j)
+            b(i,j) = b(i,jm)
+            b(i,jm) = temp
+            temp = a(i,j)
+            a(i,j) = a(i,jm)
+            a(i,jm) = temp
+  440    continue
+  450 continue
+      cperm(igh) = igh
+c
+c     compute row norms of a / those of b
+c
+      do 470 i = low,igh
+         suma = 0.0d0
+         sumb = 0.0d0
+         do 460 j = low,igh
+            suma = suma + dabs(a(i,j))
+            sumb = sumb + dabs(b(i,j))
+  460    continue
+         if (sumb .eq. 0.0d0) go to 465
+         wk(i,2) = suma / sumb
+         go to 470
+  465    continue
+         wk(i,2) = 1.0d38
+c
+c     permute rows to order them by decreasing quotients
+c
+  470 continue
+      do 500 i = low,ighm1
+         rmax = wk(i,2)
+         im = i
+         ip1 = i+1
+         do 480 k = ip1,igh
+            if (rmax .ge. wk(k,2)) go to 480
+            im = k
+            rmax = wk(k,2)
+  480    continue
+         wk(i,1) = im
+         if (im .eq. i) go to 500
+         temp = wk(i,2)
+         wk(i,2) = wk(im,2)
+         wk(im,2) = temp
+         do 490 j = low,n
+            temp = b(i,j)
+            b(i,j) = b(im,j)
+            b(im,j) = temp
+            temp = a(i,j)
+            a(i,j) = a(im,j)
+            a(im,j) = temp
+  490    continue
+  500 continue
+      wk(igh,1) = igh
+  510 continue
+      return
+c
+c     last line of gradeq
+c
+      end

new file mode 100644
--- /dev/null
+++ b/libcruft/balgen/reduce.f
@@ -0,0 +1,181 @@
+      subroutine reduce (n,ma,a,mb,b,low,igh,cscale,wk)
+c
+c     *****parameters:
+      integer igh,low,ma,mb,n
+      double precision a(ma,n),b(mb,n),cscale(n),wk(n)
+c
+c     *****local variables:
+      integer i,iflow,ii,ip1,is,j,jp1,k,l,lm1,m
+      double precision f
+c
+c     *****functions:
+c     none
+c
+c     *****subroutines called:
+c     none
+c
+c     ---------------------------------------------------------------
+c
+c     *****purpose:
+c     this subroutine reduces, if possible, the order of the
+c     generalized eigenvalue problem a*x = (lambda)*b*x by permuting
+c     the rows and columns of a and b so that they each have the
+c     form
+c                       u  x  y
+c                       0  c  z
+c                       0  0  r
+c
+c     where u and r are upper triangular and c, x, y, and z are
+c     arbitrary.  thus, the isolated eigenvalues corresponding to
+c     the triangular matrices are obtained by a division, leaving
+c     only eigenvalues corresponding to the center matrices to be
+c     computed.
+c     ref.:  ward, r. c., balancing the generalized eigenvalue
+c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
+c     141-152.
+c
+c     *****parameter description:
+c
+c     on input:
+c
+c       ma,mb   integer
+c               row dimensions of the arrays containing matrices
+c               a and b respectively, as declared in the main calling
+c               program dimension statement;
+c
+c       n       integer
+c               order of the matrices a and b;
+c
+c       a       real(ma,n)
+c               contains the a matrix of the generalized eigenproblem
+c               defined above;
+c
+c       b       real(mb,n)
+c               contains the b matrix of the generalized eigenproblem
+c               defined above.
+c
+c     on output:
+c
+c       a,b     contain the permuted a and b matrices;
+c
+c       low     integer
+c               beginning -1 of the submatrices of a and b
+c               containing the non-isolated eigenvalues;
+c
+c       igh     integer
+c               ending -1 of the submatrices of a and b
+c               containing the non-isolated eigenvalues.  if
+c               igh = 1 (low = 1 also), the permuted a and b
+c               matrices are upper triangular;
+c
+c       cscale  real(n)
+c               contains the required column permutations in its
+c               first low-1 and its igh+1 through n locations;
+c
+c       wk      real(n)
+c               contains the required row permutations in its first
+c               low-1 and its igh+1 through n locations.
+c
+c     *****algorithm notes:
+c     none
+c
+c     *****history:
+c     written by r. c. ward.......
+c
+c     ---------------------------------------------------------------
+c
+      k = 1
+      l = n
+      go to 20
+c
+c     find row with one nonzero in columns 1 through l
+c
+   10 continue
+      l = lm1
+      if (l .ne. 1) go to 20
+      wk(1) = 1
+      cscale(1) = 1
+      go to 200
+   20 continue
+      lm1 = l-1
+      do 70 ii = 1,l
+         i = l+1-ii
+         do 30 j = 1,lm1
+            jp1 = j+1
+            if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 40
+   30    continue
+         j = l
+         go to 60
+   40    continue
+         do 50 j = jp1,l
+            if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 70
+   50    continue
+         j = jp1-1
+   60    continue
+         m = l
+         iflow = 1
+         go to 150
+   70 continue
+      go to 90
+c
+c     find column with one nonzero in rows k through n
+c
+   80 continue
+      k = k+1
+   90 continue
+      do 140 j = k,l
+         do 100 i = k,lm1
+            ip1 = i+1
+            if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 110
+  100    continue
+         i = l
+         go to 130
+  110    continue
+         do 120 i = ip1,l
+            if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 140
+  120    continue
+         i = ip1-1
+  130    continue
+         m = k
+         iflow = 2
+         go to 150
+  140 continue
+      go to 200
+c
+c     permute rows m and i
+c
+  150 continue
+      wk(m) = i
+      if (i .eq. m) go to 170
+      do 160 is = k,n
+         f = a(i,is)
+         a(i,is) = a(m,is)
+         a(m,is) = f
+         f = b(i,is)
+         b(i,is) = b(m,is)
+         b(m,is) = f
+  160 continue
+c
+c     permute columns m and j
+c
+  170 continue
+      cscale(m) = j
+      if (j .eq. m) go to 190
+      do 180 is = 1,l
+         f = a(is,j)
+         a(is,j) = a(is,m)
+         a(is,m) = f
+         f = b(is,j)
+         b(is,j) = b(is,m)
+         b(is,m) = f
+  180 continue
+  190 continue
+      go to (10,80), iflow
+  200 continue
+      low = k
+      igh = l
+      return
+c
+c     last line of reduce
+c
+      end

new file mode 100644
--- /dev/null
+++ b/libcruft/balgen/scaleg.f
@@ -0,0 +1,236 @@
+      subroutine scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk)
+c
+c     *****parameters:
+      integer igh,low,ma,mb,n
+      double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6)
+c
+c     *****local variables:
+      integer i,ir,it,j,jc,kount,nr,nrp2
+      double precision alpha,basl,beta,cmax,coef,coef2,coef5,cor,
+     *                 ew,ewc,fi,fj,gamma,pgamma,sum,t,ta,tb,tc
+c
+c     *****fortran functions:
+      double precision dabs, dlog10, dsign
+c     float
+c
+c     *****subroutines called:
+c     none
+c
+c     ---------------------------------------------------------------
+c
+c     *****purpose:
+c     scales the matrices a and b in the generalized eigenvalue
+c     problem a*x = (lambda)*b*x such that the magnitudes of the
+c     elements of the submatrices of a and b (as specified by low
+c     and igh) are close to unity in the least squares sense.
+c     ref.:  ward, r. c., balancing the generalized eigenvalue
+c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
+c     141-152.
+c
+c     *****parameter description:
+c
+c     on input:
+c
+c       ma,mb   integer
+c               row dimensions of the arrays containing matrices
+c               a and b respectively, as declared in the main calling
+c               program dimension statement;
+c
+c       n       integer
+c               order of the matrices a and b;
+c
+c       a       real(ma,n)
+c               contains the a matrix of the generalized eigenproblem
+c               defined above;
+c
+c       b       real(mb,n)
+c               contains the b matrix of the generalized eigenproblem
+c               defined above;
+c
+c       low     integer
+c               specifies the beginning -1 for the rows and
+c               columns of a and b to be scaled;
+c
+c       igh     integer
+c               specifies the ending -1 for the rows and columns
+c               of a and b to be scaled;
+c
+c       cperm   real(n)
+c               work array.  only locations low through igh are
+c               referenced and altered by this subroutine;
+c
+c       wk      real(n,6)
+c               work array that must contain at least 6*n locations.
+c               only locations low through igh, n+low through n+igh,
+c               ..., 5*n+low through 5*n+igh are referenced and
+c               altered by this subroutine.
+c
+c     on output:
+c
+c       a,b     contain the scaled a and b matrices;
+c
+c       cscale  real(n)
+c               contains in its low through igh locations the integer
+c               exponents of 2 used for the column scaling factors.
+c               the other locations are not referenced;
+c
+c       wk      contains in its low through igh locations the integer
+c               exponents of 2 used for the row scaling factors.
+c
+c     *****algorithm notes:
+c     none.
+c
+c     *****history:
+c     written by r. c. ward.......
+c     modified 8/86 by bobby bodenheimer so that if
+c       sum = 0 (corresponding to the case where the matrix
+c       doesn't need to be scaled) the routine returns.
+c
+c     ---------------------------------------------------------------
+c
+      if (low .eq. igh) go to 410
+      do 210 i = low,igh
+         wk(i,1) = 0.0d0
+         wk(i,2) = 0.0d0
+         wk(i,3) = 0.0d0
+         wk(i,4) = 0.0d0
+         wk(i,5) = 0.0d0
+         wk(i,6) = 0.0d0
+         cscale(i) = 0.0d0
+         cperm(i) = 0.0d0
+  210 continue
+c
+c     compute right side vector in resulting linear equations
+c
+      basl = dlog10(2.0d0)
+      do 240 i = low,igh
+         do 240 j = low,igh
+            tb = b(i,j)
+            ta = a(i,j)
+            if (ta .eq. 0.0d0) go to 220
+            ta = dlog10(dabs(ta)) / basl
+  220       continue
+            if (tb .eq. 0.0d0) go to 230
+            tb = dlog10(dabs(tb)) / basl
+  230       continue
+            wk(i,5) = wk(i,5) - ta - tb
+            wk(j,6) = wk(j,6) - ta - tb
+  240 continue
+      nr = igh-low+1
+      coef = 1.0d0/float(2*nr)
+      coef2 = coef*coef
+      coef5 = 0.5d0*coef2
+      nrp2 = nr+2
+      beta = 0.0d0
+      it = 1
+c
+c     start generalized conjugate gradient iteration
+c
+  250 continue
+      ew = 0.0d0
+      ewc = 0.0d0
+      gamma = 0.0d0
+      do 260 i = low,igh
+         gamma = gamma + wk(i,5)*wk(i,5) + wk(i,6)*wk(i,6)
+         ew = ew + wk(i,5)
+         ewc = ewc + wk(i,6)
+  260 continue
+      gamma = coef*gamma - coef2*(ew**2 + ewc**2)
+     +        - coef5*(ew - ewc)**2
+      if (it .ne. 1) beta = gamma / pgamma
+      t = coef5*(ewc - 3.0d0*ew)
+      tc = coef5*(ew - 3.0d0*ewc)
+      do 270 i = low,igh
+         wk(i,2) = beta*wk(i,2) + coef*wk(i,5) + t
+         cperm(i) = beta*cperm(i) + coef*wk(i,6) + tc
+  270 continue
+c
+c     apply matrix to vector
+c
+      do 300 i = low,igh
+         kount = 0
+         sum = 0.0d0
+         do 290 j = low,igh
+            if (a(i,j) .eq. 0.0d0) go to 280
+            kount = kount+1
+            sum = sum + cperm(j)
+  280       continue
+            if (b(i,j) .eq. 0.0d0) go to 290
+            kount = kount+1
+            sum = sum + cperm(j)
+  290    continue
+         wk(i,3) = float(kount)*wk(i,2) + sum
+  300 continue
+      do 330 j = low,igh
+         kount = 0
+         sum = 0.0d0
+         do 320 i = low,igh
+            if (a(i,j) .eq. 0.0d0) go to 310
+            kount = kount+1
+            sum = sum + wk(i,2)
+  310       continue
+            if (b(i,j) .eq. 0.0d0) go to 320
+            kount = kount+1
+            sum = sum + wk(i,2)
+  320    continue
+         wk(j,4) = float(kount)*cperm(j) + sum
+  330 continue
+      sum = 0.0d0
+      do 340 i = low,igh
+         sum = sum + wk(i,2)*wk(i,3) + cperm(i)*wk(i,4)
+  340 continue
+      if(sum.eq.0.0d0) return
+      alpha = gamma / sum
+c
+c     determine correction to current iterate
+c
+      cmax = 0.0d0
+      do 350 i = low,igh
+         cor = alpha * wk(i,2)
+         if (dabs(cor) .gt. cmax) cmax = dabs(cor)
+         wk(i,1) = wk(i,1) + cor
+         cor = alpha * cperm(i)
+         if (dabs(cor) .gt. cmax) cmax = dabs(cor)
+         cscale(i) = cscale(i) + cor
+  350 continue
+      if (cmax .lt. 0.5d0) go to 370
+      do 360 i = low,igh
+         wk(i,5) = wk(i,5) - alpha*wk(i,3)
+         wk(i,6) = wk(i,6) - alpha*wk(i,4)
+  360 continue
+      pgamma = gamma
+      it = it+1
+      if (it .le. nrp2) go to 250
+c
+c     end generalized conjugate gradient iteration
+c
+  370 continue
+      do 380 i = low,igh
+         ir = wk(i,1) + dsign(0.5d0,wk(i,1))
+         wk(i,1) = ir
+         jc = cscale(i) + dsign(0.5d0,cscale(i))
+         cscale(i) = jc
+  380 continue
+c
+c     scale a and b
+c
+      do 400 i = 1,igh
+         ir = wk(i,1)
+         fi = 2.0d0**ir
+         if (i .lt. low) fi = 1.0d0
+         do 400 j =low,n
+            jc = cscale(j)
+            fj = 2.0d0**jc
+            if (j .le. igh) go to 390
+            if (i .lt. low) go to 400
+            fj = 1.0d0
+  390       continue
+            a(i,j) = a(i,j)*fi*fj
+            b(i,j) = b(i,j)*fi*fj
+  400 continue
+  410 continue
+      return
+c
+c     last line of scaleg
+c
+      end

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dasum.f
@@ -0,0 +1,41 @@
+      DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)
+C
+C     TAKES THE SUM OF THE ABSOLUTE VALUES.
+C     JACK DONGARRA, LINPACK, 3/11/78.
+C
+      DOUBLE PRECISION DX(1),DTEMP
+      INTEGER I,INCX,M,MP1,N,NINCX
+C
+      DASUM = 0.0D0
+      DTEMP = 0.0D0
+      IF(N.LE.0)RETURN
+      IF(INCX.EQ.1)GO TO 20
+C
+C        CODE FOR INCREMENT NOT EQUAL TO 1
+C
+      NINCX = N*INCX
+      DO 10 I = 1,NINCX,INCX
+        DTEMP = DTEMP + DABS(DX(I))
+   10 CONTINUE
+      DASUM = DTEMP
+      RETURN
+C
+C        CODE FOR INCREMENT EQUAL TO 1
+C
+C
+C        CLEAN-UP LOOP
+C
+   20 M = MOD(N,6)
+      IF( M .EQ. 0 ) GO TO 40
+      DO 30 I = 1,M
+        DTEMP = DTEMP + DABS(DX(I))
+   30 CONTINUE
+      IF( N .LT. 6 ) GO TO 60
+   40 MP1 = M + 1
+      DO 50 I = MP1,N,6
+        DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2))
+     *  + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5))
+   50 CONTINUE
+   60 DASUM = DTEMP
+      RETURN
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/daxpy.f
@@ -0,0 +1,47 @@
+      SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
+C
+C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
+C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
+C     JACK DONGARRA, LINPACK, 3/11/78.
+C
+      DOUBLE PRECISION DX(1),DY(1),DA
+      INTEGER I,INCX,INCY,IXIY,M,MP1,N
+C
+      IF(N.LE.0)RETURN
+      IF (DA .EQ. 0.0D0) RETURN
+      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
+C
+C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
+C          NOT EQUAL TO 1
+C
+      IX = 1
+      IY = 1
+      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
+      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
+      DO 10 I = 1,N
+        DY(IY) = DY(IY) + DA*DX(IX)
+        IX = IX + INCX
+        IY = IY + INCY
+   10 CONTINUE
+      RETURN
+C
+C        CODE FOR BOTH INCREMENTS EQUAL TO 1
+C
+C
+C        CLEAN-UP LOOP
+C
+   20 M = MOD(N,4)
+      IF( M .EQ. 0 ) GO TO 40
+      DO 30 I = 1,M
+        DY(I) = DY(I) + DA*DX(I)
+   30 CONTINUE
+      IF( N .LT. 4 ) RETURN
+   40 MP1 = M + 1
+      DO 50 I = MP1,N,4
+        DY(I) = DY(I) + DA*DX(I)
+        DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
+        DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
+        DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
+   50 CONTINUE
+      RETURN
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dcabs1.f
@@ -0,0 +1,8 @@
+      double precision function dcabs1(z)
+      double complex z,zz
+      double precision t(2)
+      equivalence (zz,t(1))
+      zz = z
+      dcabs1 = dabs(t(1)) + dabs(t(2))
+      return
+      end

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dcopy.f
@@ -0,0 +1,49 @@
+      SUBROUTINE  DCOPY(N,DX,INCX,DY,INCY)
+C
+C     COPIES A VECTOR, X, TO A VECTOR, Y.
+C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
+C     JACK DONGARRA, LINPACK, 3/11/78.
+C
+      DOUBLE PRECISION DX(1),DY(1)
+      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
+C
+      IF(N.LE.0)RETURN
+      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
+C
+C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
+C          NOT EQUAL TO 1
+C
+      IX = 1
+      IY = 1
+      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
+      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
+      DO 10 I = 1,N
+        DY(IY) = DX(IX)
+        IX = IX + INCX
+        IY = IY + INCY
+   10 CONTINUE
+      RETURN
+C
+C        CODE FOR BOTH INCREMENTS EQUAL TO 1
+C
+C
+C        CLEAN-UP LOOP
+C
+   20 M = MOD(N,7)
+      IF( M .EQ. 0 ) GO TO 40
+      DO 30 I = 1,M
+        DY(I) = DX(I)
+   30 CONTINUE
+      IF( N .LT. 7 ) RETURN
+   40 MP1 = M + 1
+      DO 50 I = MP1,N,7
+        DY(I) = DX(I)
+        DY(I + 1) = DX(I + 1)
+        DY(I + 2) = DX(I + 2)
+        DY(I + 3) = DX(I + 3)
+        DY(I + 4) = DX(I + 4)
+        DY(I + 5) = DX(I + 5)
+        DY(I + 6) = DX(I + 6)
+   50 CONTINUE
+      RETURN
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/ddot.f
@@ -0,0 +1,48 @@
+      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
+C
+C     FORMS THE DOT PRODUCT OF TWO VECTORS.
+C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
+C     JACK DONGARRA, LINPACK, 3/11/78.
+C
+      DOUBLE PRECISION DX(1),DY(1),DTEMP
+      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
+C
+      DDOT = 0.0D0
+      DTEMP = 0.0D0
+      IF(N.LE.0)RETURN
+      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
+C
+C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
+C          NOT EQUAL TO 1
+C
+      IX = 1
+      IY = 1
+      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
+      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
+      DO 10 I = 1,N
+        DTEMP = DTEMP + DX(IX)*DY(IY)
+        IX = IX + INCX
+        IY = IY + INCY
+   10 CONTINUE
+      DDOT = DTEMP
+      RETURN
+C
+C        CODE FOR BOTH INCREMENTS EQUAL TO 1
+C
+C
+C        CLEAN-UP LOOP
+C
+   20 M = MOD(N,5)
+      IF( M .EQ. 0 ) GO TO 40
+      DO 30 I = 1,M
+        DTEMP = DTEMP + DX(I)*DY(I)
+   30 CONTINUE
+      IF( N .LT. 5 ) GO TO 60
+   40 MP1 = M + 1
+      DO 50 I = MP1,N,5
+        DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
+     *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
+   50 CONTINUE
+   60 DDOT = DTEMP
+      RETURN
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dgemm.f
@@ -0,0 +1,347 @@
+************************************************************************
+*
+*     File of the DOUBLE PRECISION Level-3 BLAS.
+*     ==========================================
+*
+*     SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
+*    $                   BETA, C, LDC )
+*
+*     SUBROUTINE DSYMM ( SIDE,   UPLO,   M, N,    ALPHA, A, LDA, B, LDB,
+*    $                   BETA, C, LDC )
+*
+*     SUBROUTINE DSYRK ( UPLO,   TRANS,     N, K, ALPHA, A, LDA,
+*    $                   BETA, C, LDC )
+*
+*     SUBROUTINE DSYR2K( UPLO,   TRANS,     N, K, ALPHA, A, LDA, B, LDB,
+*    $                   BETA, C, LDC )
+*
+*     SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
+*    $                   B, LDB )
+*
+*     SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
+*    $                   B, LDB )
+*
+*     See:
+*
+*        Dongarra J. J.,   Du Croz J. J.,   Duff I.  and   Hammarling S.
+*        A set of  Level 3  Basic Linear Algebra Subprograms.  Technical
+*        Memorandum No.88 (Revision 1), Mathematics and Computer Science
+*        Division,  Argonne National Laboratory, 9700 South Cass Avenue,
+*        Argonne, Illinois 60439.
+*
+*
+************************************************************************
+*
+      SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
+     $                   BETA, C, LDC )
+*     .. Scalar Arguments ..
+      CHARACTER*1        TRANSA, TRANSB
+      INTEGER            M, N, K, LDA, LDB, LDC
+      DOUBLE PRECISION   ALPHA, BETA
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DGEMM  performs one of the matrix-matrix operations
+*
+*     C := alpha*op( A )*op( B ) + beta*C,
+*
+*  where  op( X ) is one of
+*
+*     op( X ) = X   or   op( X ) = X',
+*
+*  alpha and beta are scalars, and A, B and C are matrices, with op( A )
+*  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
+*
+*  Parameters
+*  ==========
+*
+*  TRANSA - CHARACTER*1.
+*           On entry, TRANSA specifies the form of op( A ) to be used in
+*           the matrix multiplication as follows:
+*
+*              TRANSA = 'N' or 'n',  op( A ) = A.
+*
+*              TRANSA = 'T' or 't',  op( A ) = A'.
+*
+*              TRANSA = 'C' or 'c',  op( A ) = A'.
+*
+*           Unchanged on exit.
+*
+*  TRANSB - CHARACTER*1.
+*           On entry, TRANSB specifies the form of op( B ) to be used in
+*           the matrix multiplication as follows:
+*
+*              TRANSB = 'N' or 'n',  op( B ) = B.
+*
+*              TRANSB = 'T' or 't',  op( B ) = B'.
+*
+*              TRANSB = 'C' or 'c',  op( B ) = B'.
+*
+*           Unchanged on exit.
+*
+*  M      - INTEGER.
+*           On entry,  M  specifies  the number  of rows  of the  matrix
+*           op( A )  and of the  matrix  C.  M  must  be at least  zero.
+*           Unchanged on exit.
+*
+*  N      - INTEGER.
+*           On entry,  N  specifies the number  of columns of the matrix
+*           op( B ) and the number of columns of the matrix C. N must be
+*           at least zero.
+*           Unchanged on exit.
+*
+*  K      - INTEGER.
+*           On entry,  K  specifies  the number of columns of the matrix
+*           op( A ) and the number of rows of the matrix op( B ). K must
+*           be at least  zero.
+*           Unchanged on exit.
+*
+*  ALPHA  - DOUBLE PRECISION.
+*           On entry, ALPHA specifies the scalar alpha.
+*           Unchanged on exit.
+*
+*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
+*           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
+*           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
+*           part of the array  A  must contain the matrix  A,  otherwise
+*           the leading  k by m  part of the array  A  must contain  the
+*           matrix A.
+*           Unchanged on exit.
+*
+*  LDA    - INTEGER.
+*           On entry, LDA specifies the first dimension of A as declared
+*           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
+*           LDA must be at least  max( 1, m ), otherwise  LDA must be at
+*           least  max( 1, k ).
+*           Unchanged on exit.
+*
+*  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
+*           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
+*           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
+*           part of the array  B  must contain the matrix  B,  otherwise
+*           the leading  n by k  part of the array  B  must contain  the
+*           matrix B.
+*           Unchanged on exit.
+*
+*  LDB    - INTEGER.
+*           On entry, LDB specifies the first dimension of B as declared
+*           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
+*           LDB must be at least  max( 1, k ), otherwise  LDB must be at
+*           least  max( 1, n ).
+*           Unchanged on exit.
+*
+*  BETA   - DOUBLE PRECISION.
+*           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
+*           supplied as zero then C need not be set on input.
+*           Unchanged on exit.
+*
+*  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
+*           Before entry, the leading  m by n  part of the array  C must
+*           contain the matrix  C,  except when  beta  is zero, in which
+*           case C need not be set on entry.
+*           On exit, the array  C  is overwritten by the  m by n  matrix
+*           ( alpha*op( A )*op( B ) + beta*C ).
+*
+*  LDC    - INTEGER.
+*           On entry, LDC specifies the first dimension of C as declared
+*           in  the  calling  (sub)  program.   LDC  must  be  at  least
+*           max( 1, m ).
+*           Unchanged on exit.
+*
+*
+*  Level 3 Blas routine.
+*
+*  -- Written on 8-February-1989.
+*     Jack Dongarra, Argonne National Laboratory.
+*     Iain Duff, AERE Harwell.
+*     Jeremy Du Croz, Numerical Algorithms Group Ltd.
+*     Sven Hammarling, Numerical Algorithms Group Ltd.
+*
+*
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     .. Local Scalars ..
+      LOGICAL            NOTA, NOTB
+      INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
+      DOUBLE PRECISION   TEMP
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE         , ZERO
+      PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Executable Statements ..
+*
+*     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
+*     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
+*     and  columns of  A  and the  number of  rows  of  B  respectively.
+*
+      NOTA  = LSAME( TRANSA, 'N' )
+      NOTB  = LSAME( TRANSB, 'N' )
+      IF( NOTA )THEN
+         NROWA = M
+         NCOLA = K
+      ELSE
+         NROWA = K
+         NCOLA = M
+      END IF
+      IF( NOTB )THEN
+         NROWB = K
+      ELSE
+         NROWB = N
+      END IF
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF(      ( .NOT.NOTA                 ).AND.
+     $         ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
+     $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
+         INFO = 1
+      ELSE IF( ( .NOT.NOTB                 ).AND.
+     $         ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
+     $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
+         INFO = 2
+      ELSE IF( M  .LT.0               )THEN
+         INFO = 3
+      ELSE IF( N  .LT.0               )THEN
+         INFO = 4
+      ELSE IF( K  .LT.0               )THEN
+         INFO = 5
+      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
+         INFO = 8
+      ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
+         INFO = 10
+      ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
+         INFO = 13
+      END IF
+      IF( INFO.NE.0 )THEN
+         CALL XERBLA( 'DGEMM ', INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
+     $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
+     $   RETURN
+*
+*     And if  alpha.eq.zero.
+*
+      IF( ALPHA.EQ.ZERO )THEN
+         IF( BETA.EQ.ZERO )THEN
+            DO 20, J = 1, N
+               DO 10, I = 1, M
+                  C( I, J ) = ZERO
+   10          CONTINUE
+   20       CONTINUE
+         ELSE
+            DO 40, J = 1, N
+               DO 30, I = 1, M
+                  C( I, J ) = BETA*C( I, J )
+   30          CONTINUE
+   40       CONTINUE
+         END IF
+         RETURN
+      END IF
+*
+*     Start the operations.
+*
+      IF( NOTB )THEN
+         IF( NOTA )THEN
+*
+*           Form  C := alpha*A*B + beta*C.
+*
+            DO 90, J = 1, N
+               IF( BETA.EQ.ZERO )THEN
+                  DO 50, I = 1, M
+                     C( I, J ) = ZERO
+   50             CONTINUE
+               ELSE IF( BETA.NE.ONE )THEN
+                  DO 60, I = 1, M
+                     C( I, J ) = BETA*C( I, J )
+   60             CONTINUE
+               END IF
+               DO 80, L = 1, K
+                  IF( B( L, J ).NE.ZERO )THEN
+                     TEMP = ALPHA*B( L, J )
+                     DO 70, I = 1, M
+                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
+   70                CONTINUE
+                  END IF
+   80          CONTINUE
+   90       CONTINUE
+         ELSE
+*
+*           Form  C := alpha*A'*B + beta*C
+*
+            DO 120, J = 1, N
+               DO 110, I = 1, M
+                  TEMP = ZERO
+                  DO 100, L = 1, K
+                     TEMP = TEMP + A( L, I )*B( L, J )
+  100             CONTINUE
+                  IF( BETA.EQ.ZERO )THEN
+                     C( I, J ) = ALPHA*TEMP
+                  ELSE
+                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
+                  END IF
+  110          CONTINUE
+  120       CONTINUE
+         END IF
+      ELSE
+         IF( NOTA )THEN
+*
+*           Form  C := alpha*A*B' + beta*C
+*
+            DO 170, J = 1, N
+               IF( BETA.EQ.ZERO )THEN
+                  DO 130, I = 1, M
+                     C( I, J ) = ZERO
+  130             CONTINUE
+               ELSE IF( BETA.NE.ONE )THEN
+                  DO 140, I = 1, M
+                     C( I, J ) = BETA*C( I, J )
+  140             CONTINUE
+               END IF
+               DO 160, L = 1, K
+                  IF( B( J, L ).NE.ZERO )THEN
+                     TEMP = ALPHA*B( J, L )
+                     DO 150, I = 1, M
+                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
+  150                CONTINUE
+                  END IF
+  160          CONTINUE
+  170       CONTINUE
+         ELSE
+*
+*           Form  C := alpha*A'*B' + beta*C
+*
+            DO 200, J = 1, N
+               DO 190, I = 1, M
+                  TEMP = ZERO
+                  DO 180, L = 1, K
+                     TEMP = TEMP + A( L, I )*B( J, L )
+  180             CONTINUE
+                  IF( BETA.EQ.ZERO )THEN
+                     C( I, J ) = ALPHA*TEMP
+                  ELSE
+                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
+                  END IF
+  190          CONTINUE
+  200       CONTINUE
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of DGEMM .
+*
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dgemv.f
@@ -0,0 +1,321 @@
+*
+************************************************************************
+*
+*     File of the DOUBLE PRECISION  Level-2 BLAS.
+*     ===========================================
+*
+*     SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
+*    $                   BETA, Y, INCY )
+*
+*     SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
+*    $                   BETA, Y, INCY )
+*
+*     SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX,
+*    $                   BETA, Y, INCY )
+*
+*     SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
+*    $                   BETA, Y, INCY )
+*
+*     SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
+*
+*     SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
+*
+*     SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
+*
+*     SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
+*
+*     SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
+*
+*     SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
+*
+*     SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
+*
+*     SUBROUTINE DGER  ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
+*
+*     SUBROUTINE DSYR  ( UPLO, N, ALPHA, X, INCX, A, LDA )
+*
+*     SUBROUTINE DSPR  ( UPLO, N, ALPHA, X, INCX, AP )
+*
+*     SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
+*
+*     SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
+*
+*     See:
+*
+*        Dongarra J. J., Du Croz J. J., Hammarling S.  and Hanson R. J..
+*        An  extended  set of Fortran  Basic Linear Algebra Subprograms.
+*
+*        Technical  Memoranda  Nos. 41 (revision 3) and 81,  Mathematics
+*        and  Computer Science  Division,  Argonne  National Laboratory,
+*        9700 South Cass Avenue, Argonne, Illinois 60439, US.
+*
+*        Or
+*
+*        NAG  Technical Reports TR3/87 and TR4/87,  Numerical Algorithms
+*        Group  Ltd.,  NAG  Central  Office,  256  Banbury  Road, Oxford
+*        OX2 7DE, UK,  and  Numerical Algorithms Group Inc.,  1101  31st
+*        Street,  Suite 100,  Downers Grove,  Illinois 60515-1263,  USA.
+*
+************************************************************************
+*
+      SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
+     $                   BETA, Y, INCY )
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   ALPHA, BETA
+      INTEGER            INCX, INCY, LDA, M, N
+      CHARACTER*1        TRANS
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DGEMV  performs one of the matrix-vector operations
+*
+*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
+*
+*  where alpha and beta are scalars, x and y are vectors and A is an
+*  m by n matrix.
+*
+*  Parameters
+*  ==========
+*
+*  TRANS  - CHARACTER*1.
+*           On entry, TRANS specifies the operation to be performed as
+*           follows:
+*
+*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
+*
+*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
+*
+*              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
+*
+*           Unchanged on exit.
+*
+*  M      - INTEGER.
+*           On entry, M specifies the number of rows of the matrix A.
+*           M must be at least zero.
+*           Unchanged on exit.
+*
+*  N      - INTEGER.
+*           On entry, N specifies the number of columns of the matrix A.
+*           N must be at least zero.
+*           Unchanged on exit.
+*
+*  ALPHA  - DOUBLE PRECISION.
+*           On entry, ALPHA specifies the scalar alpha.
+*           Unchanged on exit.
+*
+*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
+*           Before entry, the leading m by n part of the array A must
+*           contain the matrix of coefficients.
+*           Unchanged on exit.
+*
+*  LDA    - INTEGER.
+*           On entry, LDA specifies the first dimension of A as declared
+*           in the calling (sub) program. LDA must be at least
+*           max( 1, m ).
+*           Unchanged on exit.
+*
+*  X      - DOUBLE PRECISION array of DIMENSION at least
+*           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
+*           and at least
+*           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
+*           Before entry, the incremented array X must contain the
+*           vector x.
+*           Unchanged on exit.
+*
+*  INCX   - INTEGER.
+*           On entry, INCX specifies the increment for the elements of
+*           X. INCX must not be zero.
+*           Unchanged on exit.
+*
+*  BETA   - DOUBLE PRECISION.
+*           On entry, BETA specifies the scalar beta. When BETA is
+*           supplied as zero then Y need not be set on input.
+*           Unchanged on exit.
+*
+*  Y      - DOUBLE PRECISION array of DIMENSION at least
+*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
+*           and at least
+*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
+*           Before entry with BETA non-zero, the incremented array Y
+*           must contain the vector y. On exit, Y is overwritten by the
+*           updated vector y.
+*
+*  INCY   - INTEGER.
+*           On entry, INCY specifies the increment for the elements of
+*           Y. INCY must not be zero.
+*           Unchanged on exit.
+*
+*
+*  Level 2 Blas routine.
+*
+*  -- Written on 22-October-1986.
+*     Jack Dongarra, Argonne National Lab.
+*     Jeremy Du Croz, Nag Central Office.
+*     Sven Hammarling, Nag Central Office.
+*     Richard Hanson, Sandia National Labs.
+*
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE         , ZERO
+      PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     .. Local Scalars ..
+      DOUBLE PRECISION   TEMP
+      INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
+     $         .NOT.LSAME( TRANS, 'T' ).AND.
+     $         .NOT.LSAME( TRANS, 'C' )      )THEN
+         INFO = 1
+      ELSE IF( M.LT.0 )THEN
+         INFO = 2
+      ELSE IF( N.LT.0 )THEN
+         INFO = 3
+      ELSE IF( LDA.LT.MAX( 1, M ) )THEN
+         INFO = 6
+      ELSE IF( INCX.EQ.0 )THEN
+         INFO = 8
+      ELSE IF( INCY.EQ.0 )THEN
+         INFO = 11
+      END IF
+      IF( INFO.NE.0 )THEN
+         CALL XERBLA( 'DGEMV ', INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
+     $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
+     $   RETURN
+*
+*     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
+*     up the start points in  X  and  Y.
+*
+      IF( LSAME( TRANS, 'N' ) )THEN
+         LENX = N
+         LENY = M
+      ELSE
+         LENX = M
+         LENY = N
+      END IF
+      IF( INCX.GT.0 )THEN
+         KX = 1
+      ELSE
+         KX = 1 - ( LENX - 1 )*INCX
+      END IF
+      IF( INCY.GT.0 )THEN
+         KY = 1
+      ELSE
+         KY = 1 - ( LENY - 1 )*INCY
+      END IF
+*
+*     Start the operations. In this version the elements of A are
+*     accessed sequentially with one pass through A.
+*
+*     First form  y := beta*y.
+*
+      IF( BETA.NE.ONE )THEN
+         IF( INCY.EQ.1 )THEN
+            IF( BETA.EQ.ZERO )THEN
+               DO 10, I = 1, LENY
+                  Y( I ) = ZERO
+   10          CONTINUE
+            ELSE
+               DO 20, I = 1, LENY
+                  Y( I ) = BETA*Y( I )
+   20          CONTINUE
+            END IF
+         ELSE
+            IY = KY
+            IF( BETA.EQ.ZERO )THEN
+               DO 30, I = 1, LENY
+                  Y( IY ) = ZERO
+                  IY      = IY   + INCY
+   30          CONTINUE
+            ELSE
+               DO 40, I = 1, LENY
+                  Y( IY ) = BETA*Y( IY )
+                  IY      = IY           + INCY
+   40          CONTINUE
+            END IF
+         END IF
+      END IF
+      IF( ALPHA.EQ.ZERO )
+     $   RETURN
+      IF( LSAME( TRANS, 'N' ) )THEN
+*
+*        Form  y := alpha*A*x + y.
+*
+         JX = KX
+         IF( INCY.EQ.1 )THEN
+            DO 60, J = 1, N
+               IF( X( JX ).NE.ZERO )THEN
+                  TEMP = ALPHA*X( JX )
+                  DO 50, I = 1, M
+                     Y( I ) = Y( I ) + TEMP*A( I, J )
+   50             CONTINUE
+               END IF
+               JX = JX + INCX
+   60       CONTINUE
+         ELSE
+            DO 80, J = 1, N
+               IF( X( JX ).NE.ZERO )THEN
+                  TEMP = ALPHA*X( JX )
+                  IY   = KY
+                  DO 70, I = 1, M
+                     Y( IY ) = Y( IY ) + TEMP*A( I, J )
+                     IY      = IY      + INCY
+   70             CONTINUE
+               END IF
+               JX = JX + INCX
+   80       CONTINUE
+         END IF
+      ELSE
+*
+*        Form  y := alpha*A'*x + y.
+*
+         JY = KY
+         IF( INCX.EQ.1 )THEN
+            DO 100, J = 1, N
+               TEMP = ZERO
+               DO 90, I = 1, M
+                  TEMP = TEMP + A( I, J )*X( I )
+   90          CONTINUE
+               Y( JY ) = Y( JY ) + ALPHA*TEMP
+               JY      = JY      + INCY
+  100       CONTINUE
+         ELSE
+            DO 120, J = 1, N
+               TEMP = ZERO
+               IX   = KX
+               DO 110, I = 1, M
+                  TEMP = TEMP + A( I, J )*X( IX )
+                  IX   = IX   + INCX
+  110          CONTINUE
+               Y( JY ) = Y( JY ) + ALPHA*TEMP
+               JY      = JY      + INCY
+  120       CONTINUE
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of DGEMV .
+*
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dger.f
@@ -0,0 +1,160 @@
+*
+************************************************************************
+*
+      SUBROUTINE DGER  ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   ALPHA
+      INTEGER            INCX, INCY, LDA, M, N
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DGER   performs the rank 1 operation
+*
+*     A := alpha*x*y' + A,
+*
+*  where alpha is a scalar, x is an m element vector, y is an n element
+*  vector and A is an m by n matrix.
+*
+*  Parameters
+*  ==========
+*
+*  M      - INTEGER.
+*           On entry, M specifies the number of rows of the matrix A.
+*           M must be at least zero.
+*           Unchanged on exit.
+*
+*  N      - INTEGER.
+*           On entry, N specifies the number of columns of the matrix A.
+*           N must be at least zero.
+*           Unchanged on exit.
+*
+*  ALPHA  - DOUBLE PRECISION.
+*           On entry, ALPHA specifies the scalar alpha.
+*           Unchanged on exit.
+*
+*  X      - DOUBLE PRECISION array of dimension at least
+*           ( 1 + ( m - 1 )*abs( INCX ) ).
+*           Before entry, the incremented array X must contain the m
+*           element vector x.
+*           Unchanged on exit.
+*
+*  INCX   - INTEGER.
+*           On entry, INCX specifies the increment for the elements of
+*           X. INCX must not be zero.
+*           Unchanged on exit.
+*
+*  Y      - DOUBLE PRECISION array of dimension at least
+*           ( 1 + ( n - 1 )*abs( INCY ) ).
+*           Before entry, the incremented array Y must contain the n
+*           element vector y.
+*           Unchanged on exit.
+*
+*  INCY   - INTEGER.
+*           On entry, INCY specifies the increment for the elements of
+*           Y. INCY must not be zero.
+*           Unchanged on exit.
+*
+*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
+*           Before entry, the leading m by n part of the array A must
+*           contain the matrix of coefficients. On exit, A is
+*           overwritten by the updated matrix.
+*
+*  LDA    - INTEGER.
+*           On entry, LDA specifies the first dimension of A as declared
+*           in the calling (sub) program. LDA must be at least
+*           max( 1, m ).
+*           Unchanged on exit.
+*
+*
+*  Level 2 Blas routine.
+*
+*  -- Written on 22-October-1986.
+*     Jack Dongarra, Argonne National Lab.
+*     Jeremy Du Croz, Nag Central Office.
+*     Sven Hammarling, Nag Central Office.
+*     Richard Hanson, Sandia National Labs.
+*
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO
+      PARAMETER        ( ZERO = 0.0D+0 )
+*     .. Local Scalars ..
+      DOUBLE PRECISION   TEMP
+      INTEGER            I, INFO, IX, J, JY, KX
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF     ( M.LT.0 )THEN
+         INFO = 1
+      ELSE IF( N.LT.0 )THEN
+         INFO = 2
+      ELSE IF( INCX.EQ.0 )THEN
+         INFO = 5
+      ELSE IF( INCY.EQ.0 )THEN
+         INFO = 7
+      ELSE IF( LDA.LT.MAX( 1, M ) )THEN
+         INFO = 9
+      END IF
+      IF( INFO.NE.0 )THEN
+         CALL XERBLA( 'DGER  ', INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
+     $   RETURN
+*
+*     Start the operations. In this version the elements of A are
+*     accessed sequentially with one pass through A.
+*
+      IF( INCY.GT.0 )THEN
+         JY = 1
+      ELSE
+         JY = 1 - ( N - 1 )*INCY
+      END IF
+      IF( INCX.EQ.1 )THEN
+         DO 20, J = 1, N
+            IF( Y( JY ).NE.ZERO )THEN
+               TEMP = ALPHA*Y( JY )
+               DO 10, I = 1, M
+                  A( I, J ) = A( I, J ) + X( I )*TEMP
+   10          CONTINUE
+            END IF
+            JY = JY + INCY
+   20    CONTINUE
+      ELSE
+         IF( INCX.GT.0 )THEN
+            KX = 1
+         ELSE
+            KX = 1 - ( M - 1 )*INCX
+         END IF
+         DO 40, J = 1, N
+            IF( Y( JY ).NE.ZERO )THEN
+               TEMP = ALPHA*Y( JY )
+               IX   = KX
+               DO 30, I = 1, M
+                  A( I, J ) = A( I, J ) + X( IX )*TEMP
+                  IX        = IX        + INCX
+   30          CONTINUE
+            END IF
+            JY = JY + INCY
+   40    CONTINUE
+      END IF
+*
+      RETURN
+*
+*     End of DGER  .
+*
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dmach.f
@@ -0,0 +1,58 @@
+      DOUBLE PRECISION FUNCTION DMACH(JOB)
+      INTEGER JOB
+C
+C     SMACH COMPUTES MACHINE PARAMETERS OF FLOATING POINT
+C     ARITHMETIC FOR USE IN TESTING ONLY.  NOT REQUIRED BY
+C     LINPACK PROPER.
+C
+C     IF TROUBLE WITH AUTOMATIC COMPUTATION OF THESE QUANTITIES,
+C     THEY CAN BE SET BY DIRECT ASSIGNMENT STATEMENTS.
+C     ASSUME THE COMPUTER HAS
+C
+C        B = BASE OF ARITHMETIC
+C        T = NUMBER OF BASE  B  DIGITS
+C        L = SMALLEST POSSIBLE EXPONENT
+C        U = LARGEST POSSIBLE EXPONENT
+C
+C     THEN
+C
+C        EPS = B**(1-T)
+C        TINY = 100.0*B**(-L+T)
+C        HUGE = 0.01*B**(U-T)
+C
+C     DMACH SAME AS SMACH EXCEPT T, L, U APPLY TO
+C     DOUBLE PRECISION.
+C
+C     CMACH SAME AS SMACH EXCEPT IF COMPLEX DIVISION
+C     IS DONE BY
+C
+C        1/(X+I*Y) = (X-I*Y)/(X**2+Y**2)
+C
+C     THEN
+C
+C        TINY = SQRT(TINY)
+C        HUGE = SQRT(HUGE)
+C
+C
+C     JOB IS 1, 2 OR 3 FOR EPSILON, TINY AND HUGE, RESPECTIVELY.
+C
+      DOUBLE PRECISION EPS,TINY,HUGE,S
+C
+      EPS = 1.0D0
+   10 EPS = EPS/2.0D0
+      S = 1.0D0 + EPS
+      IF (S .GT. 1.0D0) GO TO 10
+      EPS = 2.0D0*EPS
+C
+      S = 1.0D0
+   20 TINY = S
+      S = S/16.0D0
+      IF (S*1.0 .NE. 0.0D0) GO TO 20
+      TINY = (TINY/EPS)*100.0
+      HUGE = 1.0D0/TINY
+C
+      IF (JOB .EQ. 1) DMACH = EPS
+      IF (JOB .EQ. 2) DMACH = TINY
+      IF (JOB .EQ. 3) DMACH = HUGE
+      RETURN
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dnrm2.f
@@ -0,0 +1,122 @@
+      DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
+      INTEGER          NEXT
+      DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
+      DATA   ZERO, ONE /0.0D0, 1.0D0/
+C
+C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
+C     INCREMENT INCX .
+C     IF    N .LE. 0 RETURN WITH RESULT = 0.
+C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
+C
+C           C.L.LAWSON, 1978 JAN 08
+C
+C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
+C     HOPEFULLY APPLICABLE TO ALL MACHINES.
+C         CUTLO = MAXIMUM OF  DSQRT(U/EPS)  OVER ALL KNOWN MACHINES.
+C         CUTHI = MINIMUM OF  DSQRT(V)      OVER ALL KNOWN MACHINES.
+C     WHERE
+C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
+C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
+C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
+C
+C     BRIEF OUTLINE OF ALGORITHM..
+C
+C     PHASE 1    SCANS ZERO COMPONENTS.
+C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
+C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
+C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
+C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
+C
+C     VALUES FOR CUTLO AND CUTHI..
+C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
+C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
+C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
+C                   UNIVAC AND DEC AT 2**(-103)
+C                   THUS CUTLO = 2**(-51) = 4.44089E-16
+C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
+C                   THUS CUTHI = 2**(63.5) = 1.30438E19
+C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
+C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
+C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
+C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
+C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
+      DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
+C
+      IF(N .GT. 0) GO TO 10
+         DNRM2  = ZERO
+         GO TO 300
+C
+   10 ASSIGN 30 TO NEXT
+      SUM = ZERO
+      NN = N * INCX
+C                                                 BEGIN MAIN LOOP
+      I = 1
+   20    GO TO NEXT,(30, 50, 70, 110)
+   30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
+      ASSIGN 50 TO NEXT
+      XMAX = ZERO
+C
+C                        PHASE 1.  SUM IS ZERO
+C
+   50 IF( DX(I) .EQ. ZERO) GO TO 200
+      IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
+C
+C                                PREPARE FOR PHASE 2.
+      ASSIGN 70 TO NEXT
+      GO TO 105
+C
+C                                PREPARE FOR PHASE 4.
+C
+  100 I = J
+      ASSIGN 110 TO NEXT
+      SUM = (SUM / DX(I)) / DX(I)
+  105 XMAX = DABS(DX(I))
+      GO TO 115
+C
+C                   PHASE 2.  SUM IS SMALL.
+C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
+C
+   70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
+C
+C                     COMMON CODE FOR PHASES 2 AND 4.
+C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
+C
+  110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
+         SUM = ONE + SUM * (XMAX / DX(I))**2
+         XMAX = DABS(DX(I))
+         GO TO 200
+C
+  115 SUM = SUM + (DX(I)/XMAX)**2
+      GO TO 200
+C
+C
+C                  PREPARE FOR PHASE 3.
+C
+   75 SUM = (SUM * XMAX) * XMAX
+C
+C
+C     FOR REAL OR D.P. SET HITEST = CUTHI/N
+C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
+C
+   85 HITEST = CUTHI/FLOAT( N )
+C
+C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
+C
+      DO 95 J =I,NN,INCX
+      IF(DABS(DX(J)) .GE. HITEST) GO TO 100
+   95    SUM = SUM + DX(J)**2
+      DNRM2 = DSQRT( SUM )
+      GO TO 300
+C
+  200 CONTINUE
+      I = I + INCX
+      IF ( I .LE. NN ) GO TO 20
+C
+C              END OF MAIN LOOP.
+C
+C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
+C
+      DNRM2 = XMAX * DSQRT(SUM)
+  300 CONTINUE
+      RETURN
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/drot.f
@@ -0,0 +1,36 @@
+      SUBROUTINE  DROT (N,DX,INCX,DY,INCY,C,S)
+C
+C     APPLIES A PLANE ROTATION.
+C     JACK DONGARRA, LINPACK, 3/11/78.
+C
+      DOUBLE PRECISION DX(1),DY(1),DTEMP,C,S
+      INTEGER I,INCX,INCY,IX,IY,N
+C
+      IF(N.LE.0)RETURN
+      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
+C
+C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
+C         TO 1
+C
+      IX = 1
+      IY = 1
+      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
+      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
+      DO 10 I = 1,N
+        DTEMP = C*DX(IX) + S*DY(IY)
+        DY(IY) = C*DY(IY) - S*DX(IX)
+        DX(IX) = DTEMP
+        IX = IX + INCX
+        IY = IY + INCY
+   10 CONTINUE
+      RETURN
+C
+C       CODE FOR BOTH INCREMENTS EQUAL TO 1
+C
+   20 DO 30 I = 1,N
+        DTEMP = C*DX(I) + S*DY(I)
+        DY(I) = C*DY(I) - S*DX(I)
+        DX(I) = DTEMP
+   30 CONTINUE
+      RETURN
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dscal.f
@@ -0,0 +1,41 @@
+      SUBROUTINE  DSCAL(N,DA,DX,INCX)
+C
+C     SCALES A VECTOR BY A CONSTANT.
+C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
+C     JACK DONGARRA, LINPACK, 3/11/78.
+C
+      DOUBLE PRECISION DA,DX(1)
+      INTEGER I,INCX,M,MP1,N,NINCX
+C
+      IF(N.LE.0)RETURN
+      IF(INCX.EQ.1)GO TO 20
+C
+C        CODE FOR INCREMENT NOT EQUAL TO 1
+C
+      NINCX = N*INCX
+      DO 10 I = 1,NINCX,INCX
+        DX(I) = DA*DX(I)
+   10 CONTINUE
+      RETURN
+C
+C        CODE FOR INCREMENT EQUAL TO 1
+C
+C
+C        CLEAN-UP LOOP
+C
+   20 M = MOD(N,5)
+      IF( M .EQ. 0 ) GO TO 40
+      DO 30 I = 1,M
+        DX(I) = DA*DX(I)
+   30 CONTINUE
+      IF( N .LT. 5 ) RETURN
+   40 MP1 = M + 1
+      DO 50 I = MP1,N,5
+        DX(I) = DA*DX(I)
+        DX(I + 1) = DA*DX(I + 1)
+        DX(I + 2) = DA*DX(I + 2)
+        DX(I + 3) = DA*DX(I + 3)
+        DX(I + 4) = DA*DX(I + 4)
+   50 CONTINUE
+      RETURN
+      END

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dswap.f
@@ -0,0 +1,55 @@
+      SUBROUTINE  DSWAP (N,DX,INCX,DY,INCY)
+C
+C     INTERCHANGES TWO VECTORS.
+C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
+C     JACK DONGARRA, LINPACK, 3/11/78.
+C
+      DOUBLE PRECISION DX(1),DY(1),DTEMP
+      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
+C
+      IF(N.LE.0)RETURN
+      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
+C
+C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
+C         TO 1
+C
+      IX = 1
+      IY = 1
+      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
+      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
+      DO 10 I = 1,N
+        DTEMP = DX(IX)
+        DX(IX) = DY(IY)
+        DY(IY) = DTEMP
+        IX = IX + INCX
+        IY = IY + INCY
+   10 CONTINUE
+      RETURN
+C
+C       CODE FOR BOTH INCREMENTS EQUAL TO 1
+C
+C
+C       CLEAN-UP LOOP
+C
+   20 M = MOD(N,3)
+      IF( M .EQ. 0 ) GO TO 40
+      DO 30 I = 1,M
+        DTEMP = DX(I)
+        DX(I) = DY(I)
+        DY(I) = DTEMP
+   30 CONTINUE
+      IF( N .LT. 3 ) RETURN
+   40 MP1 = M + 1
+      DO 50 I = MP1,N,3
+        DTEMP = DX(I)
+        DX(I) = DY(I)
+        DY(I) = DTEMP
+        DTEMP = DX(I + 1)
+        DX(I + 1) = DY(I + 1)
+        DY(I + 1) = DTEMP
+        DTEMP = DX(I + 2)
+        DX(I + 2) = DY(I + 2)
+        DY(I + 2) = DTEMP
+   50 CONTINUE
+      RETURN
+      END