(defun qrbd (ipass q e nn v mdv nrv c mdc ncc) (declare (type fixnum ipass))
 (declare (type (simple-array double-float (*)) q))
 (declare (type (simple-array double-float (*)) e)) (declare (type fixnum nn))
 (declare (type (simple-array double-float (* *)) v))
 (declare (type fixnum mdv)) (declare (type fixnum nrv))
 (declare (type (simple-array double-float (* *)) c))
 (declare (type fixnum mdc)) (declare (type fixnum ncc))
 (prog
  ((wntv nil) (havers nil) (fail nil) (lp1 0) (t_ 0.0d0) (h 0.0d0) (g 0.0d0)
   (y 0.0d0) (x 0.0d0) (z 0.0d0) (l 0) (ll 0) (f 0.0d0) (i 0) (ii 0) (sn 0.0d0)
   (cs 0.0d0) (k 0) (kk 0) (j 0) (dnorm 0.0d0) (nqrs 0) (n10 0) (n 0)
   (two 0.0d0) (one 0.0d0) (zero 0.0d0)
  )
  (declare (type t wntv)) (declare (type t havers)) (declare (type t fail))
  (declare (type fixnum lp1)) (declare (type double-float t_))
  (declare (type double-float h)) (declare (type double-float g))
  (declare (type double-float y)) (declare (type double-float x))
  (declare (type double-float z)) (declare (type fixnum l))
  (declare (type fixnum ll)) (declare (type double-float f))
  (declare (type fixnum i)) (declare (type fixnum ii))
  (declare (type double-float sn)) (declare (type double-float cs))
  (declare (type fixnum k)) (declare (type fixnum kk))
  (declare (type fixnum j)) (declare (type double-float dnorm))
  (declare (type fixnum nqrs)) (declare (type fixnum n10))
  (declare (type fixnum n)) (declare (type double-float two))
  (declare (type double-float one)) (declare (type double-float zero))
  (comment "     subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc)")
  (comment
   "     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12"
  )
  (comment
   "     to appear in @solving least squares problems@, prentice-hall, 1974"
  )
  (comment "          qr algorithm for singular values of a bidiagonal matrix."
  )
  (comment "") (comment "     the bidiagonal matrix") (comment "")
  (comment "                       (q1,e2,0...    )")
  (comment "                       (   q2,e3,0... )")
  (comment "                d=     (       .      )")
  (comment "                       (         .   0)")
  (comment "                       (           .en)")
  (comment "                       (          0,qn)") (comment "")
  (comment "                 is pre and post multiplied by")
  (comment "                 elementary rotation matrices")
  (comment "                 ri and pi so that") (comment "")
  (comment "                 rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sn)")
  (comment "") (comment "                 to within working accuracy.")
  (comment "") (comment "  1. ei and qi occupy e(i) and q(i) as input.")
  (comment "") (comment "  2. rm...r1*c replaces @c@ in storage as output.")
  (comment "")
  (comment "  3. v*p1**(t)...pm**(t) replaces @v@ in storage as output.")
  (comment "") (comment "  4. si occupies q(i) as output.") (comment "")
  (comment "  5. the si@s are nonincreasing and nonnegative.") (comment "")
  (comment "     this code is based on the paper and @algol@ code..")
  (comment " ref..")
  (comment "  1. reinsch,c.h. and golub,g.h. @singular value decomposition")
  (comment "     and least squares solutions@ (numer. math.), vol. 14,(1970).")
  (comment "") (setf zero 0.0) (setf one 1.0) (setf two 2.0) (comment "")
  (setf n nn) (setf ipass 1) (if (<= n 0) (go end_label)) (setf n10 (* 10 n))
  (setf wntv (> nrv 0)) (setf havers (> ncc 0)) (setf fail false) (setf nqrs 0)
  (fset (fref e 1) zero) (setf dnorm zero)
  (fdo ((j 1 (+ j 1))) ((> j n) nil)
   (tagbody (setf dnorm (amax1 (+ (abs (fref q j)) (abs (fref e j))) dnorm)))
  )
  (fdo ((kk 1 (+ kk 1))) ((> kk n) nil)
   (tagbody (setf k (+ (+ n 1) (- kk))) (comment "")
    (comment "     test for splitting or rank deficiencies..")
    (comment
     "         first make test for last diagonal term, q(k), being small."
    )
    label20 (if (= k 1) (go label50))
    (arithmetic-if (diff (+ dnorm (fref q k)) dnorm) (go label50) (go label25)
     (go label50)
    )
    (comment "")
    (comment "     since q(k) is small we will make a special pass to")
    (comment "     transform e(k) to zero.") (comment "") label25
    (setf cs zero) (setf sn (- one))
    (fdo ((ii 2 (+ ii 1))) ((> ii k) nil)
     (tagbody (setf i (+ (+ k 1) (- ii)))
      (setf f (* (* -1 sn) (fref e (+ i 1))))
      (fset (fref e (+ i 1)) (* cs (fref e (+ i 1))))
      (multiple-value-setq (dummy_var f cs sn dummy_var)
       (g1 (fref q i) f cs sn (fref q i))
      )
      (comment "         transformation constructed to zero position (i,k).")
      (comment "") (if (not wntv) (go label40))
      (fdo ((j 1 (+ j 1))) ((> j nrv) nil)
       (tagbody
        (multiple-value-setq (cs sn dummy_var dummy_var)
         (g2 cs sn (fref v j i) (fref v j k))
      )))
      (comment "              accumulate rt. transformations in v.")
      (comment "") label40
    ))
    (comment "")
    (comment "         the matrix is now bidiagonal, and of lower order")
    (comment "         since e(k) .eq. zero..") (comment "") label50
    (fdo ((ll 1 (+ ll 1))) ((> ll k) nil)
     (tagbody (setf l (+ (+ k 1) (- ll)))
      (arithmetic-if (diff (+ dnorm (fref e l)) dnorm) (go label55)
       (go label100) (go label55)
      )
      label55
      (arithmetic-if (diff (+ dnorm (fref q (+ l (- 1)))) dnorm) (go label60)
       (go label70) (go label60)
      )
      label60
    ))
    (comment "     this loop can@t complete since e(1) = zero.") (comment "")
    (go label100) (comment "")
    (comment "         cancellation of e(l), l.gt.1.") label70 (setf cs zero)
    (setf sn (- one))
    (fdo ((i l (+ i 1))) ((> i k) nil)
     (tagbody (setf f (* (* -1 sn) (fref e i)))
      (fset (fref e i) (* cs (fref e i)))
      (arithmetic-if (diff (+ dnorm f) dnorm) (go label75) (go label100)
       (go label75)
      )
      label75
      (multiple-value-setq (dummy_var f cs sn dummy_var)
       (g1 (fref q i) f cs sn (fref q i))
      )
      (if (not havers) (go label90))
      (fdo ((j 1 (+ j 1))) ((> j ncc) nil)
       (tagbody
        (multiple-value-setq (cs sn dummy_var dummy_var)
         (g2 cs sn (fref c i j) (fref c (+ l (- 1)) j))
      )))
      label90
    ))
    (comment "") (comment "         test for convergence..") label100
    (setf z (fref q k)) (if (= l k) (go label170)) (comment "")
    (comment "         shift from bottom 2 by 2 minor of b**(t)*b.")
    (setf x (fref q l)) (setf y (fref q (+ k (- 1))))
    (setf g (fref e (+ k (- 1)))) (setf h (fref e k))
    (setf f
     (/ (+ (* (+ y (- z)) (+ y z)) (* (+ g (- h)) (+ g h))) (* (* two h) y))
    )
    (setf g (sqrt (+ one (expt f 2)))) (if (< f zero) (go label110))
    (setf t_ (+ f g)) (go label120) label110 (setf t_ (+ f (- g))) label120
    (setf f (/ (+ (* (+ x (- z)) (+ x z)) (* h (+ (/ y t_) (- h)))) x))
    (comment "") (comment "         next qr sweep..") (setf cs one)
    (setf sn one) (setf lp1 (+ l 1))
    (fdo ((i lp1 (+ i 1))) ((> i k) nil)
     (tagbody (setf g (fref e i)) (setf y (fref q i)) (setf h (* sn g))
      (setf g (* cs g))
      (multiple-value-setq (f h cs sn dummy_var)
       (g1 f h cs sn (fref e (+ i (- 1))))
      )
      (setf f (+ (* x cs) (* g sn))) (setf g (+ (* (* -1 x) sn) (* g cs)))
      (setf h (* y sn)) (setf y (* y cs)) (if (not wntv) (go label140))
      (comment "")
      (comment "              accumulate rotations (from the right) in @v@")
      (fdo ((j 1 (+ j 1))) ((> j nrv) nil)
       (tagbody
        (multiple-value-setq (cs sn dummy_var dummy_var)
         (g2 cs sn (fref v j (+ i (- 1))) (fref v j i))
      )))
      label140
      (multiple-value-setq (f h cs sn dummy_var)
       (g1 f h cs sn (fref q (+ i (- 1))))
      )
      (setf f (+ (* cs g) (* sn y))) (setf x (+ (* (* -1 sn) g) (* cs y)))
      (if (not havers) (go label160))
      (fdo ((j 1 (+ j 1))) ((> j ncc) nil)
       (tagbody
        (multiple-value-setq (cs sn dummy_var dummy_var)
         (g2 cs sn (fref c (+ i (- 1)) j) (fref c i j))
      )))
      (comment "              apply rotations from the left to")
      (comment "              right hand sides in @c@..") (comment "") label160
    ))
    (fset (fref e l) zero) (fset (fref e k) f) (fset (fref q k) x)
    (setf nqrs (+ nqrs 1)) (if (<= nqrs n10) (go label20))
    (comment "          return to @test for splitting@.") (comment "")
    (setf fail true) (comment "     ..")
    (comment "     cutoff for convergence failure. @nqrs@ will be 2*n usually."
    )
    label170 (if (>= z zero) (go label190)) (fset (fref q k) (- z))
    (if (not wntv) (go label190))
    (fdo ((j 1 (+ j 1))) ((> j nrv) nil)
     (tagbody (fset (fref v j k) (- (fref v j k))))
    )
    label190 (comment "         convergence. q(k) is made nonnegative..")
    (comment "")
  ))
  (if (= n 1) (go end_label))
  (fdo ((i 2 (+ i 1))) ((> i n) nil)
   (tagbody (if (> (fref q i) (fref q (+ i (- 1)))) (go label220)))
  )
  (if fail (setf ipass 2)) (go end_label) (comment "     ..")
  (comment "     every singular value is in order..") label220
  (fdo ((i 2 (+ i 1))) ((> i n) nil)
   (tagbody (setf t_ (fref q (+ i (- 1)))) (setf k (+ i (- 1)))
    (fdo ((j i (+ j 1))) ((> j n) nil)
     (tagbody (if (>= t_ (fref q j)) (go label230)) (setf t_ (fref q j))
      (setf k j) label230
    ))
    (if (= k (+ i (- 1))) (go label270)) (fset (fref q k) (fref q (+ i (- 1))))
    (fset (fref q (+ i (- 1))) t_) (if (not havers) (go label250))
    (fdo ((j 1 (+ j 1))) ((> j ncc) nil)
     (tagbody (setf t_ (fref c (+ i (- 1)) j))
      (fset (fref c (+ i (- 1)) j) (fref c k j)) (fset (fref c k j) t_)
    ))
    label250 (if (not wntv) (go label270))
    (fdo ((j 1 (+ j 1))) ((> j nrv) nil)
     (tagbody (setf t_ (fref v j (+ i (- 1))))
      (fset (fref v j (+ i (- 1))) (fref v j k)) (fset (fref v j k) t_)
    ))
    label270
  ))
  (comment "         end of ordering algorithm.") (comment "")
  (if fail (setf ipass 2)) (go end_label) end_label
  (return (values ipass q e nn v mdv nrv c mdc ncc))
))