Bidiagonalization of a matrix based on Lapack interface

The code computes the bidiagonal matrix decomposition of a matrix by calling Lapack routines.

This code provides a reliable bidiagonal decomposition routine based on Lapack subroutines. Lapack is one of most reliable fortran routines in numerical analysis.
In Lapack, the bidiagonalization is used to calculate SVD of a matrix via routines zgebrd and zungbr. See http://www.netlib.org/lapack/lug/node53.html But, Matlab does not support the bidiagonal routine as a standalone style.
Thus, this submission helps people to use bidiagonalization of a matrix, A = Q * B * P', where A is non-square and complex or real, B is bidiagonal and real, and Q and P are unitary or orthogonal.
The original Lapack interface routine comes from the work by Tim Toolan in "File Exchange".
This code can be used for understanding the Lapack interface.
This code is also compatible with Octave platform.
Note that the bidiagonalization is a good tool for matrix reduction preserving singular values. Instead, if you reduce a matrix preserving eigenvalues, you may want to use similarity transforms such as Hessenberg reduction (hess) or symmetric tridiagonalization.
These reductions are pre-calculated before Schur decomposition (schur), eigenvalue decomposition (eig), or singular values decomposition (svd) for speed.
The following is the function header for the usage example.
----
function [varargout]=bidiag_lapack(varargin)
% BIDIAG_LAPACK computes the matrix bidiagonal decomposition
% by calling LAPACK subroutines.
%
% [Q, B, P] = bidiag_lapack (A)
% [Q, B, P] = bidiag_lapack (A, ECON)
% R = bidiag_lapack (A)
% B = bidiag_lapack (A, ECON)
% [D, E] = bidiag_lapack (A)
% [Q, B, P, err] = bidiag_lapack (A, ...)
% [Q, B, P, err, R] = bidiag_lapack (A, ...)
%
% Compute the bidiagonal matrix decomposition of matrix A
% A = Q*B*P'
% where A is complex or real (not necessarily square),
% Q and P are unitary or orthogonal, and
% B is always real bidiagonal.
% Also, the signs of column vectors of Q are properly chosen to have
% the non-negative real diagonal entries of B.
%
% The bidiagonalization has applications in matrix reduction
% preserving singular values. Thus, the singular values of B are
% identical to those of A. (Note eigenvalues are not same.)
%
% ECON argument makes B square and Unitary matrices compact,
% where the unnecessary rows or columns of B, Q and P are eliminated.
%
% When called with one return value, it returns R, a compressed form of B
% and Householder reflectors of P and Q, which are computed in LAPACK
% subroutine zgebrd.
%
% When called with two return values, it returns D and E, which are
% the diagonal vector and the off-diagonal vector of B.
% B = diag(D) + diag(E,1) when m>=n;
% Otherwise, B = diag(D) + diag(E,-1).
%
% When called with 4 return values or more, it returns the residual error:
% err = norm(A-Q*B*P')/norm(A).
%
% See also: svd, hess, eig, qr.
% See LAPACK subroutines zgebrd and zungbr.