function?[U?B?V]?=?bidiagonalize（A）
%?Bidiagonalize?A??such?that?U?*?A?*?V?=?B?and?U‘?*?B?*?V‘?=?A.?B?is?the
%?same?size?as?A?has?no?vector?part?and?is?upper?or?lower?bidiagonal
%?depending?on?its?shape.?U?and?V?are?unitary?quaternion?matrices.

%?Copyright??2005?2008?Stephen?J.?Sangwine?and?Nicolas?Le?Bihan.
%?See?the?file?:?Copyright.m?for?further?details.

error（nargchk（1?1?nargin））?error（nargoutchk（3?3?nargout））

%?References:
%
%?Sangwine?S.?J.?and?Le?Bihan?N.
%?Quaternion?singular?value?decomposition?based?on?bidiagonalization
%?to?a?real?or?complex?matrix?using?quaternion?Householder?transformations
%?Applied?Mathematics?and?Computation?182（1）?1?November?2006?727-738?
%?DOI:10.1016/j.amc.2006.04.032.
%
%?Sangwine?S.?J.?and?Le?Bihan?N.
%?Quaternion?Singular?Value?Decomposition?based?on?Bidiagonalization
%?to?a?Real?Matrix?using?Quaternion?Householder?Transformations
%?arXiv:math.NA/0603251?10?March?2006?available?at?http://www.arxiv.org/
%
%?Gene?H.?Golub?and?Charles?van?Loan?‘Matrix?Computations‘?3rd?edition
%?Johns?Hopkins?University?Press?1996?section?5.1.4.?ISBN?0-8018-5414-8.
%?（This?describes?how?to?implement?a?Householder?transformation?without?an
%?explicit?Householder?matrix.）

[r?c]?=?size（A）;

if?prod（[r?c]）?==?1
????error（‘Cannot?bidiagonalize?a?matrix?of?one?element.‘）;????
end

if?c?<=?r
????[U?B?V]?=?internal_bidiagonalizer（A）;?%?Upper?bidiagonal?result.
else
????%?This?requires?a?lower?bidiagonal?result.?We?handle?this?by?operating
????%?on?the?Hermitian?transpose?of?A.?The?results?require?some?swapping
????%?and?transposition?to?get?the?correct?results.
????
????[V?B?U]?=?internal_bidiagonalizer（A‘）;?U?=?U‘;?B?=?B.‘;?V?=?V‘;
end

B?=?check（B）;?%?Verify?the?result?and?convert?to?exactly?bidiagonal?form
??????????????%?with?a?real?or?complex?B.

%?-------------------------------------------------------------------------

function?[U?B?V]?=?internal_bidiagonalizer（A）

[r?c]?=?size（A）;?assert（c?<=?r）;

%?Iterate?over?the?whole?matrix?dealing?with?one?column?and?one?row?on
%?each?pass?through?the?loop.?See?Figure?1?in?the?Applied?Mathematics?and
%?Computation?paper?cited?above?for?a?diagrammatic?representation?of?this
%?process.

U?=?eyeq（r）;?%?Initialise?the?three?matrix?results?to?appropriately?sized
B?=?A;???????%?quaternion?matrices.?These?are?conformant?so?that?U?*?B?*?V
V?=?eyeq（c）;?%?has?the?same?size?as?A.

for?i?=?1:c

????%?Compute?and?apply?a?left?Householder?transformation?to?the?ith
????%?column?（part?of?the?column?for?the?second?and?subsequent?columns）.
????
????[h?zeta]?=?householder_vector（B（i:end?i）?eye（r?-?i?+?1?1））;

????T?=?B（i:end?i:end）;?B（i:end?i:end）?=?（1./zeta）?.*?（T?-?h?*?（T‘?*?h）‘）;
????T?=?U（i:end??:???）;?U（i:end??:???）?=?（1./zeta）?.*?（T?-?h?*?（T‘?*?h）‘）;
????
????if?i?==?c?return;?end?%?On?the?last?column?we?are?done?since?there
???????????????????????????%?is?no?corresponding?row.?See?Figure?1?in?the
???????????????????????????%?Applied?Mathematics?