***********************************************************************
 AUXILIARY Fortran 77 SUBROUTINES/FUNCTIONS called by DTLS,PSVD or PTLS.
           (not included in the BLAS or the LINPACK library)
------------------------------------------------------------------------

 SUBROUTINES/FUNCTIONS : BIDIAG, CANCEL, DAMIN, INIT, ESTIM, NSINGV,
                         QRQL, QRSTEP, QLSTEP, RESTOR, HOUSH, TR2.

      Sabine VAN HUFFEL
      ESAT Laboratory, KU Leuven.
      Kardinaal Mercierlaan 94, 3030 Heverlee, Belgium.

***********************************************************************

   SUBROUTINE : BIDIAG

 1 PURPOSE:

   The subroutine BIDIAG reduces an N by P matrix X by orthogonal
   Householder transformations G and H to bidiagonal form X = G' J H,
   with

            !Q(1) E(2)  0    ...   0  !
            ! 0   Q(2) E(3)        .  !
        J = ! .                    .  !     N >= P,
            ! .                   E(P)!
            ! 0   ...             Q(P)!

            !Q(1) E(2)  0    ...   0    0    !
            ! 0   Q(2) E(3)        .    .    !
        J = ! .                    .    .    !   N < P.
            ! .                   E(N)  0    !
            ! 0   ...             Q(N) E(N+1)!

 2 SPECIFICATION:

   SUBROUTINE BIDIAG(X, LDX, N, P, Q, E, WRK)
   INTEGER LDX, N, P
   DOUBLE PRECISION X(LDX,P), Q(*), E(*), WRK(N+P)

 3 ARGUMENT LIST:

   3.1 ARGUMENTS IN

       X - DOUBLE PRECISION array of DIMENSION (LDX,P)
           The leading N by P part of this array contains the N by P
           matrix X which is to be bidiagonalized.
           NOTE that this array is overwritten.

       LDX - INTEGER
           LDX is the leading dimension of the array X (LDX >= N).

       N - INTEGER
           N is the number of rows of the matrix X (N >= 1).

       P - INTEGER
           P is the number of columns of the matrix X (P >= 1).

   3.2 ARGUMENTS OUT

       X - DOUBLE PRECISION array of DIMENSION (LDX,P)
           The Householder transformations G and H are stored in the
           lower (resp. upper) triangular part of X.

       Q - DOUBLE PRECISION array of DIMENSION (min(N,P)).
           Q contains the diagonal elements of the computed bidiagonal
           matrix.

       E - DOUBLE PRECISION array of DIMENSION (min(N+1,P)).
           E contains the superdiagonal elements of the computed
           bidiagonal matrix. E(1) = 0.

   3.3 WORK SPACE

       WRK - DOUBLE PRECISION array of DIMENSION (N+P).

 5 EXTERNAL SUBROUTINES and FUNCTIONS:

   DSCAL, DNRM2, DDOT, DAXPY from BLAS.

 6 METHOD DESCRIPTION:

   BIDIAG is derived from the first part of the subroutine DSVDC of
   the LINPACK library. For more details, see [1, pp.11.6-11.9].

 7 REFERENCES:

   [1] J.J. Dongarra, J.R. Bunch, C.B. Moler and G.W. Stewart,
       LINPACK User's Guide. SIAM, Philadelphia (1979).

***********************************************************************

   1988, February 15.