Repatch (from linux) of OSX IRAF

40 files changed, 3820 insertions, 0 deletions

diff --git a/math/llsq/README b/math/llsq/README
new file mode 100644
index 00000000..bfaabb58
--- /dev/null
+++ b/math/llsq/README
@@ -0,0 +1,33 @@
+This directory contains a collection of routines for solving linear least
+squares problems by the Singular Value Decomposition (SVD) method, as
+described in "Solving Least Squares Problems", by Charles L. Lawson and
+Richard J. Hanson, Prentice Hall, 1974.  The appendix of this book contains
+full listings of the Fortran codes, as well as a users guide.
+
+The numerical subroutines are in this directory.  The directory "progs"
+contains a number of examples of the use of these routines in Fortran programs.
+
+The numerical routines have been modified to eliminate the use of Fortran
+i/o for error conditions.  The integer status return IER has been added to
+all such routines, and the Fortran write statement(s) removed.  A successful
+call returns zero in IER, while an unsucessful call returns a positive integer
+error code, identifying the error.  The original codes are in the directory
+"original_f".
+
+The affected routines and the new calling sequences are as follows:
+
+      subroutine BNDSOL (mode,g,mdg,nb,ip,ir,x,n,rnorm,ier)
+      subroutine LDP	(g,mdg,m,n,h,x,xnorm,w,index,ier)
+      subroutine NNLS	(a,mda,m,n,b,x,rnorm,w,zz,index,ier)
+      subroutine SVDRS	(a,mda,mm,nn,b,mdb,nb,s,ier)
+
+The routines SVA and MFEOUT were not installed in the library, since they
+do extensive i/o, but have been modified to reflect the changes to the
+above subroutines.
+
+See "lsq.x" and "band.x" in progs/ for examples demonstrating the use
+of these routines in IRAF spp programs.
+
+20Nov82 D.Tody
+
+Oct85	Added comma after P edit descriptor in sva.f
diff --git a/math/llsq/bndacc.f b/math/llsq/bndacc.f
new file mode 100644
index 00000000..80b9914c
--- /dev/null
+++ b/math/llsq/bndacc.f
@@ -0,0 +1,74 @@
+      subroutine bndacc (g,mdg,nb,ip,ir,mt,jt)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   sequential algorithm for banded least squares problem..
+c	   accumulation phase.	    for solution phase use bndsol.
+c
+c     the calling program must set ir=1 and ip=1 before the first call
+c     to bndacc for a new case.
+c
+c     the second subscript of g( ) must be dimensioned at least
+c     nb+1 in the calling program.
+      dimension g(mdg,1)
+      zero=0.
+c
+c	       alg. steps 1-4 are performed external to this subroutine.
+c
+      nbp1=nb+1
+      if (mt.le.0) return
+c					      alg. step 5
+      if (jt.eq.ip) go to 70
+c					      alg. steps 6-7
+      if (jt.le.ir) go to 30
+c					      alg. steps 8-9
+      do 10 i=1,mt
+	ig1=jt+mt-i
+	ig2=ir+mt-i
+	do 10 j=1,nbp1
+   10	g(ig1,j)=g(ig2,j)
+c					      alg. step 10
+      ie=jt-ir
+      do 20 i=1,ie
+	ig=ir+i-1
+	do 20 j=1,nbp1
+   20	g(ig,j)=zero
+c					      alg. step 11
+      ir=jt
+c					      alg. step 12
+   30 mu=min0(nb-1,ir-ip-1)
+      if (mu.eq.0) go to 60
+c					      alg. step 13
+      do 50 l=1,mu
+c					      alg. step 14
+	k=min0(l,jt-ip)
+c					      alg. step 15
+	lp1=l+1
+	ig=ip+l
+	do 40 i=lp1,nb
+	  jg=i-k
+   40	  g(ig,jg)=g(ig,i)
+c					      alg. step 16
+	do 50 i=1,k
+	jg=nbp1-i
+   50	g(ig,jg)=zero
+c					      alg. step 17
+   60 ip=jt
+c					      alg. steps 18-19
+   70 mh=ir+mt-ip
+      kh=min0(nbp1,mh)
+c					      alg. step 20
+      do 80 i=1,kh
+   80	call h12 (1,i,max0(i+1,ir-ip+1),mh,g(ip,i),1,rho,
+     1		  g(ip,i+1),1,mdg,nbp1-i)
+c					      alg. step 21
+      ir=ip+kh
+c					      alg. step 22
+      if (kh.lt.nbp1) go to 100
+c					      alg. step 23
+      do 90 i=1,nb
+   90	g(ir-1,i)=zero
+c					      alg. step 24
+  100 continue
+c					      alg. step 25
+      return
+      end
diff --git a/math/llsq/bndsol.f b/math/llsq/bndsol.f
new file mode 100644
index 00000000..4c4bd973
--- /dev/null
+++ b/math/llsq/bndsol.f
@@ -0,0 +1,71 @@
+      subroutine bndsol (mode,g,mdg,nb,ip,ir,x,n,rnorm,ier)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   sequential solution of a banded least squares problem..
+c	   solution phase.   for the accumulation phase use bndacc.
+c
+c     mode = 1	   solve r*x=y	 where r and y are in the g( ) array
+c		   and x will be stored in the x( ) array.
+c	     2	   solve (r**t)*x=y   where r is in g( ),
+c		   y is initially in x( ), and x replaces y in x( ),
+c	     3	   solve r*x=y	 where r is in g( ).
+c		   y is initially in x( ), and x replaces y in x( ).
+c
+c     the second subscript of g( ) must be dimensioned at least
+c     nb+1 in the calling program.
+      dimension g(mdg,1),x(n)
+      zero=0.
+      ier = 0
+c
+      rnorm=zero
+      go to (10,90,50), mode
+c				    ********************* mode = 1
+c				    algc step 26
+   10	   do 20 j=1,n
+   20	   x(j)=g(j,nb+1)
+      rsq=zero
+      np1=n+1
+      irm1=ir-1
+      if (np1.gt.irm1) go to 40
+	   do 30 j=np1,irm1
+   30	   rsq=rsq+g(j,nb+1)**2
+      rnorm=sqrt(rsq)
+   40 continue
+c				    ********************* mode = 3
+c				    alg. step 27
+   50	   do 80 ii=1,n
+	   i=n+1-ii
+c				    alg. step 28
+	   s=zero
+	   l=max0(0,i-ip)
+c				    alg. step 29
+	   if (i.eq.n) go to 70
+c				    alg. step 30
+	   ie=min0(n+1-i,nb)
+		do 60 j=2,ie
+		jg=j+l
+		ix=i-1+j
+   60		s=s+g(i,jg)*x(ix)
+c				    alg. step 31
+   70	   if (g(i,l+1)) 80,130,80
+   80	   x(i)=(x(i)-s)/g(i,l+1)
+c				    alg. step 32
+      return
+c				    ********************* mode = 2
+   90	   do 120 j=1,n
+	   s=zero
+	   if (j.eq.1) go to 110
+	   i1=max0(1,j-nb+1)
+	   i2=j-1
+		do 100 i=i1,i2
+		l=j-i+1+max0(0,i-ip)
+  100		s=s+x(i)*g(i,l)
+  110	   l=max0(0,j-ip)
+	   if (g(j,l+1)) 120,130,120
+  120	   x(j)=(x(j)-s)/g(j,l+1)
+      return
+c
+c     error return: zero diagonal term in bndsol
+  130 ier = 1
+      return
+      end
diff --git a/math/llsq/diff.f b/math/llsq/diff.f
new file mode 100644
index 00000000..79dacb92
--- /dev/null
+++ b/math/llsq/diff.f
@@ -0,0 +1,6 @@
+      function diff(x,y)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 june 7
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+      diff=x-y
+      return
+      end
diff --git a/math/llsq/g1.f b/math/llsq/g1.f
new file mode 100644
index 00000000..5d44c50b
--- /dev/null
+++ b/math/llsq/g1.f
@@ -0,0 +1,33 @@
+      subroutine g1 (a,b,cos,sin,sig)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c
+c     compute orthogonal rotation matrix..
+c     compute.. matrix	 (c, s) so that (c, s)(a) = (sqrt(a**2+b**2))
+c			 (-s,c) 	(-s,c)(b)   (	0	   )
+c     compute sig = sqrt(a**2+b**2)
+c	 sig is computed last to allow for the possibility that
+c	 sig may be in the same location as a or b .
+c
+      zero=0.
+      one=1.
+      if (abs(a).le.abs(b)) go to 10
+      xr=b/a
+      yr=sqrt(one+xr**2)
+      cos=sign(one/yr,a)
+      sin=cos*xr
+      sig=abs(a)*yr
+      return
+   10 if (b) 20,30,20
+   20 xr=a/b
+      yr=sqrt(one+xr**2)
+      sin=sign(one/yr,b)
+      cos=sin*xr
+      sig=abs(b)*yr
+      return
+   30 sig=zero
+      cos=zero
+      sin=one
+      return
+      end
diff --git a/math/llsq/g2.f b/math/llsq/g2.f
new file mode 100644
index 00000000..6f01a582
--- /dev/null
+++ b/math/llsq/g2.f
@@ -0,0 +1,9 @@
+      subroutine g2    (cos,sin,x,y)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1972 dec 15
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   apply the rotation computed by g1 to (x,y).
+      xr=cos*x+sin*y
+      y=-sin*x+cos*y
+      x=xr
+      return
+      end
diff --git a/math/llsq/gen.f b/math/llsq/gen.f
new file mode 100644
index 00000000..e7933d6d
--- /dev/null
+++ b/math/llsq/gen.f
@@ -0,0 +1,31 @@
+      function	 gen(anoise)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1972 dec 15
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   generate numbers for construction of test cases.
+      mi=891
+      mj=457
+      i=5
+      j=7
+      aj=0.
+c
+      if (anoise) 10,30,20
+   10 gen=0.
+      return
+c
+c     the sequence of values of j  is bounded between 1 and 996
+c     if initial j = 1,2,3,4,5,6,7,8, or 9, the period is 332
+c
+   20 j=j*mj
+      j=j-997*(j/997)
+      aj=j-498
+c     the sequence of values of i  is bounded between 1 and 999
+c     if initial i = 1,2,3,6,7, or 9,  the period will be 50
+c     if initial i = 4 or 8   the period will be 25
+c     if initial i = 5	      the period will be 10
+c
+   30 i=i*mi
+      i=i-1000*(i/1000)
+      ai=i-500
+      gen=ai+aj*anoise
+      return
+      end
diff --git a/math/llsq/h12.f b/math/llsq/h12.f
new file mode 100644
index 00000000..81daa629
--- /dev/null
+++ b/math/llsq/h12.f
@@ -0,0 +1,80 @@
+c     subroutine h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c     construction and/or application of a single
+c     householder transformation..     q = i + u*(u**t)/b
+c
+c     mode    = 1 or 2	 to select algorithm  h1  or  h2 .
+c     lpivot is the index of the pivot element.
+c     l1,m   if l1 .le. m   the transformation will be constructed to
+c	     zero elements indexed from l1 through m.	if l1 gt. m
+c	     the subroutine does an identity transformation.
+c     u(),iue,up    on entry to h1 u() contains the pivot vector.
+c		    iue is the storage increment between elements.
+c					on exit from h1 u() and up
+c		    contain quantities defining the vector u of the
+c		    householder transformation.   on entry to h2 u()
+c		    and up should contain quantities previously computed
+c		    by h1.  these will not be modified by h2.
+c     c()    on entry to h1 or h2 c() contains a matrix which will be
+c	     regarded as a set of vectors to which the householder
+c	     transformation is to be applied.  on exit c() contains the
+c	     set of transformed vectors.
+c     ice    storage increment between elements of vectors in c().
+c     icv    storage increment between vectors in c().
+c     ncv    number of vectors in c() to be transformed. if ncv .le. 0
+c	     no operations will be done on c().
+c
+      subroutine h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv)
+      dimension u(iue,m), c(1)
+      double precision sm,b
+      one=1.
+c
+      if (0.ge.lpivot.or.lpivot.ge.l1.or.l1.gt.m) return
+      cl=abs(u(1,lpivot))
+      if (mode.eq.2) go to 60
+c			     ****** construct the transformation. ******
+	  do 10 j=l1,m
+   10	  cl=amax1(abs(u(1,j)),cl)
+      if (cl) 130,130,20
+   20 clinv=one/cl
+      sm=(dble(u(1,lpivot))*clinv)**2
+	  do 30 j=l1,m
+   30	  sm=sm+(dble(u(1,j))*clinv)**2
+c			       convert dble. prec. sm to sngl. prec. sm1
+      sm1=sm
+      cl=cl*sqrt(sm1)
+      if (u(1,lpivot)) 50,50,40
+   40 cl=-cl
+   50 up=u(1,lpivot)-cl
+      u(1,lpivot)=cl
+      go to 70
+c	     ****** apply the transformation  i+u*(u**t)/b  to c. ******
+c
+   60 if (cl) 130,130,70
+   70 if (ncv.le.0) return
+      b=dble(up)*u(1,lpivot)
+c			b  must be nonpositive here.  if b = 0., return.
+c
+      if (b) 80,130,130
+   80 b=one/b
+      i2=1-icv+ice*(lpivot-1)
+      incr=ice*(l1-lpivot)
+	  do 120 j=1,ncv
+	  i2=i2+icv
+	  i3=i2+incr
+	  i4=i3
+	  sm=c(i2)*dble(up)
+	      do 90 i=l1,m
+	      sm=sm+c(i3)*dble(u(1,i))
+   90	      i3=i3+ice
+	  if (sm) 100,120,100
+  100	  sm=sm*b
+	  c(i2)=c(i2)+sm*dble(up)
+	      do 110 i=l1,m
+	      c(i4)=c(i4)+sm*dble(u(1,i))
+  110	      i4=i4+ice
+  120	  continue
+  130 return
+      end
diff --git a/math/llsq/hfti.f b/math/llsq/hfti.f
new file mode 100644
index 00000000..62dcebe4
--- /dev/null
+++ b/math/llsq/hfti.f
@@ -0,0 +1,136 @@
+      subroutine hfti (a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   solve least squares problem using algorithm, hfti.
+c
+      dimension a(mda,n),b(mdb,nb),h(n),g(n),rnorm(nb)
+      integer ip(n)
+      double precision sm,dzero
+      szero=0.
+      dzero=0.d0
+      factor=0.001
+c
+      k=0
+      ldiag=min0(m,n)
+      if (ldiag.le.0) go to 270
+	  do 80 j=1,ldiag
+	  if (j.eq.1) go to 20
+c
+c     update squared column lengths and find lmax
+c    ..
+	  lmax=j
+	      do 10 l=j,n
+	      h(l)=h(l)-a(j-1,l)**2
+	      if (h(l).gt.h(lmax)) lmax=l
+   10	      continue
+	  if(diff(hmax+factor*h(lmax),hmax)) 20,20,50
+c
+c     compute squared column lengths and find lmax
+c    ..
+   20	  lmax=j
+	      do 40 l=j,n
+	      h(l)=0.
+		  do 30 i=j,m
+   30		  h(l)=h(l)+a(i,l)**2
+	      if (h(l).gt.h(lmax)) lmax=l
+   40	      continue
+	  hmax=h(lmax)
+c    ..
+c     lmax has been determined
+c
+c     do column interchanges if needed.
+c    ..
+   50	  continue
+	  ip(j)=lmax
+	  if (ip(j).eq.j) go to 70
+	      do 60 i=1,m
+	      tmp=a(i,j)
+	      a(i,j)=a(i,lmax)
+   60	      a(i,lmax)=tmp
+	  h(lmax)=h(j)
+c
+c     compute the j-th transformation and apply it to a and b.
+c    ..
+   70	  call h12 (1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j)
+   80	  call h12 (2,j,j+1,m,a(1,j),1,h(j),b,1,mdb,nb)
+c
+c     determine the pseudorank, k, using the tolerance, tau.
+c    ..
+	  do 90 j=1,ldiag
+	  if (abs(a(j,j)).le.tau) go to 100
+   90	  continue
+      k=ldiag
+      go to 110
+  100 k=j-1
+  110 kp1=k+1
+c
+c     compute the norms of the residual vectors.
+c
+      if (nb.le.0) go to 140
+	  do 130 jb=1,nb
+	  tmp=szero
+	  if (kp1.gt.m) go to 130
+	      do 120 i=kp1,m
+  120	      tmp=tmp+b(i,jb)**2
+  130	  rnorm(jb)=sqrt(tmp)
+  140 continue
+c					    special for pseudorank = 0
+      if (k.gt.0) go to 160
+      if (nb.le.0) go to 270
+	  do 150 jb=1,nb
+	      do 150 i=1,n
+  150	      b(i,jb)=szero
+      go to 270
+c
+c     if the pseudorank is less than n compute householder
+c     decomposition of first k rows.
+c    ..
+  160 if (k.eq.n) go to 180
+	  do 170 ii=1,k
+	  i=kp1-ii
+  170	  call h12 (1,i,kp1,n,a(i,1),mda,g(i),a,mda,1,i-1)
+  180 continue
+c
+c
+      if (nb.le.0) go to 270
+	  do 260 jb=1,nb
+c
+c     solve the k by k triangular system.
+c    ..
+	      do 210 l=1,k
+	      sm=dzero
+	      i=kp1-l
+	      if (i.eq.k) go to 200
+	      ip1=i+1
+		  do 190 j=ip1,k
+  190		  sm=sm+a(i,j)*dble(b(j,jb))
+  200	      sm1=sm
+  210	      b(i,jb)=(b(i,jb)-sm1)/a(i,i)
+c
+c     complete computation of solution vector.
+c    ..
+	  if (k.eq.n) go to 240
+	      do 220 j=kp1,n
+  220	      b(j,jb)=szero
+	      do 230 i=1,k
+  230	      call h12 (2,i,kp1,n,a(i,1),mda,g(i),b(1,jb),1,mdb,1)
+c
+c      re-order the solution vector to compensate for the
+c      column interchanges.
+c    ..
+  240	      do 250 jj=1,ldiag
+	      j=ldiag+1-jj
+	      if (ip(j).eq.j) go to 250
+	      l=ip(j)
+	      tmp=b(l,jb)
+	      b(l,jb)=b(j,jb)
+	      b(j,jb)=tmp
+  250	      continue
+  260	  continue
+c    ..
+c     the solution vectors, x, are now
+c     in the first  n  rows of the array b(,).
+c
+  270 krank=k
+      return
+      end
diff --git a/math/llsq/ldp.f b/math/llsq/ldp.f
new file mode 100644
index 00000000..a083f35b
--- /dev/null
+++ b/math/llsq/ldp.f
@@ -0,0 +1,79 @@
+      subroutine ldp (g,mdg,m,n,h,x,xnorm,w,index,ier)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1974 mar 1
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c	     **********  least distance programming  **********
+c
+      integer index(m)
+      dimension g(mdg,n), h(m), x(n), w(1)
+      zero=0.
+      one=1.
+      if (n.le.0) go to 120
+	  do 10 j=1,n
+   10	  x(j)=zero
+      xnorm=zero
+      if (m.le.0) go to 110
+c
+c     the declared dimension of w() must be at least (n+1)*(m+2)+2*m.
+c
+c      first (n+1)*m locs of w()   =  matrix e for problem nnls.
+c	next	 n+1 locs of w()   =  vector f for problem nnls.
+c	next	 n+1 locs of w()   =  vector z for problem nnls.
+c	next	   m locs of w()   =  vector y for problem nnls.
+c	next	   m locs of w()   =  vector wdual for problem nnls.
+c     copy g**t into first n rows and m columns of e.
+c     copy h**t into row n+1 of e.
+c
+      iw=0
+	  do 30 j=1,m
+	      do 20 i=1,n
+	      iw=iw+1
+   20	      w(iw)=g(j,i)
+	  iw=iw+1
+   30	  w(iw)=h(j)
+      if=iw+1
+c				 store n zeros followed by a one into f.
+	  do 40 i=1,n
+	  iw=iw+1
+   40	  w(iw)=zero
+      w(iw+1)=one
+c
+      np1=n+1
+      iz=iw+2
+      iy=iz+np1
+      iwdual=iy+m
+c
+      call nnls (w,np1,np1,m,w(if),w(iy),rnorm,w(iwdual),w(iz),index,
+     *		 ier)
+c		       use the following return if unsuccessful in nnls.
+      if (ier.ne.0) return
+      if (rnorm) 130,130,50
+   50 fac=one
+      iw=iy-1
+	  do 60 i=1,m
+	  iw=iw+1
+c				here we are using the solution vector y.
+   60	  fac=fac-h(i)*w(iw)
+c
+      if (diff(one+fac,one)) 130,130,70
+   70 fac=one/fac
+	  do 90 j=1,n
+	  iw=iy-1
+	      do 80 i=1,m
+	      iw=iw+1
+c				here we are using the solution vector y.
+   80	      x(j)=x(j)+g(i,j)*w(iw)
+   90	  x(j)=x(j)*fac
+	  do 100 j=1,n
+  100	  xnorm=xnorm+x(j)**2
+      xnorm=sqrt(xnorm)
+c			      successful return.
+  110 ier=0
+      return
+c			      error return.	  n .le. 0.
+  120 ier=2
+      return
+c			      returning with constraints not compatible.
+  130 ier=4
+      return
+      end
diff --git a/math/llsq/mfeout.f b/math/llsq/mfeout.f
new file mode 100644
index 00000000..fb439e40
--- /dev/null
+++ b/math/llsq/mfeout.f
@@ -0,0 +1,64 @@
+      subroutine mfeout (a,mda,m,n,names,mode)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   subroutine for matrix output with labeling.
+c
+c     a( )	   matrix to be output
+c		   mda	   first dimension of a array
+c		   m	     no. of rows in a matrix
+c		   n	     no. of cols in a matrix
+c     names()	   array of names.  if names(1) = 1h , the rest
+c		   of the names() array will be ignored.
+c     mode	   =1	for   4p8f15.0	format	for v matrix.
+c		   =2	for   8e15.8  format  for candidate solutions.
+c
+      dimension    a(mda,01)
+      integer names(m),ihead(2)
+      logical	notblk
+      data  maxcol/8/, iblank/1h /,ihead(1)/4h col/,ihead(2)/4hsoln/
+c
+      notblk=names(1).ne.iblank
+      if (m.le.0.or.n.le.0) return
+c
+      if (mode.eq.2) go to 10
+      write (6,70)
+      go to 20
+   10 write (6,80)
+   20 continue
+c
+      nblock=n/maxcol
+      last=n-nblock*maxcol
+      ncol=maxcol
+      j1=1
+c
+c			     main loop starts here
+c
+   30 if (nblock.gt.0) go to 40
+      if (last.le.0) return
+      ncol=last
+      last=0
+c
+   40 j2=j1+ncol-1
+      write (6,90) (ihead(mode),j,j=j1,j2)
+c
+	   do 60 i=1,m
+	   name=iblank
+	   if (notblk) name=names(i)
+c
+	   if (mode.eq.2) go to 50
+	   write (6,100) i,name,(a(i,j),j=j1,j2)
+	   go to 60
+   50	   write (6,110) i,name,(a(i,j),j=j1,j2)
+   60	   continue
+c
+      j1=j1+maxcol
+      nblock=nblock-1
+      go to 30
+c
+   70 format (45h0v-matrix of the singular value decomposition,
+     * 8h of a*d./47h (elements of v scaled up by a factor of 10**4))
+   80 format (35h0sequence of candidate solutions, x)
+   90 format (1h0,11x,8(6x,a4,i4,1x)/1x)
+  100 format (1x,i3,1x,a6,1x,4p8f15.0)
+  110 format (1x,i3,1x,a6,1x,8e15.8)
+      end
diff --git a/math/llsq/mkpkg b/math/llsq/mkpkg
new file mode 100644
index 00000000..1d2cf4a5
--- /dev/null
+++ b/math/llsq/mkpkg
@@ -0,0 +1,23 @@
+# Update Lawson's and Hanson's linear least squares package (LIBLLSQ).
+
+$checkout libllsq.a lib$
+$update   libllsq.a
+$checkin  libllsq.a lib$
+$exit
+
+libllsq.a:
+	bndacc.f
+	bndsol.f
+	diff.f
+	g1.f
+	g2.f
+	gen.f
+	h12.f
+	hfti.f
+	ldp.f
+	mfeout.f
+	nnls.f
+	qrbd.f
+	sva.f
+	svdrs.f
+	;
diff --git a/math/llsq/nnls.f b/math/llsq/nnls.f
new file mode 100644
index 00000000..58337147
--- /dev/null
+++ b/math/llsq/nnls.f
@@ -0,0 +1,276 @@
+c     subroutine nnls  (a,mda,m,n,b,x,rnorm,w,zz,index,ier)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 june 15
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c	  **********   nonnegative least squares   **********
+c
+c     given an m by n matrix, a, and an m-vector, b,  compute an
+c     n-vector, x, which solves the least squares problem
+c
+c		       a * x = b  subject to x .ge. 0
+c
+c     a(),mda,m,n     mda is the first dimensioning parameter for the
+c		      array, a().   on entry a() contains the m by n
+c		      matrix, a.	   on exit a() contains
+c		      the product matrix, q*a , where q is an
+c		      m by m orthogonal matrix generated implicitly by
+c		      this subroutine.
+c     b()     on entry b() contains the m-vector, b.   on exit b() con-
+c	      tains q*b.
+c     x()     on entry x() need not be initialized.  on exit x() will
+c	      contain the solution vector.
+c     rnorm   on exit rnorm contains the euclidean norm of the
+c	      residual vector.
+c     w()     an n-array of working space.  on exit w() will contain
+c	      the dual solution vector.   w will satisfy w(i) = 0.
+c	      for all i in set p  and w(i) .le. 0. for all i in set z
+c     zz()     an m-array of working space.
+c     index()	  an integer working array of length at least n.
+c		  on exit the contents of this array define the sets
+c		  p and z as follows..
+c
+c		  index(1)   thru index(nsetp) = set p.
+c		  index(iz1) thru index(iz2)   = set z.
+c		  iz1 = nsetp + 1 = npp1
+c		  iz2 = n
+c     ier     this is a success-failure flag with the following
+c	      meanings.
+c	      0     the solution has been computed successfully.
+c	      2     the dimensions of the problem are bad.
+c		    either m .le. 0 or n .le. 0.
+c	      3     iteration count exceeded.  more than 3*n iterations.
+c
+      subroutine nnls (a,mda,m,n,b,x,rnorm,w,zz,index,ier)
+      dimension a(mda,n), b(m), x(n), w(n), zz(m)
+      integer index(n)
+      zero=0.
+      one=1.
+      two=2.
+      factor=0.01
+c
+      ier=0
+      if (m.gt.0.and.n.gt.0) go to 10
+      ier=2
+      return
+   10 iter=0
+      itmax=3*n
+c
+c		     initialize the arrays index() and x().
+c
+	  do 20 i=1,n
+	  x(i)=zero
+   20	  index(i)=i
+c
+      iz2=n
+      iz1=1
+      nsetp=0
+      npp1=1
+c			      ******  main loop begins here  ******
+   30 continue
+c		   quit if all coefficients are already in the solution.
+c			 or if m cols of a have been triangularized.
+c
+      if (iz1.gt.iz2.or.nsetp.ge.m) go to 350
+c
+c	  compute components of the dual (negative gradient) vector w().
+c
+	  do 50 iz=iz1,iz2
+	  j=index(iz)
+	  sm=zero
+	      do 40 l=npp1,m
+   40	      sm=sm+a(l,j)*b(l)
+   50	  w(j)=sm
+c				    find largest positive w(j).
+   60 wmax=zero
+	  do 70 iz=iz1,iz2
+	  j=index(iz)
+	  if (w(j).le.wmax) go to 70
+	  wmax=w(j)
+	  izmax=iz
+   70	  continue
+c
+c	      if wmax .le. 0. go to termination.
+c	      this indicates satisfaction of the kuhn-tucker conditions.
+c
+      if (wmax) 350,350,80
+   80 iz=izmax
+      j=index(iz)
+c
+c     the sign of w(j) is ok for j to be moved to set p.
+c     begin the transformation and check new diagonal element to avoid
+c     near linear dependence.
+c
+      asave=a(npp1,j)
+      call h12 (1,npp1,npp1+1,m,a(1,j),1,up,dummy,1,1,0)
+      unorm=zero
+      if (nsetp.eq.0) go to 100
+	  do 90 l=1,nsetp
+   90	  unorm=unorm+a(l,j)**2
+  100 unorm=sqrt(unorm)
+      if (diff(unorm+abs(a(npp1,j))*factor,unorm)) 130,130,110
+c
+c     col j is sufficiently independent.  copy b into zz, update zz and
+c   > solve for ztest ( = proposed new value for x(j) ).
+c
+  110	  do 120 l=1,m
+  120	  zz(l)=b(l)
+      call h12 (2,npp1,npp1+1,m,a(1,j),1,up,zz,1,1,1)
+      ztest=zz(npp1)/a(npp1,j)
+c
+c				      see if ztest is positive
+c     reject j as a candidate to be moved from set z to set p.
+c     restore a(npp1,j), set w(j)=0., and loop back to test dual
+c
+      if (ztest) 130,130,140
+c
+c     coeffs again.
+c
+  130 a(npp1,j)=asave
+      w(j)=zero
+      go to 60
+c
+c     the index  j=index(iz)  has been selected to be moved from
+c     set z to set p.	 update b,  update indices,  apply householder
+c     transformations to cols in new set z,  zero subdiagonal elts in
+c     col j,  set w(j)=0.
+c
+  140	  do 150 l=1,m
+  150	  b(l)=zz(l)
+c
+      index(iz)=index(iz1)
+      index(iz1)=j
+      iz1=iz1+1
+      nsetp=npp1
+      npp1=npp1+1
+c
+      if (iz1.gt.iz2) go to 170
+	  do 160 jz=iz1,iz2
+	  jj=index(jz)
+  160	  call h12 (2,nsetp,npp1,m,a(1,j),1,up,a(1,jj),1,mda,1)
+  170 continue
+c
+      if (nsetp.eq.m) go to 190
+	  do 180 l=npp1,m
+  180	  a(l,j)=zero
+  190 continue
+c
+      w(j)=zero
+c				 solve the triangular system.
+c				 store the solution temporarily in zz().
+      assign 200 to next
+      go to 400
+  200 continue
+c
+c			******	secondary loop begins here ******
+c
+c			   iteration counter.
+c
+  210 iter=iter+1
+      if (iter.le.itmax) go to 220
+         ier=3
+         go to 350
+  220 continue
+c
+c		     see if all new constrained coeffs are feasible.
+c				   if not compute alpha.
+c
+      alpha=two
+	  do 240 ip=1,nsetp
+	  l=index(ip)
+	  if (zz(ip)) 230,230,240
+c
+  230	  t=-x(l)/(zz(ip)-x(l))
+	  if (alpha.le.t) go to 240
+	  alpha=t
+	  jj=ip
+  240	  continue
+c
+c	   if all new constrained coeffs are feasible then alpha will
+c	   still = 2.	 if so exit from secondary loop to main loop.
+c
+      if (alpha.eq.two) go to 330
+c
+c	   otherwise use alpha which will be between 0. and 1. to
+c	   interpolate between the old x and the new zz.
+c
+	  do 250 ip=1,nsetp
+	  l=index(ip)
+  250	  x(l)=x(l)+alpha*(zz(ip)-x(l))
+c
+c	 modify a and b and the index arrays to move coefficient i
+c	 from set p to set z.
+c
+      i=index(jj)
+  260 x(i)=zero
+c
+      if (jj.eq.nsetp) go to 290
+      jj=jj+1
+	  do 280 j=jj,nsetp
+	  ii=index(j)
+	  index(j-1)=ii
+	  call g1 (a(j-1,ii),a(j,ii),cc,ss,a(j-1,ii))
+	  a(j,ii)=zero
+	      do 270 l=1,n
+	      if (l.ne.ii) call g2 (cc,ss,a(j-1,l),a(j,l))
+  270	      continue
+  280	  call g2 (cc,ss,b(j-1),b(j))
+  290 npp1=nsetp
+      nsetp=nsetp-1
+      iz1=iz1-1
+      index(iz1)=i
+c
+c	 see if the remaining coeffs in set p are feasible.  they should
+c	 be because of the way alpha was determined.
+c	 if any are infeasible it is due to round-off error.  any
+c	 that are nonpositive will be set to zero
+c	 and moved from set p to set z.
+c
+	  do 300 jj=1,nsetp
+	  i=index(jj)
+	  if (x(i)) 260,260,300
+  300	  continue
+c
+c	  copy b( ) into zz( ).  then solve again and loop back.
+c
+
+	  do 310 i=1,m
+  310	  zz(i)=b(i)
+      assign 320 to next
+      go to 400
+  320 continue
+      go to 210
+c		       ******  end of secondary loop  ******
+c
+  330	  do 340 ip=1,nsetp
+	  i=index(ip)
+  340	  x(i)=zz(ip)
+c	 all new coeffs are positive.  loop back to beginning.
+      go to 30
+c
+c			 ******  end of main loop  ******
+c
+c			 come to here for termination.
+c		      compute the norm of the final residual vector.
+c
+  350 sm=zero
+      if (npp1.gt.m) go to 370
+	  do 360 i=npp1,m
+  360	  sm=sm+b(i)**2
+      go to 390
+  370	  do 380 j=1,n
+  380	  w(j)=zero
+  390 rnorm=sqrt(sm)
+      return
+c
+c     the following block of code is used as an internal subroutine
+c     to solve the triangular system, putting the solution in zz().
+c
+  400	  do 430 l=1,nsetp
+	  ip=nsetp+1-l
+	  if (l.eq.1) go to 420
+	      do 410 ii=1,ip
+  410	      zz(ii)=zz(ii)-a(ii,jj)*zz(ip+1)
+  420	  jj=index(ip)
+  430	  zz(ip)=zz(ip)/a(ip,jj)
+      go to next, (200,320)
+      end
diff --git a/math/llsq/original_f/bndacc.f b/math/llsq/original_f/bndacc.f
new file mode 100644
index 00000000..80b9914c
--- /dev/null
+++ b/math/llsq/original_f/bndacc.f
@@ -0,0 +1,74 @@
+      subroutine bndacc (g,mdg,nb,ip,ir,mt,jt)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   sequential algorithm for banded least squares problem..
+c	   accumulation phase.	    for solution phase use bndsol.
+c
+c     the calling program must set ir=1 and ip=1 before the first call
+c     to bndacc for a new case.
+c
+c     the second subscript of g( ) must be dimensioned at least
+c     nb+1 in the calling program.
+      dimension g(mdg,1)
+      zero=0.
+c
+c	       alg. steps 1-4 are performed external to this subroutine.
+c
+      nbp1=nb+1
+      if (mt.le.0) return
+c					      alg. step 5
+      if (jt.eq.ip) go to 70
+c					      alg. steps 6-7
+      if (jt.le.ir) go to 30
+c					      alg. steps 8-9
+      do 10 i=1,mt
+	ig1=jt+mt-i
+	ig2=ir+mt-i
+	do 10 j=1,nbp1
+   10	g(ig1,j)=g(ig2,j)
+c					      alg. step 10
+      ie=jt-ir
+      do 20 i=1,ie
+	ig=ir+i-1
+	do 20 j=1,nbp1
+   20	g(ig,j)=zero
+c					      alg. step 11
+      ir=jt
+c					      alg. step 12
+   30 mu=min0(nb-1,ir-ip-1)
+      if (mu.eq.0) go to 60
+c					      alg. step 13
+      do 50 l=1,mu
+c					      alg. step 14
+	k=min0(l,jt-ip)
+c					      alg. step 15
+	lp1=l+1
+	ig=ip+l
+	do 40 i=lp1,nb
+	  jg=i-k
+   40	  g(ig,jg)=g(ig,i)
+c					      alg. step 16
+	do 50 i=1,k
+	jg=nbp1-i
+   50	g(ig,jg)=zero
+c					      alg. step 17
+   60 ip=jt
+c					      alg. steps 18-19
+   70 mh=ir+mt-ip
+      kh=min0(nbp1,mh)
+c					      alg. step 20
+      do 80 i=1,kh
+   80	call h12 (1,i,max0(i+1,ir-ip+1),mh,g(ip,i),1,rho,
+     1		  g(ip,i+1),1,mdg,nbp1-i)
+c					      alg. step 21
+      ir=ip+kh
+c					      alg. step 22
+      if (kh.lt.nbp1) go to 100
+c					      alg. step 23
+      do 90 i=1,nb
+   90	g(ir-1,i)=zero
+c					      alg. step 24
+  100 continue
+c					      alg. step 25
+      return
+      end
diff --git a/math/llsq/original_f/bndsol.f b/math/llsq/original_f/bndsol.f
new file mode 100644
index 00000000..05dd35d7
--- /dev/null
+++ b/math/llsq/original_f/bndsol.f
@@ -0,0 +1,70 @@
+      subroutine bndsol (mode,g,mdg,nb,ip,ir,x,n,rnorm)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   sequential solution of a banded least squares problem..
+c	   solution phase.   for the accumulation phase use bndacc.
+c
+c     mode = 1	   solve r*x=y	 where r and y are in the g( ) array
+c		   and x will be stored in the x( ) array.
+c	     2	   solve (r**t)*x=y   where r is in g( ),
+c		   y is initially in x( ), and x replaces y in x( ),
+c	     3	   solve r*x=y	 where r is in g( ).
+c		   y is initially in x( ), and x replaces y in x( ).
+c
+c     the second subscript of g( ) must be dimensioned at least
+c     nb+1 in the calling program.
+      dimension g(mdg,1),x(n)
+      zero=0.
+c
+      rnorm=zero
+      go to (10,90,50), mode
+c				    ********************* mode = 1
+c				    algc step 26
+   10	   do 20 j=1,n
+   20	   x(j)=g(j,nb+1)
+      rsq=zero
+      np1=n+1
+      irm1=ir-1
+      if (np1.gt.irm1) go to 40
+	   do 30 j=np1,irm1
+   30	   rsq=rsq+g(j,nb+1)**2
+      rnorm=sqrt(rsq)
+   40 continue
+c				    ********************* mode = 3
+c				    alg. step 27
+   50	   do 80 ii=1,n
+	   i=n+1-ii
+c				    alg. step 28
+	   s=zero
+	   l=max0(0,i-ip)
+c				    alg. step 29
+	   if (i.eq.n) go to 70
+c				    alg. step 30
+	   ie=min0(n+1-i,nb)
+		do 60 j=2,ie
+		jg=j+l
+		ix=i-1+j
+   60		s=s+g(i,jg)*x(ix)
+c				    alg. step 31
+   70	   if (g(i,l+1)) 80,130,80
+   80	   x(i)=(x(i)-s)/g(i,l+1)
+c				    alg. step 32
+      return
+c				    ********************* mode = 2
+   90	   do 120 j=1,n
+	   s=zero
+	   if (j.eq.1) go to 110
+	   i1=max0(1,j-nb+1)
+	   i2=j-1
+		do 100 i=i1,i2
+		l=j-i+1+max0(0,i-ip)
+  100		s=s+x(i)*g(i,l)
+  110	   l=max0(0,j-ip)
+	   if (g(j,l+1)) 120,130,120
+  120	   x(j)=(x(j)-s)/g(j,l+1)
+      return
+c
+  130 write (6,140) mode,i,j,l
+      stop
+  140 format (30h0zero diagonal term in bndsol.,12h mode,i,j,l=,4i6)
+      end
diff --git a/math/llsq/original_f/diff.f b/math/llsq/original_f/diff.f
new file mode 100644
index 00000000..79dacb92
--- /dev/null
+++ b/math/llsq/original_f/diff.f
@@ -0,0 +1,6 @@
+      function diff(x,y)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 june 7
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+      diff=x-y
+      return
+      end
diff --git a/math/llsq/original_f/g1.f b/math/llsq/original_f/g1.f
new file mode 100644
index 00000000..5d44c50b
--- /dev/null
+++ b/math/llsq/original_f/g1.f
@@ -0,0 +1,33 @@
+      subroutine g1 (a,b,cos,sin,sig)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c
+c     compute orthogonal rotation matrix..
+c     compute.. matrix	 (c, s) so that (c, s)(a) = (sqrt(a**2+b**2))
+c			 (-s,c) 	(-s,c)(b)   (	0	   )
+c     compute sig = sqrt(a**2+b**2)
+c	 sig is computed last to allow for the possibility that
+c	 sig may be in the same location as a or b .
+c
+      zero=0.
+      one=1.
+      if (abs(a).le.abs(b)) go to 10
+      xr=b/a
+      yr=sqrt(one+xr**2)
+      cos=sign(one/yr,a)
+      sin=cos*xr
+      sig=abs(a)*yr
+      return
+   10 if (b) 20,30,20
+   20 xr=a/b
+      yr=sqrt(one+xr**2)
+      sin=sign(one/yr,b)
+      cos=sin*xr
+      sig=abs(b)*yr
+      return
+   30 sig=zero
+      cos=zero
+      sin=one
+      return
+      end
diff --git a/math/llsq/original_f/g2.f b/math/llsq/original_f/g2.f
new file mode 100644
index 00000000..6f01a582
--- /dev/null
+++ b/math/llsq/original_f/g2.f
@@ -0,0 +1,9 @@
+      subroutine g2    (cos,sin,x,y)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1972 dec 15
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   apply the rotation computed by g1 to (x,y).
+      xr=cos*x+sin*y
+      y=-sin*x+cos*y
+      x=xr
+      return
+      end
diff --git a/math/llsq/original_f/gen.f b/math/llsq/original_f/gen.f
new file mode 100644
index 00000000..98181a93
--- /dev/null
+++ b/math/llsq/original_f/gen.f
@@ -0,0 +1,28 @@
+      function	 gen(anoise)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1972 dec 15
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   generate numbers for construction of test cases.
+      if (anoise) 10,30,20
+   10 mi=891
+      mj=457
+      i=5
+      j=7
+      aj=0.
+      gen=0.
+      return
+c
+c     the sequence of values of j  is bounded between 1 and 996
+c     if initial j = 1,2,3,4,5,6,7,8, or 9, the period is 332
+   20 j=j*mj
+      j=j-997*(j/997)
+      aj=j-498
+c     the sequence of values of i  is bounded between 1 and 999
+c     if initial i = 1,2,3,6,7, or 9,  the period will be 50
+c     if initial i = 4 or 8   the period will be 25
+c     if initial i = 5	      the period will be 10
+   30 i=i*mi
+      i=i-1000*(i/1000)
+      ai=i-500
+      gen=ai+aj*anoise
+      return
+      end
diff --git a/math/llsq/original_f/h12.f b/math/llsq/original_f/h12.f
new file mode 100644
index 00000000..81daa629
--- /dev/null
+++ b/math/llsq/original_f/h12.f
@@ -0,0 +1,80 @@
+c     subroutine h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c     construction and/or application of a single
+c     householder transformation..     q = i + u*(u**t)/b
+c
+c     mode    = 1 or 2	 to select algorithm  h1  or  h2 .
+c     lpivot is the index of the pivot element.
+c     l1,m   if l1 .le. m   the transformation will be constructed to
+c	     zero elements indexed from l1 through m.	if l1 gt. m
+c	     the subroutine does an identity transformation.
+c     u(),iue,up    on entry to h1 u() contains the pivot vector.
+c		    iue is the storage increment between elements.
+c					on exit from h1 u() and up
+c		    contain quantities defining the vector u of the
+c		    householder transformation.   on entry to h2 u()
+c		    and up should contain quantities previously computed
+c		    by h1.  these will not be modified by h2.
+c     c()    on entry to h1 or h2 c() contains a matrix which will be
+c	     regarded as a set of vectors to which the householder
+c	     transformation is to be applied.  on exit c() contains the
+c	     set of transformed vectors.
+c     ice    storage increment between elements of vectors in c().
+c     icv    storage increment between vectors in c().
+c     ncv    number of vectors in c() to be transformed. if ncv .le. 0
+c	     no operations will be done on c().
+c
+      subroutine h12 (mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv)
+      dimension u(iue,m), c(1)
+      double precision sm,b
+      one=1.
+c
+      if (0.ge.lpivot.or.lpivot.ge.l1.or.l1.gt.m) return
+      cl=abs(u(1,lpivot))
+      if (mode.eq.2) go to 60
+c			     ****** construct the transformation. ******
+	  do 10 j=l1,m
+   10	  cl=amax1(abs(u(1,j)),cl)
+      if (cl) 130,130,20
+   20 clinv=one/cl
+      sm=(dble(u(1,lpivot))*clinv)**2
+	  do 30 j=l1,m
+   30	  sm=sm+(dble(u(1,j))*clinv)**2
+c			       convert dble. prec. sm to sngl. prec. sm1
+      sm1=sm
+      cl=cl*sqrt(sm1)
+      if (u(1,lpivot)) 50,50,40
+   40 cl=-cl
+   50 up=u(1,lpivot)-cl
+      u(1,lpivot)=cl
+      go to 70
+c	     ****** apply the transformation  i+u*(u**t)/b  to c. ******
+c
+   60 if (cl) 130,130,70
+   70 if (ncv.le.0) return
+      b=dble(up)*u(1,lpivot)
+c			b  must be nonpositive here.  if b = 0., return.
+c
+      if (b) 80,130,130
+   80 b=one/b
+      i2=1-icv+ice*(lpivot-1)
+      incr=ice*(l1-lpivot)
+	  do 120 j=1,ncv
+	  i2=i2+icv
+	  i3=i2+incr
+	  i4=i3
+	  sm=c(i2)*dble(up)
+	      do 90 i=l1,m
+	      sm=sm+c(i3)*dble(u(1,i))
+   90	      i3=i3+ice
+	  if (sm) 100,120,100
+  100	  sm=sm*b
+	  c(i2)=c(i2)+sm*dble(up)
+	      do 110 i=l1,m
+	      c(i4)=c(i4)+sm*dble(u(1,i))
+  110	      i4=i4+ice
+  120	  continue
+  130 return
+      end
diff --git a/math/llsq/original_f/hfti.f b/math/llsq/original_f/hfti.f
new file mode 100644
index 00000000..62dcebe4
--- /dev/null
+++ b/math/llsq/original_f/hfti.f
@@ -0,0 +1,136 @@
+      subroutine hfti (a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   solve least squares problem using algorithm, hfti.
+c
+      dimension a(mda,n),b(mdb,nb),h(n),g(n),rnorm(nb)
+      integer ip(n)
+      double precision sm,dzero
+      szero=0.
+      dzero=0.d0
+      factor=0.001
+c
+      k=0
+      ldiag=min0(m,n)
+      if (ldiag.le.0) go to 270
+	  do 80 j=1,ldiag
+	  if (j.eq.1) go to 20
+c
+c     update squared column lengths and find lmax
+c    ..
+	  lmax=j
+	      do 10 l=j,n
+	      h(l)=h(l)-a(j-1,l)**2
+	      if (h(l).gt.h(lmax)) lmax=l
+   10	      continue
+	  if(diff(hmax+factor*h(lmax),hmax)) 20,20,50
+c
+c     compute squared column lengths and find lmax
+c    ..
+   20	  lmax=j
+	      do 40 l=j,n
+	      h(l)=0.
+		  do 30 i=j,m
+   30		  h(l)=h(l)+a(i,l)**2
+	      if (h(l).gt.h(lmax)) lmax=l
+   40	      continue
+	  hmax=h(lmax)
+c    ..
+c     lmax has been determined
+c
+c     do column interchanges if needed.
+c    ..
+   50	  continue
+	  ip(j)=lmax
+	  if (ip(j).eq.j) go to 70
+	      do 60 i=1,m
+	      tmp=a(i,j)
+	      a(i,j)=a(i,lmax)
+   60	      a(i,lmax)=tmp
+	  h(lmax)=h(j)
+c
+c     compute the j-th transformation and apply it to a and b.
+c    ..
+   70	  call h12 (1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j)
+   80	  call h12 (2,j,j+1,m,a(1,j),1,h(j),b,1,mdb,nb)
+c
+c     determine the pseudorank, k, using the tolerance, tau.
+c    ..
+	  do 90 j=1,ldiag
+	  if (abs(a(j,j)).le.tau) go to 100
+   90	  continue
+      k=ldiag
+      go to 110
+  100 k=j-1
+  110 kp1=k+1
+c
+c     compute the norms of the residual vectors.
+c
+      if (nb.le.0) go to 140
+	  do 130 jb=1,nb
+	  tmp=szero
+	  if (kp1.gt.m) go to 130
+	      do 120 i=kp1,m
+  120	      tmp=tmp+b(i,jb)**2
+  130	  rnorm(jb)=sqrt(tmp)
+  140 continue
+c					    special for pseudorank = 0
+      if (k.gt.0) go to 160
+      if (nb.le.0) go to 270
+	  do 150 jb=1,nb
+	      do 150 i=1,n
+  150	      b(i,jb)=szero
+      go to 270
+c
+c     if the pseudorank is less than n compute householder
+c     decomposition of first k rows.
+c    ..
+  160 if (k.eq.n) go to 180
+	  do 170 ii=1,k
+	  i=kp1-ii
+  170	  call h12 (1,i,kp1,n,a(i,1),mda,g(i),a,mda,1,i-1)
+  180 continue
+c
+c
+      if (nb.le.0) go to 270
+	  do 260 jb=1,nb
+c
+c     solve the k by k triangular system.
+c    ..
+	      do 210 l=1,k
+	      sm=dzero
+	      i=kp1-l
+	      if (i.eq.k) go to 200
+	      ip1=i+1
+		  do 190 j=ip1,k
+  190		  sm=sm+a(i,j)*dble(b(j,jb))
+  200	      sm1=sm
+  210	      b(i,jb)=(b(i,jb)-sm1)/a(i,i)
+c
+c     complete computation of solution vector.
+c    ..
+	  if (k.eq.n) go to 240
+	      do 220 j=kp1,n
+  220	      b(j,jb)=szero
+	      do 230 i=1,k
+  230	      call h12 (2,i,kp1,n,a(i,1),mda,g(i),b(1,jb),1,mdb,1)
+c
+c      re-order the solution vector to compensate for the
+c      column interchanges.
+c    ..
+  240	      do 250 jj=1,ldiag
+	      j=ldiag+1-jj
+	      if (ip(j).eq.j) go to 250
+	      l=ip(j)
+	      tmp=b(l,jb)
+	      b(l,jb)=b(j,jb)
+	      b(j,jb)=tmp
+  250	      continue
+  260	  continue
+c    ..
+c     the solution vectors, x, are now
+c     in the first  n  rows of the array b(,).
+c
+  270 krank=k
+      return
+      end
diff --git a/math/llsq/original_f/ldp.f b/math/llsq/original_f/ldp.f
new file mode 100644
index 00000000..65d4c77e
--- /dev/null
+++ b/math/llsq/original_f/ldp.f
@@ -0,0 +1,79 @@
+      subroutine ldp (g,mdg,m,n,h,x,xnorm,w,index,mode)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1974 mar 1
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c	     **********  least distance programming  **********
+c
+      integer index(m)
+      dimension g(mdg,n), h(m), x(n), w(1)
+      zero=0.
+      one=1.
+      if (n.le.0) go to 120
+	  do 10 j=1,n
+   10	  x(j)=zero
+      xnorm=zero
+      if (m.le.0) go to 110
+c
+c     the declared dimension of w() must be at least (n+1)*(m+2)+2*m.
+c
+c      first (n+1)*m locs of w()   =  matrix e for problem nnls.
+c	next	 n+1 locs of w()   =  vector f for problem nnls.
+c	next	 n+1 locs of w()   =  vector z for problem nnls.
+c	next	   m locs of w()   =  vector y for problem nnls.
+c	next	   m locs of w()   =  vector wdual for problem nnls.
+c     copy g**t into first n rows and m columns of e.
+c     copy h**t into row n+1 of e.
+c
+      iw=0
+	  do 30 j=1,m
+	      do 20 i=1,n
+	      iw=iw+1
+   20	      w(iw)=g(j,i)
+	  iw=iw+1
+   30	  w(iw)=h(j)
+      if=iw+1
+c				 store n zeros followed by a one into f.
+	  do 40 i=1,n
+	  iw=iw+1
+   40	  w(iw)=zero
+      w(iw+1)=one
+c
+      np1=n+1
+      iz=iw+2
+      iy=iz+np1
+      iwdual=iy+m
+c
+      call nnls (w,np1,np1,m,w(if),w(iy),rnorm,w(iwdual),w(iz),index,
+     *		 mode)
+c		       use the following return if unsuccessful in nnls.
+      if (mode.ne.1) return
+      if (rnorm) 130,130,50
+   50 fac=one
+      iw=iy-1
+	  do 60 i=1,m
+	  iw=iw+1
+c				here we are using the solution vector y.
+   60	  fac=fac-h(i)*w(iw)
+c
+      if (diff(one+fac,one)) 130,130,70
+   70 fac=one/fac
+	  do 90 j=1,n
+	  iw=iy-1
+	      do 80 i=1,m
+	      iw=iw+1
+c				here we are using the solution vector y.
+   80	      x(j)=x(j)+g(i,j)*w(iw)
+   90	  x(j)=x(j)*fac
+	  do 100 j=1,n
+  100	  xnorm=xnorm+x(j)**2
+      xnorm=sqrt(xnorm)
+c			      successful return.
+  110 mode=1
+      return
+c			      error return.	  n .le. 0.
+  120 mode=2
+      return
+c			      returning with constraints not compatible.
+  130 mode=4
+      return
+      end
diff --git a/math/llsq/original_f/nnls.f b/math/llsq/original_f/nnls.f
new file mode 100644
index 00000000..fb5aa767
--- /dev/null
+++ b/math/llsq/original_f/nnls.f
@@ -0,0 +1,278 @@
+c     subroutine nnls  (a,mda,m,n,b,x,rnorm,w,zz,index,mode)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 june 15
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c
+c	  **********   nonnegative least squares   **********
+c
+c     given an m by n matrix, a, and an m-vector, b,  compute an
+c     n-vector, x, which solves the least squares problem
+c
+c		       a * x = b  subject to x .ge. 0
+c
+c     a(),mda,m,n     mda is the first dimensioning parameter for the
+c		      array, a().   on entry a() contains the m by n
+c		      matrix, a.	   on exit a() contains
+c		      the product matrix, q*a , where q is an
+c		      m by m orthogonal matrix generated implicitly by
+c		      this subroutine.
+c     b()     on entry b() contains the m-vector, b.   on exit b() con-
+c	      tains q*b.
+c     x()     on entry x() need not be initialized.  on exit x() will
+c	      contain the solution vector.
+c     rnorm   on exit rnorm contains the euclidean norm of the
+c	      residual vector.
+c     w()     an n-array of working space.  on exit w() will contain
+c	      the dual solution vector.   w will satisfy w(i) = 0.
+c	      for all i in set p  and w(i) .le. 0. for all i in set z
+c     zz()     an m-array of working space.
+c     index()	  an integer working array of length at least n.
+c		  on exit the contents of this array define the sets
+c		  p and z as follows..
+c
+c		  index(1)   thru index(nsetp) = set p.
+c		  index(iz1) thru index(iz2)   = set z.
+c		  iz1 = nsetp + 1 = npp1
+c		  iz2 = n
+c     mode    this is a success-failure flag with the following
+c	      meanings.
+c	      1     the solution has been computed successfully.
+c	      2     the dimensions of the problem are bad.
+c		    either m .le. 0 or n .le. 0.
+c	      3    iteration count exceeded.  more than 3*n iterations.
+c
+      subroutine nnls (a,mda,m,n,b,x,rnorm,w,zz,index,mode)
+      dimension a(mda,n), b(m), x(n), w(n), zz(m)
+      integer index(n)
+      zero=0.
+      one=1.
+      two=2.
+      factor=0.01
+c
+      mode=1
+      if (m.gt.0.and.n.gt.0) go to 10
+      mode=2
+      return
+   10 iter=0
+      itmax=3*n
+c
+c		     initialize the arrays index() and x().
+c
+	  do 20 i=1,n
+	  x(i)=zero
+   20	  index(i)=i
+c
+      iz2=n
+      iz1=1
+      nsetp=0
+      npp1=1
+c			      ******  main loop begins here  ******
+   30 continue
+c		   quit if all coefficients are already in the solution.
+c			 or if m cols of a have been triangularized.
+c
+      if (iz1.gt.iz2.or.nsetp.ge.m) go to 350
+c
+c	  compute components of the dual (negative gradient) vector w().
+c
+	  do 50 iz=iz1,iz2
+	  j=index(iz)
+	  sm=zero
+	      do 40 l=npp1,m
+   40	      sm=sm+a(l,j)*b(l)
+   50	  w(j)=sm
+c				    find largest positive w(j).
+   60 wmax=zero
+	  do 70 iz=iz1,iz2
+	  j=index(iz)
+	  if (w(j).le.wmax) go to 70
+	  wmax=w(j)
+	  izmax=iz
+   70	  continue
+c
+c	      if wmax .le. 0. go to termination.
+c	      this indicates satisfaction of the kuhn-tucker conditions.
+c
+      if (wmax) 350,350,80
+   80 iz=izmax
+      j=index(iz)
+c
+c     the sign of w(j) is ok for j to be moved to set p.
+c     begin the transformation and check new diagonal element to avoid
+c     near linear dependence.
+c
+      asave=a(npp1,j)
+      call h12 (1,npp1,npp1+1,m,a(1,j),1,up,dummy,1,1,0)
+      unorm=zero
+      if (nsetp.eq.0) go to 100
+	  do 90 l=1,nsetp
+   90	  unorm=unorm+a(l,j)**2
+  100 unorm=sqrt(unorm)
+      if (diff(unorm+abs(a(npp1,j))*factor,unorm)) 130,130,110
+c
+c     col j is sufficiently independent.  copy b into zz, update zz and
+c   > solve for ztest ( = proposed new value for x(j) ).
+c
+  110	  do 120 l=1,m
+  120	  zz(l)=b(l)
+      call h12 (2,npp1,npp1+1,m,a(1,j),1,up,zz,1,1,1)
+      ztest=zz(npp1)/a(npp1,j)
+c
+c				      see if ztest is positive
+c     reject j as a candidate to be moved from set z to set p.
+c     restore a(npp1,j), set w(j)=0., and loop back to test dual
+c
+      if (ztest) 130,130,140
+c
+c     coeffs again.
+c
+  130 a(npp1,j)=asave
+      w(j)=zero
+      go to 60
+c
+c     the index  j=index(iz)  has been selected to be moved from
+c     set z to set p.	 update b,  update indices,  apply householder
+c     transformations to cols in new set z,  zero subdiagonal elts in
+c     col j,  set w(j)=0.
+c
+  140	  do 150 l=1,m
+  150	  b(l)=zz(l)
+c
+      index(iz)=index(iz1)
+      index(iz1)=j
+      iz1=iz1+1
+      nsetp=npp1
+      npp1=npp1+1
+c
+      if (iz1.gt.iz2) go to 170
+	  do 160 jz=iz1,iz2
+	  jj=index(jz)
+  160	  call h12 (2,nsetp,npp1,m,a(1,j),1,up,a(1,jj),1,mda,1)
+  170 continue
+c
+      if (nsetp.eq.m) go to 190
+	  do 180 l=npp1,m
+  180	  a(l,j)=zero
+  190 continue
+c
+      w(j)=zero
+c				 solve the triangular system.
+c				 store the solution temporarily in zz().
+      assign 200 to next
+      go to 400
+  200 continue
+c
+c			******	secondary loop begins here ******
+c
+c			   iteration counter.
+c
+  210 iter=iter+1
+      if (iter.le.itmax) go to 220
+      mode=3
+      write (6,440)
+      go to 350
+  220 continue
+c
+c		     see if all new constrained coeffs are feasible.
+c				   if not compute alpha.
+c
+      alpha=two
+	  do 240 ip=1,nsetp
+	  l=index(ip)
+	  if (zz(ip)) 230,230,240
+c
+  230	  t=-x(l)/(zz(ip)-x(l))
+	  if (alpha.le.t) go to 240
+	  alpha=t
+	  jj=ip
+  240	  continue
+c
+c	   if all new constrained coeffs are feasible then alpha will
+c	   still = 2.	 if so exit from secondary loop to main loop.
+c
+      if (alpha.eq.two) go to 330
+c
+c	   otherwise use alpha which will be between 0. and 1. to
+c	   interpolate between the old x and the new zz.
+c
+	  do 250 ip=1,nsetp
+	  l=index(ip)
+  250	  x(l)=x(l)+alpha*(zz(ip)-x(l))
+c
+c	 modify a and b and the index arrays to move coefficient i
+c	 from set p to set z.
+c
+      i=index(jj)
+  260 x(i)=zero
+c
+      if (jj.eq.nsetp) go to 290
+      jj=jj+1
+	  do 280 j=jj,nsetp
+	  ii=index(j)
+	  index(j-1)=ii
+	  call g1 (a(j-1,ii),a(j,ii),cc,ss,a(j-1,ii))
+	  a(j,ii)=zero
+	      do 270 l=1,n
+	      if (l.ne.ii) call g2 (cc,ss,a(j-1,l),a(j,l))
+  270	      continue
+  280	  call g2 (cc,ss,b(j-1),b(j))
+  290 npp1=nsetp
+      nsetp=nsetp-1
+      iz1=iz1-1
+      index(iz1)=i
+c
+c	 see if the remaining coeffs in set p are feasible.  they should
+c	 be because of the way alpha was determined.
+c	 if any are infeasible it is due to round-off error.  any
+c	 that are nonpositive will be set to zero
+c	 and moved from set p to set z.
+c
+	  do 300 jj=1,nsetp
+	  i=index(jj)
+	  if (x(i)) 260,260,300
+  300	  continue
+c
+c	  copy b( ) into zz( ).  then solve again and loop back.
+c
+
+	  do 310 i=1,m
+  310	  zz(i)=b(i)
+      assign 320 to next
+      go to 400
+  320 continue
+      go to 210
+c		       ******  end of secondary loop  ******
+c
+  330	  do 340 ip=1,nsetp
+	  i=index(ip)
+  340	  x(i)=zz(ip)
+c	 all new coeffs are positive.  loop back to beginning.
+      go to 30
+c
+c			 ******  end of main loop  ******
+c
+c			 come to here for termination.
+c		      compute the norm of the final residual vector.
+c
+  350 sm=zero
+      if (npp1.gt.m) go to 370
+	  do 360 i=npp1,m
+  360	  sm=sm+b(i)**2
+      go to 390
+  370	  do 380 j=1,n
+  380	  w(j)=zero
+  390 rnorm=sqrt(sm)
+      return
+c
+c     the following block of code is used as an internal subroutine
+c     to solve the triangular system, putting the solution in zz().
+c
+  400	  do 430 l=1,nsetp
+	  ip=nsetp+1-l
+	  if (l.eq.1) go to 420
+	      do 410 ii=1,ip
+  410	      zz(ii)=zz(ii)-a(ii,jj)*zz(ip+1)
+  420	  jj=index(ip)
+  430	  zz(ip)=zz(ip)/a(ip,jj)
+      go to next, (200,320)
+  440 format (35h0 nnls quitting on iteration count.)
+      end
diff --git a/math/llsq/original_f/qrbd.f b/math/llsq/original_f/qrbd.f
new file mode 100644
index 00000000..d9a12a24
--- /dev/null
+++ b/math/llsq/original_f/qrbd.f
@@ -0,0 +1,208 @@
+c     subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   qr algorithm for singular values of a bidiagonal matrix.
+c
+c     the bidiagonal matrix
+c
+c			(q1,e2,0...    )
+c			(   q2,e3,0... )
+c		 d=	(	.      )
+c			(	  .   0)
+c			(	    .en)
+c			(	   0,qn)
+c
+c		  is pre and post multiplied by
+c		  elementary rotation matrices
+c		  ri and pi so that
+c
+c		  rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sn)
+c
+c		  to within working accuracy.
+c
+c  1. ei and qi occupy e(i) and q(i) as input.
+c
+c  2. rm...r1*c replaces 'c' in storage as output.
+c
+c  3. v*p1**(t)...pm**(t) replaces 'v' in storage as output.
+c
+c  4. si occupies q(i) as output.
+c
+c  5. the si's are nonincreasing and nonnegative.
+c
+c     this code is based on the paper and 'algol' code..
+c ref..
+c  1. reinsch,c.h. and golub,g.h. 'singular value decomposition
+c     and least squares solutions' (numer. math.), vol. 14,(1970).
+c
+      subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc)
+      logical wntv ,havers,fail
+      dimension q(nn),e(nn),v(mdv,nn),c(mdc,ncc)
+      zero=0.
+      one=1.
+      two=2.
+c
+      n=nn
+      ipass=1
+      if (n.le.0) return
+      n10=10*n
+      wntv=nrv.gt.0
+      havers=ncc.gt.0
+      fail=.false.
+      nqrs=0
+      e(1)=zero
+      dnorm=zero
+	   do 10 j=1,n
+   10	   dnorm=amax1(abs(q(j))+abs(e(j)),dnorm)
+	   do 200 kk=1,n
+	   k=n+1-kk
+c
+c     test for splitting or rank deficiencies..
+c	  first make test for last diagonal term, q(k), being small.
+   20	    if(k.eq.1) go to 50
+	    if(diff(dnorm+q(k),dnorm)) 50,25,50
+c
+c     since q(k) is small we will make a special pass to
+c     transform e(k) to zero.
+c
+   25	   cs=zero
+	   sn=-one
+		do 40 ii=2,k
+		i=k+1-ii
+		f=-sn*e(i+1)
+		e(i+1)=cs*e(i+1)
+		call g1 (q(i),f,cs,sn,q(i))
+c	  transformation constructed to zero position (i,k).
+c
+		if (.not.wntv) go to 40
+		     do 30 j=1,nrv
+   30		     call g2 (cs,sn,v(j,i),v(j,k))
+c	       accumulate rt. transformations in v.
+c
+   40		continue
+c
+c	  the matrix is now bidiagonal, and of lower order
+c	  since e(k) .eq. zero..
+c
+   50		do 60 ll=1,k
+		l=k+1-ll
+		if(diff(dnorm+e(l),dnorm)) 55,100,55
+   55		if(diff(dnorm+q(l-1),dnorm)) 60,70,60
+   60		continue
+c     this loop can't complete since e(1) = zero.
+c
+	   go to 100
+c
+c	  cancellation of e(l), l.gt.1.
+   70	   cs=zero
+	   sn=-one
+		do 90 i=l,k
+		f=-sn*e(i)
+		e(i)=cs*e(i)
+		if(diff(dnorm+f,dnorm)) 75,100,75
+   75		call g1 (q(i),f,cs,sn,q(i))
+		if (.not.havers) go to 90
+		     do 80 j=1,ncc
+   80		     call g2 (cs,sn,c(i,j),c(l-1,j))
+   90		continue
+c
+c	  test for convergence..
+  100	   z=q(k)
+	   if (l.eq.k) go to 170
+c
+c	  shift from bottom 2 by 2 minor of b**(t)*b.
+	   x=q(l)
+	   y=q(k-1)
+	   g=e(k-1)
+	   h=e(k)
+	   f=((y-z)*(y+z)+(g-h)*(g+h))/(two*h*y)
+	   g=sqrt(one+f**2)
+	   if (f.lt.zero) go to 110
+	   t=f+g
+	   go to 120
+  110	   t=f-g
+  120	   f=((x-z)*(x+z)+h*(y/t-h))/x
+c
+c	  next qr sweep..
+	   cs=one
+	   sn=one
+	   lp1=l+1
+		do 160 i=lp1,k
+		g=e(i)
+		y=q(i)
+		h=sn*g
+		g=cs*g
+		call g1 (f,h,cs,sn,e(i-1))
+		f=x*cs+g*sn
+		g=-x*sn+g*cs
+		h=y*sn
+		y=y*cs
+		if (.not.wntv) go to 140
+c
+c	       accumulate rotations (from the right) in 'v'
+		     do 130 j=1,nrv
+  130		     call g2 (cs,sn,v(j,i-1),v(j,i))
+  140		call g1 (f,h,cs,sn,q(i-1))
+		f=cs*g+sn*y
+		x=-sn*g+cs*y
+		if (.not.havers) go to 160
+		     do 150 j=1,ncc
+  150		     call g2 (cs,sn,c(i-1,j),c(i,j))
+c	       apply rotations from the left to
+c	       right hand sides in 'c'..
+c
+  160		continue
+	   e(l)=zero
+	   e(k)=f
+	   q(k)=x
+	   nqrs=nqrs+1
+	   if (nqrs.le.n10) go to 20
+c	   return to 'test for splitting'.
+c
+	   fail=.true.
+c     ..
+c     cutoff for convergence failure. 'nqrs' will be 2*n usually.
+  170	   if (z.ge.zero) go to 190
+	   q(k)=-z
+	   if (.not.wntv) go to 190
+		do 180 j=1,nrv
+  180		v(j,k)=-v(j,k)
+  190	   continue
+c	  convergence. q(k) is made nonnegative..
+c
+  200	   continue
+      if (n.eq.1) return
+	   do 210 i=2,n
+	   if (q(i).gt.q(i-1)) go to 220
+  210	   continue
+      if (fail) ipass=2
+      return
+c     ..
+c     every singular value is in order..
+  220	   do 270 i=2,n
+	   t=q(i-1)
+	   k=i-1
+		do 230 j=i,n
+		if (t.ge.q(j)) go to 230
+		t=q(j)
+		k=j
+  230		continue
+	   if (k.eq.i-1) go to 270
+	   q(k)=q(i-1)
+	   q(i-1)=t
+	   if (.not.havers) go to 250
+		do 240 j=1,ncc
+		t=c(i-1,j)
+		c(i-1,j)=c(k,j)
+  240		c(k,j)=t
+  250	   if (.not.wntv) go to 270
+		do 260 j=1,nrv
+		t=v(j,i-1)
+		v(j,i-1)=v(j,k)
+  260		v(j,k)=t
+  270	   continue
+c	  end of ordering algorithm.
+c
+      if (fail) ipass=2
+      return
+      end
diff --git a/math/llsq/original_f/sfeout.f b/math/llsq/original_f/sfeout.f
new file mode 100644
index 00000000..fb439e40
--- /dev/null
+++ b/math/llsq/original_f/sfeout.f
@@ -0,0 +1,64 @@
+      subroutine mfeout (a,mda,m,n,names,mode)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   subroutine for matrix output with labeling.
+c
+c     a( )	   matrix to be output
+c		   mda	   first dimension of a array
+c		   m	     no. of rows in a matrix
+c		   n	     no. of cols in a matrix
+c     names()	   array of names.  if names(1) = 1h , the rest
+c		   of the names() array will be ignored.
+c     mode	   =1	for   4p8f15.0	format	for v matrix.
+c		   =2	for   8e15.8  format  for candidate solutions.
+c
+      dimension    a(mda,01)
+      integer names(m),ihead(2)
+      logical	notblk
+      data  maxcol/8/, iblank/1h /,ihead(1)/4h col/,ihead(2)/4hsoln/
+c
+      notblk=names(1).ne.iblank
+      if (m.le.0.or.n.le.0) return
+c
+      if (mode.eq.2) go to 10
+      write (6,70)
+      go to 20
+   10 write (6,80)
+   20 continue
+c
+      nblock=n/maxcol
+      last=n-nblock*maxcol
+      ncol=maxcol
+      j1=1
+c
+c			     main loop starts here
+c
+   30 if (nblock.gt.0) go to 40
+      if (last.le.0) return
+      ncol=last
+      last=0
+c
+   40 j2=j1+ncol-1
+      write (6,90) (ihead(mode),j,j=j1,j2)
+c
+	   do 60 i=1,m
+	   name=iblank
+	   if (notblk) name=names(i)
+c
+	   if (mode.eq.2) go to 50
+	   write (6,100) i,name,(a(i,j),j=j1,j2)
+	   go to 60
+   50	   write (6,110) i,name,(a(i,j),j=j1,j2)
+   60	   continue
+c
+      j1=j1+maxcol
+      nblock=nblock-1
+      go to 30
+c
+   70 format (45h0v-matrix of the singular value decomposition,
+     * 8h of a*d./47h (elements of v scaled up by a factor of 10**4))
+   80 format (35h0sequence of candidate solutions, x)
+   90 format (1h0,11x,8(6x,a4,i4,1x)/1x)
+  100 format (1x,i3,1x,a6,1x,4p8f15.0)
+  110 format (1x,i3,1x,a6,1x,8e15.8)
+      end
diff --git a/math/llsq/original_f/sva.f b/math/llsq/original_f/sva.f
new file mode 100644
index 00000000..469a90ef
--- /dev/null
+++ b/math/llsq/original_f/sva.f
@@ -0,0 +1,193 @@
+      subroutine sva (a,mda,m,n,mdata,b,sing,names,iscale,d)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c		singular value analysis printout
+c
+c     iscale	   set by user to 1, 2, or 3 to select column scaling
+c		   option.
+c		   1   subr will use identity scaling and ignore the d()
+c		       array.
+c		   2   subr will scale nonzero cols to have unit euclide
+c		       length and will store reciprocal lengths of
+c		       original nonzero cols in d().
+c		   3   user supplies col scale factors in d(). subr
+c		       will mult col j by d(j) and remove the scaling
+c		       from the soln at the end.
+c
+      dimension  a(mda,n), b(m), sing(01),d(n)
+c     sing(3*n)
+      dimension names(n)
+      logical test
+      double precision	 sb, dzero
+      dzero=0.d0
+      one=1.
+      zero=0.
+      if (m.le.0 .or. n.le.0) return
+      np1=n+1
+      write (6,270)
+      write (6,260) m,n,mdata
+      go to (60,10,10), iscale
+c
+c     apply column scaling to a
+c
+   10	   do 50 j=1,n
+	   a1=d(j)
+	   go to (20,20,40), iscale
+   20	   sb=dzero
+		do 30 i=1,m
+   30		sb=sb+dble(a(i,j))**2
+	   a1=dsqrt(sb)
+	   if (a1.eq.zero) a1=one
+	   a1=one/a1
+	   d(j)=a1
+   40		do 50 i=1,m
+   50		a(i,j)=a(i,j)*a1
+      write (6,280) iscale,(d(j),j=1,n)
+      go to 70
+   60 continue
+      write (6,290)
+   70 continue
+c
+c     obtain  sing. value decomp. of scaled matrix.
+c
+c**********************************************************
+      call svdrs (a,mda,m,n,b,1,1,sing)
+c**********************************************************
+c
+c     print the v matrix.
+c
+      call mfeout (a,mda,n,n,names,1)
+      if (iscale.eq.1) go to 90
+c
+c     replace v by d*v in the array a()
+c
+	   do 80 i=1,n
+		do 80 j=1,n
+   80		a(i,j)=d(i)*a(i,j)
+c
+c     g  now in  b array.  v now in a array.
+c
+c
+c     obtain summary output.
+c
+   90 continue
+      write (6,220)
+c
+c     compute cumulative sums of squares of components of
+c     g and store them in sing(i), i=minmn+1,...,2*minmn+1
+c
+      sb=dzero
+      minmn=min0(m,n)
+      minmn1=minmn+1
+      if (m.eq.minmn) go to 110
+	   do 100 i=minmn1,m
+  100	   sb=sb+dble(b(i))**2
+  110 sing(2*minmn+1)=sb
+	   do 120 jj=1,minmn
+	   j=minmn+1-jj
+	   sb=sb+dble(b(j))**2
+	   js=minmn+j
+  120	   sing(js)=sb
+      a3=sing(minmn+1)
+      a4=sqrt(a3/float(max0(1,mdata)))
+      write (6,230) a3,a4
+c
+      nsol=0
+c
+c
+c
+	   do 160 k=1,minmn
+	   if (sing(k).eq.zero) go to 130
+	   nsol=k
+	   pi=b(k)/sing(k)
+	   a1=one/sing(k)
+	   a2=b(k)**2
+	   a3=sing(minmn1+k)
+	   a4=sqrt(a3/float(max0(1,mdata-k)))
+	   test=sing(k).ge.100..or.sing(k).lt..001
+	   if (test) write (6,240) k,sing(k),pi,a1,b(k),a2,a3,a4
+	   if (.not.test) write (6,250) k,sing(k),pi,a1,b(k),a2,a3,a4
+	   go to 140
+  130	   write (6,240) k,sing(k)
+	   pi=zero
+  140		do 150 i=1,n
+		a(i,k)=a(i,k)*pi
+  150		if (k.gt.1) a(i,k)=a(i,k)+a(i,k-1)
+  160	   continue
+c
+c     compute and print values of ynorm and rnorm.
+c
+      write (6,300)
+      j=0
+      ysq=zero
+      go to 180
+  170 j=j+1
+      ysq=ysq+(b(j)/sing(j))**2
+  180 ynorm=sqrt(ysq)
+      js=minmn1+j
+      rnorm=sqrt(sing(js))
+      yl=-1000.
+      if (ynorm.gt.0.) yl=alog10(ynorm)
+      rl=-1000.
+      if (rnorm.gt.0.) rl=alog10(rnorm)
+      write (6,310) j,ynorm,rnorm,yl,rl
+      if (j.lt.nsol) go to 170
+c
+c     compute values of xnorm and rnorm for a sequence of values of
+c     the levenberg-marquardt parameter.
+c
+      if (sing(1).eq.zero) go to 210
+      el=alog10(sing(1))+one
+      el2=alog10(sing(nsol))-one
+      del=(el2-el)/20.
+      ten=10.
+      aln10=alog(ten)
+      write (6,320)
+	   do 200 ie=1,21
+c     compute	     alamb=10.**el
+	   alamb=exp(aln10*el)
+	   ys=0.
+	   js=minmn1+nsol
+	   rs=sing(js)
+		do 190 i=1,minmn
+		sl=sing(i)**2+alamb**2
+		ys=ys+(b(i)*sing(i)/sl)**2
+		rs=rs+(b(i)*(alamb**2)/sl)**2
+  190		continue
+	   ynorm=sqrt(ys)
+	   rnorm=sqrt(rs)
+	   rl=-1000.
+	   if (rnorm.gt.zero) rl=alog10(rnorm)
+	   yl=-1000.
+	   if (ynorm.gt.zero) yl=alog10(ynorm)
+	   write (6,330) alamb,ynorm,rnorm,el,yl,rl
+	   el=el+del
+  200	   continue
+c
+c     print candidate solutions.
+c
+  210 if (nsol.ge.1) call mfeout (a,mda,n,nsol,names,2)
+      return
+  220 format (42h0 index  sing. value       p coef         ,48hrecip. s.
+     1v.             g coef            g**2  ,39h            c.s.s.
+     2    n.s.r.c.s.s.)
+  230 format (1h ,5x,1h0,88x,1pe15.4,1pe17.4)
+  240 format (1h ,i6,e12.4,1p(e15.4,4x,e15.4,4x,e15.4,4x,e15.4,4x,e15.4,
+     12x,e15.4))
+  250 format (1h ,i6,f12.4,1p(e15.4,4x,e15.4,4x,e15.4,4x,e15.4,4x,e15.4,
+     12x,e15.4))
+  260 format (5h0m = ,i6,8h,   n = ,i4,12h,   mdata = ,i8)
+  270 format (45h0singular value analysis of the least squares,42h probl
+     1em,  a*x=b,  scaled as  (a*d)*y=b  .)
+  280 format (19h0scaling option no.,i2,18h.  d is a diagonal,46h matrix
+     1 with the following diagonal elements../(5x,10e12.4))
+  290 format (50h0scaling option no. 1.   d is the identity matrix./1x)
+  300 format (6h0index,13x,28h         ynorm         rnorm,14x,28h  log1
+     10(ynorm)  log10(rnorm)/1x)
+  310 format (1h ,i4,14x,2e14.5,14x,2f14.5)
+  320 format (54h0norms of solution and residual vectors for a range of,
+     144h values of the levenberg-marquardt parameter,9h, lambda./1h0,4x
+     2,42h        lambda         ynorm         rnorm,42h log10(lambda)
+     3log10(ynorm)  log10(rnorm))
+  330 format (5x,3e14.5,3f14.5)
+      end
diff --git a/math/llsq/original_f/svdrs.f b/math/llsq/original_f/svdrs.f
new file mode 100644
index 00000000..66e488dc
--- /dev/null
+++ b/math/llsq/original_f/svdrs.f
@@ -0,0 +1,205 @@
+c     subroutine svdrs (a,mda,mm,nn,b,mdb,nb,s)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1974 mar 1
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   singular value decomposition also treating right side vector.
+c
+c     the array s occupies 3*n cells.
+c     a occupies m*n cells
+c     b occupies m*nb cells.
+c
+c     special singular value decomposition subroutine.
+c     we have the m x n matrix a and the system a*x=b to solve.
+c     either m .ge. n  or  m .lt. n is permitted.
+c			the singular value decomposition
+c     a = u*s*v**(t) is made in such a way that one gets
+c	 (1) the matrix v in the first n rows and columns of a.
+c	 (2) the diagonal matrix of ordered singular values in
+c	     the first n cells of the array s(ip), ip .ge. 3*n.
+c	 (3) the matrix product u**(t)*b=g gets placed back in b.
+c	 (4) the user must complete the solution and do his own
+c	     singular value analysis.
+c     *******
+c     give special
+c     treatment to rows and columns which are entirely zero.  this
+c     causes certain zero sing. vals. to appear as exact zeros rather
+c     than as about eta times the largest sing. val.   it similarly
+c     cleans up the associated columns of u and v.
+c     method..
+c      1. exchange cols of a to pack nonzero cols to the left.
+c	  set n = no. of nonzero cols.
+c	  use locations a(1,nn),a(1,nn-1),...,a(1,n+1) to record the
+c	  col permutations.
+c      2. exchange rows of a to pack nonzero rows to the top.
+c	  quit packing if find n nonzero rows.	make same row exchanges
+c	  in b.  set m so that all nonzero rows of the permuted a
+c	  are in first m rows.	if m .le. n then all m rows are
+c	  nonzero.  if m .gt. n then	  the first n rows are known
+c	  to be nonzero,and rows n+1 thru m may be zero or nonzero.
+c      3. apply original algorithm to the m by n problem.
+c      4. move permutation record from a(,) to s(i),i=n+1,...,nn.
+c      5. build v up from  n by n  to  nn by nn  by placing ones on
+c	  the diagonal and zeros elsewhere.  this is only partly done
+c	  explicitly.  it is completed during step 6.
+c      6. exchange rows of v to compensate for col exchanges of step 2.
+c      7. place zeros in  s(i),i=n+1,nn  to represent zero sing vals.
+c
+      subroutine svdrs (a,mda,mm,nn,b,mdb,nb,s)
+      dimension a(mda,nn),b(mdb,nb),s(nn,3)
+      zero=0.
+      one=1.
+c
+c			      begin.. special for zero rows and cols.
+c
+c			      pack the nonzero cols to the left
+c
+      n=nn
+      if (n.le.0.or.mm.le.0) return
+      j=n
+   10 continue
+	  do 20 i=1,mm
+	  if (a(i,j)) 50,20,50
+   20	  continue
+c
+c	  col j  is zero. exchange it with col n.
+c
+      if (j.eq.n) go to 40
+	  do 30 i=1,mm
+   30	  a(i,j)=a(i,n)
+   40 continue
+      a(1,n)=j
+      n=n-1
+   50 continue
+      j=j-1
+      if (j.ge.1) go to 10
+c			      if n=0 then a is entirely zero and svd
+c			      computation can be skipped
+      ns=0
+      if (n.eq.0) go to 240
+c			      pack nonzero rows to the top
+c			      quit packing if find n nonzero rows
+      i=1
+      m=mm
+   60 if (i.gt.n.or.i.ge.m) go to 150
+      if (a(i,i)) 90,70,90
+   70	  do 80 j=1,n
+	  if (a(i,j)) 90,80,90
+   80	  continue
+      go to 100
+   90 i=i+1
+      go to 60
+c			      row i is zero
+c			      exchange rows i and m
+  100 if(nb.le.0) go to 115
+	  do 110 j=1,nb
+	  t=b(i,j)
+	  b(i,j)=b(m,j)
+  110	  b(m,j)=t
+  115	  do 120 j=1,n
+  120	  a(i,j)=a(m,j)
+      if (m.gt.n) go to 140
+	  do 130 j=1,n
+  130	  a(m,j)=zero
+  140 continue
+c			      exchange is finished
+      m=m-1
+      go to 60
+c
+  150 continue
+c			      end.. special for zero rows and columns
+c			      begin.. svd algorithm
+c     method..
+c     (1)     reduce the matrix to upper bidiagonal form with
+c     householder transformations.
+c	   h(n)...h(1)aq(1)...q(n-2) = (d**t,0)**t
+c     where d is upper bidiagonal.
+c
+c     (2)     apply h(n)...h(1) to b.  here h(n)...h(1)*b replaces b
+c     in storage.
+c
+c     (3)     the matrix product w= q(1)...q(n-2) overwrites the first
+c     n rows of a in storage.
+c
+c     (4)     an svd for d is computed.  here k rotations ri and pi are
+c     computed so that
+c	   rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sm)
+c     to working accuracy.  the si are nonnegative and nonincreasing.
+c     here rk...r1*b overwrites b in storage while
+c     a*p1**(t)...pk**(t)  overwrites a in storage.
+c
+c     (5)     it follows that,with the proper definitions,
+c     u**(t)*b overwrites b, while v overwrites the first n row and
+c     columns of a.
+c
+      l=min0(m,n)
+c	      the following loop reduces a to upper bidiagonal and
+c	      also applies the premultiplying transformations to b.
+c
+	  do 170 j=1,l
+	  if (j.ge.m) go to 160
+	  call h12 (1,j,j+1,m,a(1,j),1,t,a(1,j+1),1,mda,n-j)
+	  call h12 (2,j,j+1,m,a(1,j),1,t,b,1,mdb,nb)
+  160	  if (j.ge.n-1) go to 170
+	  call h12 (1,j+1,j+2,n,a(j,1),mda,s(j,3),a(j+1,1),mda,1,m-j)
+  170	  continue
+c
+c     copy the bidiagonal matrix into the array s() for qrbd.
+c
+      if (n.eq.1) go to 190
+	  do 180 j=2,n
+	  s(j,1)=a(j,j)
+  180	  s(j,2)=a(j-1,j)
+  190 s(1,1)=a(1,1)
+c
+      ns=n
+      if (m.ge.n) go to 200
+      ns=m+1
+      s(ns,1)=zero
+      s(ns,2)=a(m,m+1)
+  200 continue
+c
+c     construct the explicit n by n product matrix, w=q1*q2*...*ql*i
+c     in the array a().
+c
+	  do 230 k=1,n
+	  i=n+1-k
+	  if(i.gt.min0(m,n-2)) go to 210
+	  call h12 (2,i+1,i+2,n,a(i,1),mda,s(i,3),a(1,i+1),1,mda,n-i)
+  210	      do 220 j=1,n
+  220	      a(i,j)=zero
+  230	  a(i,i)=one
+c
+c	   compute the svd of the bidiagonal matrix
+c
+      call qrbd (ipass,s(1,1),s(1,2),ns,a,mda,n,b,mdb,nb)
+c
+      go to (240,310), ipass
+  240 continue
+      if (ns.ge.n) go to 260
+      nsp1=ns+1
+	  do 250 j=nsp1,n
+  250	  s(j,1)=zero
+  260 continue
+      if (n.eq.nn) return
+      np1=n+1
+c				   move record of permutations
+c				   and store zeros
+	  do 280 j=np1,nn
+	  s(j,1)=a(1,j)
+	      do 270 i=1,n
+  270	      a(i,j)=zero
+  280	  continue
+c			      permute rows and set zero singular values.
+	  do 300 k=np1,nn
+	  i=s(k,1)
+	  s(k,1)=zero
+	      do 290 j=1,nn
+	      a(k,j)=a(i,j)
+  290	      a(i,j)=zero
+	  a(i,k)=one
+  300	  continue
+c			      end.. special for zero rows and columns
+      return
+  310 write (6,320)
+      stop
+  320 format (49h convergence failure in qr bidiagonal svd routine)
+      end
diff --git a/math/llsq/progs/README b/math/llsq/progs/README
new file mode 100644
index 00000000..e0c5a142
--- /dev/null
+++ b/math/llsq/progs/README
@@ -0,0 +1,5 @@
+
+This directory contains both Fortran and IRAF preprocessor programs
+demonstrating the use of the Lawson-Hanson llsq procedures.  The Fortran
+programs have not been modified to reflect the removal of Fortran i/o
+from the library procedures as installed in the IRAF llsq.a library.
diff --git a/math/llsq/progs/band.x b/math/llsq/progs/band.x
new file mode 100644
index 00000000..f65290ba
--- /dev/null
+++ b/math/llsq/progs/band.x
@@ -0,0 +1,70 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+
+define	MAXPTS		1024
+
+task	band
+
+# This procedure solves for the natural cubic spline which interpolates
+# the curve y[i] = 100., for i = 1 to n.  The matrix is of length n+2,
+# and has a bandwidth of 3.  Lawsons and Hansons routines BNDACC and 
+# BNDSOL (modified to return an error code IER) are used to solve the
+# system.  The computations are carried out in single precision.
+
+
+procedure band()
+
+real	g[MAXPTS,4], c[MAXPTS], rnorm
+int	ip, ir, i, n, ier, geti()
+real	marktime, cptime()
+
+begin
+	n = min (MAXPTS, geti ("npts"))
+	marktime = cptime()
+
+	ip = 1
+	ir = 1
+
+	g[ir,1] = 6.			# first row
+	g[ir,2] = -12.
+	g[ir,3] = 6.
+	g[ir,4] = 0.
+	call bndacc (g, MAXPTS, 3, ip, ir, 1, 1)
+
+	do i = 2, n-1 {			# tridiagonal elements
+	    g[ir,1] = 1.
+	    g[ir,2] = 4.
+	    g[ir,3] = 1.
+	    g[ir,4] = 100.
+	    call bndacc (g, MAXPTS, 3, ip, ir, 1, i-1)
+	}
+
+	g[ir,1] = 6.			# last row
+	g[ir,2] = -12.
+	g[ir,3] = 6.
+	g[ir,4] = 0.
+	call bndacc (g, MAXPTS, 3, ip, ir, 1, n-2)
+
+	call printf ("matrix accumulation took %8.2f cpu seconds\n")
+	    call pargr (cptime() - marktime)
+	marktime = cptime()
+
+					# solve for the b-spline coeff C
+	call bndsol (1, g, MAXPTS, 3, ip, ir, c, n, rnorm, ier)
+
+	call printf ("solution took %8.2f cpu seconds (rnorm=%g, ier=%d)\n")
+	    call pargr (cptime() - marktime)
+	    call pargr (rnorm)
+	    call pargi (ier)
+
+	call printf ("selected coefficients:\n")
+	for (i=1;  i <= n;) {		# print the first and last 4 coeff
+	    call printf ("%8d%15.6f\n")
+		call pargi (i)
+		call pargr (c[i])
+	    if (i == 4)
+		i = max(i+1, n-3)
+	    else
+		i = i + 1
+	}
+end
diff --git a/math/llsq/progs/data4 b/math/llsq/progs/data4
new file mode 100644
index 00000000..33bf1b08
--- /dev/null
+++ b/math/llsq/progs/data4
@@ -0,0 +1,15 @@
+  -.13405547  -.20162827  -.16930778  -.18971990  -.17387234	  -.4361
+  -.10379475  -.15766336  -.13346256  -.14848550  -.13597690	  -.3437
+  -.08779597  -.12883867  -.10683007  -.12011796  -.10932972	  -.2657
+   .02058554   .00335331  -.01641270   .00078606   .00271659	  -.0392
+  -.03248093  -.01876799   .00410639  -.01405894  -.01384391	   .0193
+   .05967662   .06667714   .04352153   .05740438   .05024962	   .0747
+   .06712457   .07352437   .04489770   .06471862   .05876455	   .0935
+   .08687186   .09368296   .05672327   .08141043   .07302320	   .1079
+   .02149662   .06222662   .07213486   .06200069   .05570931	   .1930
+   .06687407   .10344506   .09153849   .09508223   .08393667	   .2058
+   .15879069   .18088339   .11540692   .16160727   .14796479	   .2606
+   .17642887   .20361830   .13057860   .18385729   .17005549	   .3142
+   .11414080   .17259611   .14816471   .16007466   .14374096	   .3529
+   .07846038   .14669563   .14365800   .14003842   .12571177	   .3615
+   .10803175   .16994623   .14971519   .15885312   .14301547	   .3647
diff --git a/math/llsq/progs/lsq.x b/math/llsq/progs/lsq.x
new file mode 100644
index 00000000..ae002a16
--- /dev/null
+++ b/math/llsq/progs/lsq.x
@@ -0,0 +1,70 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+
+define	M	256
+define	N	256
+
+task	lsq
+
+# This procedure fits a natural cubic spline to an array of n data points.
+# The system being solved is a tridiagonal matrix of n+2 rows.  The system
+# is solved by Lawsons and Hansons routine HFTI, which solves a general
+# m by n linear system of equations.  This is enormous overkill for this
+# problem (see "band.x"), but serves to give timing estimates for the code.
+
+
+procedure lsq()
+
+real	a[M,N], b[M], tau, rnorm, h[N], g[N]
+int	krank, ip[N]
+int	i, j, m, n, geti()
+real	marktime, cptime()
+
+begin
+	m = min (M, geti ("npts"))	# size of matrix
+	n = min (N, m)
+	tau = 1e-6
+
+	do j = 1, n			# set up b-spline matrix
+	    do i = 1, m
+		a[i,j] = 0.
+
+	a[1,1] = 6.			# first row
+	a[1,2] = -12.
+	a[1,3] = 6.
+
+	a[m,n] = 6.			# last row
+	a[m,n-1] = -12.
+	a[m,n-2] = 6.
+
+	do j = 2, m-1 {			# tridiagonal elements
+	    a[j,j-1] = 1.
+	    a[j,j] = 4.
+	    a[j,j+1] = 1.
+	}
+
+	b[1] = 0.			# natural spline bndry conditions
+	b[m] = 0.
+
+	do i = 2, m-1			# set up data vector
+	    b[i] = 100.
+
+	marktime = cptime()
+	call hfti (a,M,m,n, b,1,1, tau, krank,rnorm, h,g,ip)
+
+	call printf ("took %8.2f cpu seconds (krank=%d, rnorm=%g)\n")
+	    call pargr (cptime() - marktime)
+	    call pargi (krank)
+	    call pargr (rnorm)
+
+	call printf ("selected coefficients:\n")
+	for (i=1;  i <= m;) {		# print first, last 4 coeff
+	    call printf ("%8d%15.5f\n")
+		call pargi (i)
+		call pargr (b[i])
+	    if (i == 4)
+		i = max(i+1, m-3)
+	    else
+		i = i + 1
+	}
+end
diff --git a/math/llsq/progs/prog1.f b/math/llsq/progs/prog1.f
new file mode 100644
index 00000000..13f64929
--- /dev/null
+++ b/math/llsq/progs/prog1.f
@@ -0,0 +1,124 @@
+c     prog1
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	 demonstrate algorithms hft and hs1 for solving least squares
+c     problems and algorithm cov for omputing the associated covariance
+c     matrices.
+c
+      dimension a(8,8),h(8),b(8)
+      real gen,anoise
+      double precision sm
+      data mda/8/
+c
+	   do 180 noise=1,2
+	   anoise=0.
+	   if (noise.eq.2) anoise=1.e-4
+	   write (6,230)
+	   write (6,240) anoise
+c     initialize the data generation function
+c    ..
+	   dummy=gen(-1.)
+		do 180 mn1=1,6,5
+		mn2=mn1+2
+		     do 180 m=mn1,mn2
+		     do 180 n=mn1,mn2
+		     np1=n+1
+		     write (6,250) m,n
+c     generate data
+c    ..
+			  do 10 i=1,m
+			       do 10 j=1,n
+   10			       a(i,j)=gen(anoise)
+			  do 20 i=1,m
+   20			  b(i)=gen(anoise)
+		     if(m .lt. n) go to 180
+c
+c     ******  begin algorithm hft  ******
+c    ..
+			  do 30 j=1,n
+   30		   call h12 (1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j)
+c    ..
+c     the algorithm 'hft' is completed.
+c
+c     ******  begin algorithm hs1  ******
+c     apply the transformations  q(n)...q(1)=q to b
+c     replacing the previous contents of the array, b .
+c    ..
+			  do 40 j=1,n
+   40			      call h12 (2,j,j+1,m,a(1,j),1,h(j),b,1,1,1)
+c     solve the triangular system for the solution x.
+c     store x in the array, b .
+c    ..
+			  do 80 k=1,n
+			  i=np1-k
+			  sm=0.d0
+			  if (i.eq.n) go to 60
+			  ip1=i+1
+			       do 50 j=ip1,n
+   50			       sm=sm+a(i,j)*dble(b(j))
+   60			  if (a(i,i)) 80,70,80
+   70			  write (6,260)
+			  go to 180
+   80			  b(i)=(b(i)-sm)/a(i,i)
+c      compute length of residual vector.
+c     ..
+		     srsmsq=0.
+		     if (n.eq.m) go to 100
+		     mmn=m-n
+			  do 90 j=1,mmn
+			  npj=n+j
+   90			  srsmsq=srsmsq+b(npj)**2
+		     srsmsq=sqrt(srsmsq)
+c     ******  begin algorithm  cov  ******
+c     compute unscaled covariance matrix   ((a**t)*a)**(-1)
+c    ..
+  100			  do 110 j=1,n
+  110			  a(j,j)=1./a(j,j)
+		     if (n.eq.1) go to 140
+		     nm1=n-1
+			  do 130 i=1,nm1
+			  ip1=i+1
+			       do 130 j=ip1,n
+			       jm1=j-1
+			       sm=0.d0
+				    do 120 l=i,jm1
+  120				    sm=sm+a(i,l)*dble(a(l,j))
+  130			       a(i,j)=-sm*a(j,j)
+c    ..
+c     the upper triangle of a has been inverted
+c     upon itself.
+  140			  do 160 i=1,n
+			       do 160 j=i,n
+			       sm=0.d0
+				    do 150 l=j,n
+  150				    sm=sm+a(i,l)*dble(a(j,l))
+  160			       a(i,j)=sm
+c    ..
+c     the upper triangular part of the
+c     symmetric matrix (a**t*a)**(-1) has
+c     replaced the upper triangular part of
+c     the a array.
+		     write (6,200) (i,b(i),i=1,n)
+		     write (6,190) srsmsq
+		     write (6,210)
+			  do 170 i=1,n
+  170			  write (6,220) (i,j,a(i,j),j=i,n)
+  180		     continue
+      stop
+  190 format (1h0,8x,17hresidual length =,e12.4)
+  200 format (1h0,8x,34hestimated parameters,  x=a**(+)*b,,22h computed
+     1by 'hft,hs1'//(9x,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8))
+  210 format (1h0,8x,31hcovariance matrix (unscaled) of,22h estimated pa
+     1rameters.,19h computed by 'cov'./1x)
+  220 format (9x,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8)
+  230 format (52h1 prog1.    this program demonstrates the algorithms,19
+     1h hft, hs1, and cov.//40h caution.. 'prog1' does no checking for ,
+     252hnear rank deficient matrices.  results in such cases,20h may be
+     3 meaningless.,/34h such cases are handled by 'prog2',11h or 'pro
+     4')
+  240 format (1h0,54hthe relative noise level of the generated data will
+     1 be,e16.4)
+  250 format (1h0////9h0   m   n/1x,2i4)
+  260 format (1h0,8x,36h******  terminating this case due to,37h a divis
+     1or being exactly zero. ******)
+      end
diff --git a/math/llsq/progs/prog2.f b/math/llsq/progs/prog2.f
new file mode 100644
index 00000000..5fb601cc
--- /dev/null
+++ b/math/llsq/progs/prog2.f
@@ -0,0 +1,125 @@
+c     prog2
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	 demonstrate algorithm	hfti  for solving least squares problems
+c     and algorithm  cov  for computing the associated unscaled
+c     covariance matrix.
+c
+      dimension a(8,8),h(8),b(8),g(8)
+      real gen,anoise
+      integer ip(8)
+      double precision sm
+      data mda /8/
+c
+	  do 180 noise=1,2
+	  anorm=500.
+	  anoise=0.
+	  tau=.5
+	  if (noise.eq.1) go to 10
+	  anoise=1.e-4
+	  tau=anorm*anoise*10.
+   10	  continue
+c     initialize the data generation function
+c    ..
+	  dummy=gen(-1.)
+	  write (6,230)
+	  write (6,240) anoise,anorm,tau
+c
+	      do 180 mn1=1,6,5
+	      mn2=mn1+2
+		  do 180 m=mn1,mn2
+		      do 180 n=mn1,mn2
+		      write (6,250) m,n
+c     generate data
+c    ..
+			  do 20 i=1,m
+			      do 20 j=1,n
+   20			      a(i,j)=gen(anoise)
+			  do 30 i=1,m
+   30			  b(i)=gen(anoise)
+c
+c     ****** call hfti	 ******
+c
+		      call hfti(a,mda,m,n,b,1,1,tau,krank,srsmsq,h,g,ip)
+c
+c
+		      write (6,260) krank
+		      write (6,200) (i,b(i),i=1,n)
+		      write (6,190) srsmsq
+		      if (krank.lt.n) go to 180
+c     ******  algorithm cov bigins here  ******
+c    ..
+			  do 40 j=1,n
+   40			  a(j,j)=1./a(j,j)
+		      if (n.eq.1) go to 70
+		      nm1=n-1
+			  do 60 i=1,nm1
+			  ip1=i+1
+			      do 60 j=ip1,n
+			      jm1=j-1
+			      sm=0.d0
+				  do 50 l=i,jm1
+   50				  sm=sm+a(i,l)*dble(a(l,j))
+   60			      a(i,j)=-sm*a(j,j)
+c    ..
+c     the upper triangle of a has been inverted
+c     upon itself.
+   70			  do 90 i=1,n
+			      do 90 j=i,n
+			      sm=0.d0
+				  do 80 l=j,n
+   80				  sm=sm+a(i,l)*dble(a(j,l))
+   90			      a(i,j)=sm
+		      if (n.lt.2) go to 160
+			  do 150 ii=2,n
+			  i=n+1-ii
+			  if (ip(i).eq.i) go to 150
+			  k=ip(i)
+			  tmp=a(i,i)
+			  a(i,i)=a(k,k)
+			  a(k,k)=tmp
+			  if (i.eq.1) go to 110
+			      do 100 l=2,i
+			      tmp=a(l-1,i)
+			      a(l-1,i)=a(l-1,k)
+  100			      a(l-1,k)=tmp
+  110			  ip1=i+1
+			  km1=k-1
+			  if (ip1.gt.km1) go to 130
+			      do 120 l=ip1,km1
+			      tmp=a(i,l)
+			      a(i,l)=a(l,k)
+  120			      a(l,k)=tmp
+  130			  if (k.eq.n) go to 150
+			  kp1=k+1
+			      do 140 l=kp1,n
+			      tmp=a(i,l)
+			      a(i,l)=a(k,l)
+  140			      a(k,l)=tmp
+  150			  continue
+  160		      continue
+c    ..
+c     covariance has been computed and repermuted.
+c     the upper triangular part of the
+c     symmetric matrix (a**t*a)**(-1) has
+c     replaced the upper triangular part of
+c     the a array.
+		      write (6,210)
+			  do 170 i=1,n
+  170			  write (6,220) (i,j,a(i,j),j=i,n)
+  180		      continue
+      stop
+  190 format (1h0,8x,17hresidual length =,e12.4)
+  200 format (1h0,8x,34hestimated parameters,  x=a**(+)*b,,22h computed
+     1by 'hfti'   //(9x,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8))
+  210 format (1h0,8x,31hcovariance matrix (unscaled) of,22h estimated pa
+     1rameters.,19h computed by 'cov'./1x)
+  220 format (9x,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8)
+  230 format (52h1 prog2.     this program demonstates the algorithms,16
+     1h hfti  and  cov.)
+  240 format (1h0,54hthe relative noise level of the generated data will
+     1 be,e16.4/33h0the matrix norm is approximately,e12.4/43h0the absol
+     2ute pseudorank tolerance, tau, is,e12.4)
+  250 format (1h0////9h0   m   n/1x,2i4)
+  260 format (1h0,8x,12hpseudorank =,i4)
+      end
diff --git a/math/llsq/progs/prog3.f b/math/llsq/progs/prog3.f
new file mode 100644
index 00000000..290f8161
--- /dev/null
+++ b/math/llsq/progs/prog3.f
@@ -0,0 +1,138 @@
+c     prog3
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   demonstrate the use of subroutine   svdrs  to compute the
+c     singular value decomposition of a matrix, a, and solve a least
+c     squares problem,	a*x=b.
+c
+c     the s.v.d.  a= u*(s,0)**t*v**t is
+c     computed so that..
+c     (1)  u**t*b replaces the m+1 st. col. of	b.
+c
+c     (2)  u**t   replaces the m by m identity in
+c     the first m cols. of  b.
+c
+c     (3) v replaces the first n rows and cols.
+c     of a.
+c
+c     (4) the diagonal entries of the s matrix
+c     replace the first n entries of the array s.
+c
+c     the array s( ) must be dimensioned at least 3*n .
+c
+      dimension a(8,8),b(8,9),s(24),x(8),aa(8,8)
+      real gen,anoise
+      double precision sm
+      data mda,mdb/8,8/
+c
+      do 150 noise=1,2
+      anoise = 0.
+      rho = 1.e-3
+      if(noise .eq. 1) go to 5
+      anoise = 1.e-4
+      rho = 10. * anoise
+    5 continue
+      write(6,230)
+      write(6,240) anoise,rho
+c     initialize data generation function
+c    ..
+      dummy= gen(-1.)
+c
+	   do 150 mn1=1,6,5
+	   mn2=mn1+2
+		do 150 m=mn1,mn2
+		do 150 n=mn1,mn2
+		write (6,160) m,n
+		     do 20 i=1,m
+			  do 10 j=1,m
+   10			  b(i,j)=0.
+		     b(i,i)=1.
+			  do 20 j=1,n
+			  a(i,j)= gen(anoise)
+   20			  aa(i,j)=a(i,j)
+		     do 30 i=1,m
+   30		     b(i,m+1)= gen(anoise)
+c
+c     the arrays are now filled in..
+c     compute the s.v.d.
+c     ******************************************************
+		call svdrs (a,mda,m,n,b,mdb,m+1,s)
+c     ******************************************************
+c
+		write (6,170)
+		write (6,220) (i,s(i),i=1,n)
+		write (6,180)
+		write (6,220) (i,b(i,m+1),i=1,m)
+c
+c     test for disparity of ratio of singular values.
+c    ..
+		krank=n
+		tau=rho*s(1)
+		     do 40 i=1,n
+		     if (s(i).le.tau) go to 50
+   40		     continue
+		go to 55
+   50		krank=i-1
+   55 write(6,250) tau, krank
+c     compute solution vector assuming pseudorank is krank
+c     ..
+   60		     do 70 i=1,krank
+   70		     b(i,m+1)=b(i,m+1)/s(i)
+		     do 90 i=1,n
+		     sm=0.d0
+			  do 80 j=1,krank
+   80			  sm=sm+a(i,j)*dble(b(j,m+1))
+   90		     x(i)=sm
+c     compute predicted norm of residual vector.
+c     ..
+		srsmsq=0.
+		if (krank.eq.m) go to 110
+		kp1=krank+1
+		     do 100 i=kp1,m
+  100		     srsmsq=srsmsq+b(i,m+1)**2
+		srsmsq=sqrt(srsmsq)
+c
+  110		continue
+		write (6,190)
+		write (6,220) (i,x(i),i=1,n)
+		write (6,200) srsmsq
+c     compute the frobenius norm of  a**t- v*(s,0)*u**t.
+c
+c     compute  v*s  first.
+c
+		minmn=min0(m,n)
+		     do 120 j=1,minmn
+			  do 120 i=1,n
+  120			  a(i,j)=a(i,j)*s(j)
+		dn=0.
+		     do 140 j=1,m
+			  do 140 i=1,n
+			  sm=0.d0
+			       do 130 l=1,minmn
+  130			       sm=sm+a(i,l)*dble(b(l,j))
+c     computed difference of (i,j) th entry
+c     of  a**t-v*(s,0)*u**t.
+c     ..
+			  t=aa(j,i)-sm
+  140			  dn=dn+t**2
+		dn=sqrt(dn)/(sqrt(float(n))*s(1))
+		write (6,210) dn
+  150		continue
+      stop
+  160 format (1h0////9h0   m   n/1x,2i4)
+  170 format (1h0,8x,25hsingular values of matrix)
+  180 format (1h0,8x,30htransformed right side, u**t*b)
+  190 format (1h0,8x,33hestimated parameters,  x=a**(+)*b,21h computed b
+     1y 'svdrs' )
+  200 format (1h0,8x,24hresidual vector length =,e12.4)
+  210 format (1h0,8x,43hfrobenius norm (a-u*(s,0)**t*v**t)/(sqrt(n),
+     *22h*spectral norm of a) =,e12.4)
+  220 format (9x,i5,e16.8,i5,e16.8,i5,e16.8,i5,e16.8,i5,e16.8)
+  230 format(51h1prog3.     this program demonstrates the algorithm,
+     * 9h, svdrs. )
+  240 format(55h0the relative noise level of the generated data will be,
+     *e16.4/44h0the relative tolerance, rho, for pseudorank,
+     *17h determination is,e16.4)
+  250 format(1h0,8x,36habsolute pseudorank tolerance, tau =,
+     *e12.4,10x,12hpseudorank =,i5)
+      end
diff --git a/math/llsq/progs/prog4.f b/math/llsq/progs/prog4.f
new file mode 100644
index 00000000..090b1865
--- /dev/null
+++ b/math/llsq/progs/prog4.f
@@ -0,0 +1,22 @@
+c     prog4
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   demonstrate singular value analysis.
+c
+      dimension a(15,5),b(15),sing(15)
+      data names/1h /
+c
+      read (5,10) ((a(i,j),j=1,5),b(i),i=1,15)
+      write (6,20)
+      write (6,30) ((a(i,j),j=1,5),b(i),i=1,15)
+      write (6,40)
+c
+      call sva (a,15,15,5,15,b,sing,names,1,d)
+c
+      stop
+   10 format (6f12.0)
+   20 format (46h1prog4.    demonstrate singular value analysis/53h list
+     1ing of input matrix, a, and vector, b, follows..)
+   30 format (1h /(5f12.8,f20.4))
+   40 format (1h1)
+      end
diff --git a/math/llsq/progs/prog5.f b/math/llsq/progs/prog5.f
new file mode 100644
index 00000000..5484b846
--- /dev/null
+++ b/math/llsq/progs/prog5.f
@@ -0,0 +1,146 @@
+c     prog5
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   example of the use of subroutines bndacc and bndsol to solve
+c     sequentially the banded least squares problem which arises in
+c     spline curve fitting.
+c
+      dimension x(12),y(12),b(10),g(12,5),c(12),q(4),cov(12)
+c
+c	 define functions p1 and p2 to be used in constructing
+c	 cubic splines over uniformly spaced breakpoints.
+c
+      p1(t)=.25*t**2*t
+      p2(t)=-(1.-t)**2*(1.+t)*.75+1.
+      zero=0.
+c
+      write (6,230)
+      mdg=12
+      nband=4
+      m=12
+c				   set ordinates of data
+      y(1)=2.2
+      y(2)=4.0
+      y(3)=5.0
+      y(4)=4.6
+      y(5)=2.8
+      y(6)=2.7
+      y(7)=3.8
+      y(8)=5.1
+      y(9)=6.1
+      y(10)=6.3
+      y(11)=5.0
+      y(12)=2.0
+c				   set abcissas of data
+	   do 10 i=1,m
+   10	   x(i)=2*i
+c
+c     begin loop thru cases using increasing nos of breakpoints.
+c
+	   do 150 nbp=5,10
+	   nc=nbp+2
+c				   set breakpoints
+	   b(1)=x(1)
+	   b(nbp)=x(m)
+	   h=(b(nbp)-b(1))/float(nbp-1)
+	   if (nbp.le.2) go to 30
+		do 20 i=3,nbp
+   20		b(i-1)=b(i-2)+h
+   30	   continue
+	   write (6,240) nbp,(b(i),i=1,nbp)
+c
+c     initialize ir and ip before first call to bndacc.
+c
+	   ir=1
+	   ip=1
+	   i=1
+	   jt=1
+   40	   mt=0
+   50	   continue
+	   if (x(i).gt.b(jt+1)) go to 60
+c
+c			 set row  for ith data point
+c
+	   u=(x(i)-b(jt))/h
+	   ig=ir+mt
+	   g(ig,1)=p1(1.-u)
+	   g(ig,2)=p2(1.-u)
+	   g(ig,3)=p2(u)
+	   g(ig,4)=p1(u)
+	   g(ig,5)=y(i)
+	   mt=mt+1
+	   if (i.eq.m) go to 60
+	   i=i+1
+	   go to 50
+c
+c		    send block of data to processor
+c
+   60	   continue
+	   call bndacc (g,mdg,nband,ip,ir,mt,jt)
+	   if (i.eq.m) go to 70
+	   jt=jt+1
+	   go to 40
+c		    compute solution c()
+   70	   continue
+	   call bndsol (1,g,mdg,nband,ip,ir,c,nc,rnorm)
+c		   write solution coefficients
+	   write (6,180) (c(l),l=1,nc)
+	   write (6,210) rnorm
+c
+c	       compute and print x,y,yfit,r=y-yfit
+c
+	   write (6,160)
+	   jt=1
+		do 110 i=1,m
+   80		if (x(i).le.b(jt+1)) go to 90
+		jt=jt+1
+		go to 80
+c
+   90		u=(x(i)-b(jt))/h
+		q(1)=p1(1.-u)
+		q(2)=p2(1.-u)
+		q(3)=p2(u)
+		q(4)=p1(u)
+		yfit=zero
+		     do 100 l=1,4
+		     ic=jt-1+l
+  100		     yfit=yfit+c(ic)*q(l)
+		r=y(i)-yfit
+		write (6,170) i,x(i),y(i),yfit,r
+  110		continue
+c
+c     compute residual vector norm.
+c
+	   if (m.le.nc) go to 150
+	   sigsq=(rnorm**2)/float(m-nc)
+	   sigfac=sqrt(sigsq)
+	   write (6,220) sigfac
+	   write (6,200)
+c
+c     compute and print cols. of covariance.
+c
+		do 140 j=1,nc
+		     do 120 i=1,nc
+  120		     cov(i)=zero
+		cov(j)=1.
+		call bndsol (2,g,mdg,nband,ip,ir,cov,nc,rdummy)
+		call bndsol (3,g,mdg,nband,ip,ir,cov,nc,rdummy)
+c
+c    compute the jth col. of the covariance matrix.
+		     do 130 i=1,nc
+  130		     cov(i)=cov(i)*sigsq
+  140		write (6,190) (l,j,cov(l),l=1,nc)
+  150	   continue
+      stop
+  160 format (4h0  i,8x,1hx,10x,1hy,6x,4hyfit,4x,8hr=y-yfit/1x)
+  170 format (1x,i3,4x,f6.0,4x,f6.2,4x,f6.2,4x,f8.4)
+  180 format (4h0c =,10f10.5/(4x,10f10.5))
+  190 format (4(02x,2i5,e15.7))
+  200 format (46h0covariance matrix of the spline coefficients.)
+  210 format (9h0rnorm	=,e15.8)
+  220 format (9h0sigfac =,e15.8)
+  230 format (50h1prog5.    execute a sequence of cubic spline fits,27h
+     1to a discrete set of data.)
+  240 format (1x////,11h0new case..,/47h0number of breakpoints, includin
+     1g endpoints, is,i5/14h0breakpoints =,10f10.5,/(14x,10f10.5))
+      end
diff --git a/math/llsq/progs/prog6.f b/math/llsq/progs/prog6.f
new file mode 100644
index 00000000..71364ec4
--- /dev/null
+++ b/math/llsq/progs/prog6.f
@@ -0,0 +1,116 @@
+c     prog6
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 15
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   demonstrate the use of the subroutine   ldp	for least
+c     distance programming by solving the constrained line data fitting
+c     problem of chapter 23.
+c
+      dimension e(4,2), f(4), g(3,2), h(3), g2(3,2), h2(3), x(2), z(2),
+     1w(4)
+      dimension wldp(21), s(6), t(4)
+      integer index(3)
+c
+      write (6,110)
+      mde=4
+      mdgh=3
+c
+      me=4
+      mg=3
+      n=2
+c		 define the least squares and constraint matrices.
+      t(1)=0.25
+      t(2)=0.5
+      t(3)=0.5
+      t(4)=0.8
+c
+      w(1)=0.5
+      w(2)=0.6
+      w(3)=0.7
+      w(4)=1.2
+c
+	  do 10 i=1,me
+	  e(i,1)=t(i)
+	  e(i,2)=1.
+   10	  f(i)=w(i)
+c
+      g(1,1)=1.
+      g(1,2)=0.
+      g(2,1)=0.
+      g(2,2)=1.
+      g(3,1)=-1.
+      g(3,2)=-1.
+c
+      h(1)=0.
+      h(2)=0.
+      h(3)=-1.
+c
+c     compute the singular value decomposition of the matrix, e.
+c
+      call svdrs (e,mde,me,n,f,1,1,s)
+c
+      write (6,120) ((e(i,j),j=1,n),i=1,n)
+      write (6,130) f,(s(j),j=1,n)
+c
+c     generally rank determination and levenberg-marquardt
+c     stabilization could be inserted here.
+c
+c	 define the constraint matrix for the z coordinate system.
+	  do 30 i=1,mg
+	      do 30 j=1,n
+	      sm=0.
+		  do 20 l=1,n
+   20		  sm=sm+g(i,l)*e(l,j)
+   30	      g2(i,j)=sm/s(j)
+c	  define constraint rt side for the z coordinate system.
+	  do 50 i=1,mg
+	  sm=0.
+	      do 40 j=1,n
+   40	      sm=sm+g2(i,j)*f(j)
+   50	  h2(i)=h(i)-sm
+c
+      write (6,140) ((g2(i,j),j=1,n),i=1,mg)
+      write (6,150) h2
+c
+c			 solve the constrained problem in z-coordinates.
+c
+      call ldp (g2,mdgh,mg,n,h2,z,znorm,wldp,index,mode)
+c
+      write (6,200) mode,znorm
+      write (6,160) z
+c
+c		     transform back from z-coordinates to x-coordinates.
+	  do 60 j=1,n
+   60	  z(j)=(z(j)+f(j))/s(j)
+	  do 80 i=1,n
+	  sm=0.
+	      do 70 j=1,n
+   70	      sm=sm+e(i,j)*z(j)
+   80	  x(i)=sm
+      res=znorm**2
+      np1=n+1
+	  do 90 i=np1,me
+   90	  res=res+f(i)**2
+      res=sqrt(res)
+c				     compute the residuals.
+	  do 100 i=1,me
+  100	  f(i)=w(i)-x(1)*t(i)-x(2)
+      write (6,170) (x(j),j=1,n)
+      write (6,180) (i,f(i),i=1,me)
+      write (6,190) res
+      stop
+c
+  110 format (46h0prog6..  example of constrained curve fitting,26h usin
+     1g the subroutine ldp.,/43h0related intermediate quantities are giv
+     2en.)
+  120 format (10h0v =	   ,2f10.5/(10x,2f10.5))
+  130 format (10h0f tilda =,4f10.5/10h0s =	,2f10.5)
+  140 format (10h0g tilda =,2f10.5/(10x,2f10.5))
+  150 format (10h0h tilda =,3f10.5)
+  160 format (10h0z =	   ,2f10.5)
+  170 format (52h0the coeficients of the fitted line f(t)=x(1)*t+x(2),12
+     1h are x(1) = ,f10.5,14h	and x(2) = ,f10.5)
+  180 format (30h0the consecutive residuals are/1x,4(i10,f10.5))
+  190 format (23h0the residuals norm is ,f10.5)
+  200 format (18h0mode (from ldp) =,i3,2x,7hznorm =,f10.5)
+c
+      end
diff --git a/math/llsq/qrbd.f b/math/llsq/qrbd.f
new file mode 100644
index 00000000..d9a12a24
--- /dev/null
+++ b/math/llsq/qrbd.f
@@ -0,0 +1,208 @@
+c     subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   qr algorithm for singular values of a bidiagonal matrix.
+c
+c     the bidiagonal matrix
+c
+c			(q1,e2,0...    )
+c			(   q2,e3,0... )
+c		 d=	(	.      )
+c			(	  .   0)
+c			(	    .en)
+c			(	   0,qn)
+c
+c		  is pre and post multiplied by
+c		  elementary rotation matrices
+c		  ri and pi so that
+c
+c		  rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sn)
+c
+c		  to within working accuracy.
+c
+c  1. ei and qi occupy e(i) and q(i) as input.
+c
+c  2. rm...r1*c replaces 'c' in storage as output.
+c
+c  3. v*p1**(t)...pm**(t) replaces 'v' in storage as output.
+c
+c  4. si occupies q(i) as output.
+c
+c  5. the si's are nonincreasing and nonnegative.
+c
+c     this code is based on the paper and 'algol' code..
+c ref..
+c  1. reinsch,c.h. and golub,g.h. 'singular value decomposition
+c     and least squares solutions' (numer. math.), vol. 14,(1970).
+c
+      subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc)
+      logical wntv ,havers,fail
+      dimension q(nn),e(nn),v(mdv,nn),c(mdc,ncc)
+      zero=0.
+      one=1.
+      two=2.
+c
+      n=nn
+      ipass=1
+      if (n.le.0) return
+      n10=10*n
+      wntv=nrv.gt.0
+      havers=ncc.gt.0
+      fail=.false.
+      nqrs=0
+      e(1)=zero
+      dnorm=zero
+	   do 10 j=1,n
+   10	   dnorm=amax1(abs(q(j))+abs(e(j)),dnorm)
+	   do 200 kk=1,n
+	   k=n+1-kk
+c
+c     test for splitting or rank deficiencies..
+c	  first make test for last diagonal term, q(k), being small.
+   20	    if(k.eq.1) go to 50
+	    if(diff(dnorm+q(k),dnorm)) 50,25,50
+c
+c     since q(k) is small we will make a special pass to
+c     transform e(k) to zero.
+c
+   25	   cs=zero
+	   sn=-one
+		do 40 ii=2,k
+		i=k+1-ii
+		f=-sn*e(i+1)
+		e(i+1)=cs*e(i+1)
+		call g1 (q(i),f,cs,sn,q(i))
+c	  transformation constructed to zero position (i,k).
+c
+		if (.not.wntv) go to 40
+		     do 30 j=1,nrv
+   30		     call g2 (cs,sn,v(j,i),v(j,k))
+c	       accumulate rt. transformations in v.
+c
+   40		continue
+c
+c	  the matrix is now bidiagonal, and of lower order
+c	  since e(k) .eq. zero..
+c
+   50		do 60 ll=1,k
+		l=k+1-ll
+		if(diff(dnorm+e(l),dnorm)) 55,100,55
+   55		if(diff(dnorm+q(l-1),dnorm)) 60,70,60
+   60		continue
+c     this loop can't complete since e(1) = zero.
+c
+	   go to 100
+c
+c	  cancellation of e(l), l.gt.1.
+   70	   cs=zero
+	   sn=-one
+		do 90 i=l,k
+		f=-sn*e(i)
+		e(i)=cs*e(i)
+		if(diff(dnorm+f,dnorm)) 75,100,75
+   75		call g1 (q(i),f,cs,sn,q(i))
+		if (.not.havers) go to 90
+		     do 80 j=1,ncc
+   80		     call g2 (cs,sn,c(i,j),c(l-1,j))
+   90		continue
+c
+c	  test for convergence..
+  100	   z=q(k)
+	   if (l.eq.k) go to 170
+c
+c	  shift from bottom 2 by 2 minor of b**(t)*b.
+	   x=q(l)
+	   y=q(k-1)
+	   g=e(k-1)
+	   h=e(k)
+	   f=((y-z)*(y+z)+(g-h)*(g+h))/(two*h*y)
+	   g=sqrt(one+f**2)
+	   if (f.lt.zero) go to 110
+	   t=f+g
+	   go to 120
+  110	   t=f-g
+  120	   f=((x-z)*(x+z)+h*(y/t-h))/x
+c
+c	  next qr sweep..
+	   cs=one
+	   sn=one
+	   lp1=l+1
+		do 160 i=lp1,k
+		g=e(i)
+		y=q(i)
+		h=sn*g
+		g=cs*g
+		call g1 (f,h,cs,sn,e(i-1))
+		f=x*cs+g*sn
+		g=-x*sn+g*cs
+		h=y*sn
+		y=y*cs
+		if (.not.wntv) go to 140
+c
+c	       accumulate rotations (from the right) in 'v'
+		     do 130 j=1,nrv
+  130		     call g2 (cs,sn,v(j,i-1),v(j,i))
+  140		call g1 (f,h,cs,sn,q(i-1))
+		f=cs*g+sn*y
+		x=-sn*g+cs*y
+		if (.not.havers) go to 160
+		     do 150 j=1,ncc
+  150		     call g2 (cs,sn,c(i-1,j),c(i,j))
+c	       apply rotations from the left to
+c	       right hand sides in 'c'..
+c
+  160		continue
+	   e(l)=zero
+	   e(k)=f
+	   q(k)=x
+	   nqrs=nqrs+1
+	   if (nqrs.le.n10) go to 20
+c	   return to 'test for splitting'.
+c
+	   fail=.true.
+c     ..
+c     cutoff for convergence failure. 'nqrs' will be 2*n usually.
+  170	   if (z.ge.zero) go to 190
+	   q(k)=-z
+	   if (.not.wntv) go to 190
+		do 180 j=1,nrv
+  180		v(j,k)=-v(j,k)
+  190	   continue
+c	  convergence. q(k) is made nonnegative..
+c
+  200	   continue
+      if (n.eq.1) return
+	   do 210 i=2,n
+	   if (q(i).gt.q(i-1)) go to 220
+  210	   continue
+      if (fail) ipass=2
+      return
+c     ..
+c     every singular value is in order..
+  220	   do 270 i=2,n
+	   t=q(i-1)
+	   k=i-1
+		do 230 j=i,n
+		if (t.ge.q(j)) go to 230
+		t=q(j)
+		k=j
+  230		continue
+	   if (k.eq.i-1) go to 270
+	   q(k)=q(i-1)
+	   q(i-1)=t
+	   if (.not.havers) go to 250
+		do 240 j=1,ncc
+		t=c(i-1,j)
+		c(i-1,j)=c(k,j)
+  240		c(k,j)=t
+  250	   if (.not.wntv) go to 270
+		do 260 j=1,nrv
+		t=v(j,i-1)
+		v(j,i-1)=v(j,k)
+  260		v(j,k)=t
+  270	   continue
+c	  end of ordering algorithm.
+c
+      if (fail) ipass=2
+      return
+      end
diff --git a/math/llsq/sva.f b/math/llsq/sva.f
new file mode 100644
index 00000000..ee7f893d
--- /dev/null
+++ b/math/llsq/sva.f
@@ -0,0 +1,198 @@
+      subroutine sva (a,mda,m,n,mdata,b,sing,names,iscale,d)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c               singular value analysis printout
+c
+c     iscale       set by user to 1, 2, or 3 to select column scaling
+c                  option.
+c                  1   subr will use identity scaling and ignore the d()
+c                      array.
+c                  2   subr will scale nonzero cols to have unit euclide
+c                      length and will store reciprocal lengths of
+c                      original nonzero cols in d().
+c                  3   user supplies col scale factors in d(). subr
+c                      will mult col j by d(j) and remove the scaling
+c                      from the soln at the end.
+c
+      dimension  a(mda,n), b(m), sing(01),d(n)
+c     sing(3*n)
+      dimension names(n)
+      logical test
+      double precision   sb, dzero
+      dzero=0.d0
+      one=1.
+      zero=0.
+      if (m.le.0 .or. n.le.0) return
+      np1=n+1
+      write (6,270)
+      write (6,260) m,n,mdata
+      go to (60,10,10), iscale
+c
+c     apply column scaling to a
+c
+   10      do 50 j=1,n
+           a1=d(j)
+           go to (20,20,40), iscale
+   20      sb=dzero
+                do 30 i=1,m
+   30           sb=sb+dble(a(i,j))**2
+           a1=dsqrt(sb)
+           if (a1.eq.zero) a1=one
+           a1=one/a1
+           d(j)=a1
+   40           do 50 i=1,m
+   50           a(i,j)=a(i,j)*a1
+      write (6,280) iscale,(d(j),j=1,n)
+      go to 70
+   60 continue
+      write (6,290)
+   70 continue
+c
+c     obtain  sing. value decomp. of scaled matrix.
+c
+c**********************************************************
+      call svdrs (a,mda,m,n,b,1,1,sing,ier)
+      if (ier.eq.0) goto 75
+         write (6,295)
+         return
+   75 continue
+c**********************************************************
+c
+c     print the v matrix.
+c
+      call mfeout (a,mda,n,n,names,1)
+      if (iscale.eq.1) go to 90
+c
+c     replace v by d*v in the array a()
+c
+           do 80 i=1,n
+                do 80 j=1,n
+   80           a(i,j)=d(i)*a(i,j)
+c
+c     g  now in  b array.  v now in a array.
+c
+c
+c     obtain summary output.
+c
+   90 continue
+      write (6,220)
+c
+c     compute cumulative sums of squares of components of
+c     g and store them in sing(i), i=minmn+1,...,2*minmn+1
+c
+      sb=dzero
+      minmn=min0(m,n)
+      minmn1=minmn+1
+      if (m.eq.minmn) go to 110
+           do 100 i=minmn1,m
+  100      sb=sb+dble(b(i))**2
+  110 sing(2*minmn+1)=sb
+           do 120 jj=1,minmn
+           j=minmn+1-jj
+           sb=sb+dble(b(j))**2
+           js=minmn+j
+  120      sing(js)=sb
+      a3=sing(minmn+1)
+      a4=sqrt(a3/float(max0(1,mdata)))
+      write (6,230) a3,a4
+c
+      nsol=0
+c
+c
+c
+           do 160 k=1,minmn
+           if (sing(k).eq.zero) go to 130
+           nsol=k
+           pi=b(k)/sing(k)
+           a1=one/sing(k)
+           a2=b(k)**2
+           a3=sing(minmn1+k)
+           a4=sqrt(a3/float(max0(1,mdata-k)))
+           test=sing(k).ge.100..or.sing(k).lt..001
+           if (test) write (6,240) k,sing(k),pi,a1,b(k),a2,a3,a4
+           if (.not.test) write (6,250) k,sing(k),pi,a1,b(k),a2,a3,a4
+           go to 140
+  130      write (6,240) k,sing(k)
+           pi=zero
+  140           do 150 i=1,n
+                a(i,k)=a(i,k)*pi
+  150           if (k.gt.1) a(i,k)=a(i,k)+a(i,k-1)
+  160      continue
+c
+c     compute and print values of ynorm and rnorm.
+c
+      write (6,300)
+      j=0
+      ysq=zero
+      go to 180
+  170 j=j+1
+      ysq=ysq+(b(j)/sing(j))**2
+  180 ynorm=sqrt(ysq)
+      js=minmn1+j
+      rnorm=sqrt(sing(js))
+      yl=-1000.
+      if (ynorm.gt.0.) yl=alog10(ynorm)
+      rl=-1000.
+      if (rnorm.gt.0.) rl=alog10(rnorm)
+      write (6,310) j,ynorm,rnorm,yl,rl
+      if (j.lt.nsol) go to 170
+c
+c     compute values of xnorm and rnorm for a sequence of values of
+c     the levenberg-marquardt parameter.
+c
+      if (sing(1).eq.zero) go to 210
+      el=alog10(sing(1))+one
+      el2=alog10(sing(nsol))-one
+      del=(el2-el)/20.
+      ten=10.
+      aln10=alog(ten)
+      write (6,320)
+           do 200 ie=1,21
+c     compute        alamb=10.**el
+           alamb=exp(aln10*el)
+           ys=0.
+           js=minmn1+nsol
+           rs=sing(js)
+                do 190 i=1,minmn
+                sl=sing(i)**2+alamb**2
+                ys=ys+(b(i)*sing(i)/sl)**2
+                rs=rs+(b(i)*(alamb**2)/sl)**2
+  190           continue
+           ynorm=sqrt(ys)
+           rnorm=sqrt(rs)
+           rl=-1000.
+           if (rnorm.gt.zero) rl=alog10(rnorm)
+           yl=-1000.
+           if (ynorm.gt.zero) yl=alog10(ynorm)
+           write (6,330) alamb,ynorm,rnorm,el,yl,rl
+           el=el+del
+  200      continue
+c
+c     print candidate solutions.
+c
+  210 if (nsol.ge.1) call mfeout (a,mda,n,nsol,names,2)
+      return
+  220 format (42h0 index  sing. value       p coef         ,48hrecip. s.
+     1v.             g coef            g**2  ,
+     234h            c.s.s.    n.s.r.c.s.s.)
+  230 format (1h ,5x,1h0,88x,1pe15.4,1pe17.4)
+  240 format (1h ,i6,e12.4,1p,(e15.4,4x,e15.4,4x,e15.4,4x,e15.4,4x,e15.4
+     1,2x,e15.4))
+  250 format (1h ,i6,f12.4,1p,(e15.4,4x,e15.4,4x,e15.4,4x,e15.4,4x,e15.4
+     1,2x,e15.4))
+  260 format (5h0m = ,i6,8h,   n = ,i4,12h,   mdata = ,i8)
+  270 format (45h0singular value analysis of the least squares,42h probl
+     1em,  a*x=b,  scaled as  (a*d)*y=b  .)
+  280 format (19h0scaling option no.,i2,18h.  d is a diagonal,46h matrix
+     1 with the following diagonal elements../(5x,10e12.4))
+  290 format (50h0scaling option no. 1.   d is the identity matrix./1x)
+  295 format (49h convergence failure in qr bidiagonal svd routine)
+  300 format (6h0index,13x,28h         ynorm         rnorm,14x,28h  log1
+     10(ynorm)  log10(rnorm)/1x)
+  310 format (1h ,i4,14x,2e14.5,14x,2f14.5)
+  320 format (54h0norms of solution and residual vectors for a range of,
+     144h values of the levenberg-marquardt parameter,9h, lambda./1h0,4x
+     2,42h        lambda         ynorm         rnorm,
+     345h log10(lambda)     log10(ynorm)  log10(rnorm))
+  330 format (5x,3e14.5,3f14.5)
+      end
diff --git a/math/llsq/svdrs.f b/math/llsq/svdrs.f
new file mode 100644
index 00000000..46dba02d
--- /dev/null
+++ b/math/llsq/svdrs.f
@@ -0,0 +1,205 @@
+c     subroutine svdrs (a,mda,mm,nn,b,mdb,nb,s,ier)
+c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1974 mar 1
+c     to appear in 'solving least squares problems', prentice-hall, 1974
+c	   singular value decomposition also treating right side vector.
+c
+c     the array s occupies 3*n cells.
+c     a occupies m*n cells
+c     b occupies m*nb cells.
+c
+c     special singular value decomposition subroutine.
+c     we have the m x n matrix a and the system a*x=b to solve.
+c     either m .ge. n  or  m .lt. n is permitted.
+c			the singular value decomposition
+c     a = u*s*v**(t) is made in such a way that one gets
+c	 (1) the matrix v in the first n rows and columns of a.
+c	 (2) the diagonal matrix of ordered singular values in
+c	     the first n cells of the array s(ip), ip .ge. 3*n.
+c	 (3) the matrix product u**(t)*b=g gets placed back in b.
+c	 (4) the user must complete the solution and do his own
+c	     singular value analysis.
+c     *******
+c     give special
+c     treatment to rows and columns which are entirely zero.  this
+c     causes certain zero sing. vals. to appear as exact zeros rather
+c     than as about eta times the largest sing. val.   it similarly
+c     cleans up the associated columns of u and v.
+c     method..
+c      1. exchange cols of a to pack nonzero cols to the left.
+c	  set n = no. of nonzero cols.
+c	  use locations a(1,nn),a(1,nn-1),...,a(1,n+1) to record the
+c	  col permutations.
+c      2. exchange rows of a to pack nonzero rows to the top.
+c	  quit packing if find n nonzero rows.	make same row exchanges
+c	  in b.  set m so that all nonzero rows of the permuted a
+c	  are in first m rows.	if m .le. n then all m rows are
+c	  nonzero.  if m .gt. n then	  the first n rows are known
+c	  to be nonzero,and rows n+1 thru m may be zero or nonzero.
+c      3. apply original algorithm to the m by n problem.
+c      4. move permutation record from a(,) to s(i),i=n+1,...,nn.
+c      5. build v up from  n by n  to  nn by nn  by placing ones on
+c	  the diagonal and zeros elsewhere.  this is only partly done
+c	  explicitly.  it is completed during step 6.
+c      6. exchange rows of v to compensate for col exchanges of step 2.
+c      7. place zeros in  s(i),i=n+1,nn  to represent zero sing vals.
+c
+      subroutine svdrs (a,mda,mm,nn,b,mdb,nb,s,ier)
+      dimension a(mda,nn),b(mdb,nb),s(nn,3)
+      zero=0.
+      ier = 0
+      one=1.
+c
+c			      begin.. special for zero rows and cols.
+c
+c			      pack the nonzero cols to the left
+c
+      n=nn
+      if (n.le.0.or.mm.le.0) return
+      j=n
+   10 continue
+	  do 20 i=1,mm
+	  if (a(i,j)) 50,20,50
+   20	  continue
+c
+c	  col j  is zero. exchange it with col n.
+c
+      if (j.eq.n) go to 40
+	  do 30 i=1,mm
+   30	  a(i,j)=a(i,n)
+   40 continue
+      a(1,n)=j
+      n=n-1
+   50 continue
+      j=j-1
+      if (j.ge.1) go to 10
+c			      if n=0 then a is entirely zero and svd
+c			      computation can be skipped
+      ns=0
+      if (n.eq.0) go to 240
+c			      pack nonzero rows to the top
+c			      quit packing if find n nonzero rows
+      i=1
+      m=mm
+   60 if (i.gt.n.or.i.ge.m) go to 150
+      if (a(i,i)) 90,70,90
+   70	  do 80 j=1,n
+	  if (a(i,j)) 90,80,90
+   80	  continue
+      go to 100
+   90 i=i+1
+      go to 60
+c			      row i is zero
+c			      exchange rows i and m
+  100 if(nb.le.0) go to 115
+	  do 110 j=1,nb
+	  t=b(i,j)
+	  b(i,j)=b(m,j)
+  110	  b(m,j)=t
+  115	  do 120 j=1,n
+  120	  a(i,j)=a(m,j)
+      if (m.gt.n) go to 140
+	  do 130 j=1,n
+  130	  a(m,j)=zero
+  140 continue
+c			      exchange is finished
+      m=m-1
+      go to 60
+c
+  150 continue
+c			      end.. special for zero rows and columns
+c			      begin.. svd algorithm
+c     method..
+c     (1)     reduce the matrix to upper bidiagonal form with
+c     householder transformations.
+c	   h(n)...h(1)aq(1)...q(n-2) = (d**t,0)**t
+c     where d is upper bidiagonal.
+c
+c     (2)     apply h(n)...h(1) to b.  here h(n)...h(1)*b replaces b
+c     in storage.
+c
+c     (3)     the matrix product w= q(1)...q(n-2) overwrites the first
+c     n rows of a in storage.
+c
+c     (4)     an svd for d is computed.  here k rotations ri and pi are
+c     computed so that
+c	   rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sm)
+c     to working accuracy.  the si are nonnegative and nonincreasing.
+c     here rk...r1*b overwrites b in storage while
+c     a*p1**(t)...pk**(t)  overwrites a in storage.
+c
+c     (5)     it follows that,with the proper definitions,
+c     u**(t)*b overwrites b, while v overwrites the first n row and
+c     columns of a.
+c
+      l=min0(m,n)
+c	      the following loop reduces a to upper bidiagonal and
+c	      also applies the premultiplying transformations to b.
+c
+	  do 170 j=1,l
+	  if (j.ge.m) go to 160
+	  call h12 (1,j,j+1,m,a(1,j),1,t,a(1,j+1),1,mda,n-j)
+	  call h12 (2,j,j+1,m,a(1,j),1,t,b,1,mdb,nb)
+  160	  if (j.ge.n-1) go to 170
+	  call h12 (1,j+1,j+2,n,a(j,1),mda,s(j,3),a(j+1,1),mda,1,m-j)
+  170	  continue
+c
+c     copy the bidiagonal matrix into the array s() for qrbd.
+c
+      if (n.eq.1) go to 190
+	  do 180 j=2,n
+	  s(j,1)=a(j,j)
+  180	  s(j,2)=a(j-1,j)
+  190 s(1,1)=a(1,1)
+c
+      ns=n
+      if (m.ge.n) go to 200
+      ns=m+1
+      s(ns,1)=zero
+      s(ns,2)=a(m,m+1)
+  200 continue
+c
+c     construct the explicit n by n product matrix, w=q1*q2*...*ql*i
+c     in the array a().
+c
+	  do 230 k=1,n
+	  i=n+1-k
+	  if(i.gt.min0(m,n-2)) go to 210
+	  call h12 (2,i+1,i+2,n,a(i,1),mda,s(i,3),a(1,i+1),1,mda,n-i)
+  210	      do 220 j=1,n
+  220	      a(i,j)=zero
+  230	  a(i,i)=one
+c
+c	   compute the svd of the bidiagonal matrix
+c
+      call qrbd (ipass,s(1,1),s(1,2),ns,a,mda,n,b,mdb,nb)
+c
+      go to (240,310), ipass
+  240 continue
+      if (ns.ge.n) go to 260
+      nsp1=ns+1
+	  do 250 j=nsp1,n
+  250	  s(j,1)=zero
+  260 continue
+      if (n.eq.nn) return
+      np1=n+1
+c				   move record of permutations
+c				   and store zeros
+	  do 280 j=np1,nn
+	  s(j,1)=a(1,j)
+	      do 270 i=1,n
+  270	      a(i,j)=zero
+  280	  continue
+c			      permute rows and set zero singular values.
+	  do 300 k=np1,nn
+	  i=s(k,1)
+	  s(k,1)=zero
+	      do 290 j=1,nn
+	      a(k,j)=a(i,j)
+  290	      a(i,j)=zero
+	  a(i,k)=one
+  300	  continue
+c			      end.. special for zero rows and columns
+      return
+c     convergence failure in qr bidiagonal svd routine
+  310 ier = 1
+      end