INTERFACE SVD;

PROCEDURE Bidiagonalize(
    VAR A: ARRAY OF ARRAY OF REAL;
    VAR D: ARRAY OF LONGREAL;
    VAR E: ARRAY OF LONGREAL;
  );
  (*
    Converts "A" to a bidiagonal matrix "B".
    
    Let "A" have "M" rows by "N" columns, and let "L" be "MIN(M,N)".
    The resulting matrix is lower bidiagonal if "M<=N", and upper
    bidiagonal if "M>N".  The diagonal elements will be returned IN
    "D[0..L-1]", and the paradiagonal ones will be in "E[0..L-1]". 
    
    The transformations that take "A" to "B" are a sequence of
    Householder rotations, acting alternatively on the row vectors and
    column vectors of "A".  A rotation that acts on the row vectors is
    said to have "direction [0,1]".  It will modify only columns
    "p..N-1" for some index "p" (the "pivot index"), and has the efect
    of nulling out all elements "A[r, p+1..N-1] for some row "r" (the 
    "reference index"). The element "A[r,p]" is called the "pivot"
    of the operation, and elements "A[r, p+1..N-1]" to the right
    of it are the "tail".  
    
    A rotation that acts on the column vectors is said to have
    "direction [1,0]", and has a completely symmetrical effect: it
    modifies only rows "p..M-1" for some "pivot index" "p", and anihilates
    elements "A[p+1..M-1, r]" for some "reference index" index "r".
    Element "A[p,r]" is called the "pivot", and the elements
    "A[p+1..M-1, r]" below it are the "tail".
    
    The first rotation has pivot element "[0,0]", and its
    direction points towards the longest dimension of "A"; namely,
    "[0,1]" if "M<=N", and "[1,0]" if "M>N".  If a rotation has 
    pivot "A[i,j]" and direction "[di,dj]", then the next one has 
    pivot "A[i+dj, j+di]" and direction "[dj,di]".
    
    Each rotation consists of two reflections.  The first reflection
    replaces each affected vector "u" (row or column) by the vector "v
    = u - (u.h) h", were the vector "h" is a parameter that defines
    the rotation.  Elements "h[0..p-1]" are all zero, so "v[0..p-1] =
    u[0..p-1]".  The second reflection merely negates "v[p]".
    
    The vectors "h" are returned in the array "A", which is therefore
    destroyed in the process.  For each rotation, "h[p]" 
    will be stored in the pivot element, and "h[p+1], h[p+2], ..."
    will be stored in the tail elements of that rotation.
  *)


END SVD.