X-Git-Url: http://gitweb.scilab.org/?p=scilab.git;a=blobdiff_plain;f=scilab%2Fmodules%2Flinear_algebra%2Fmacros%2Fhouseholder.sci;h=d480506c6748578a6f62f8debaeb4d6d45263974;hp=9d0f5447e3f9747922afb8357842328cc7cc1f6c;hb=f0ad1547b5fd8108c7ccc0779d19216160dfdacf;hpb=ed00e2c731e6ced1bfb964b4ca30f11020fc3320 diff --git a/scilab/modules/linear_algebra/macros/householder.sci b/scilab/modules/linear_algebra/macros/householder.sci index 9d0f544..d480506 100644 --- a/scilab/modules/linear_algebra/macros/householder.sci +++ b/scilab/modules/linear_algebra/macros/householder.sci @@ -1,5 +1,6 @@ // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab // Copyright (C) 2008 - INRIA +// Copyright (C) 2015 - Samuel GOUGEON : http://bugzilla.scilab.org/13810 // // Copyright (C) 2012 - 2016 - Scilab Enterprises // @@ -10,19 +11,99 @@ // For more information, see the COPYING file which you should have received // along with this program. -function u=householder(v,w) - //Syntax - //u=householder(v [,w]) - //Description - //given 2 column vectors v w of same size householder(v,w) returns a unitary - //column vector u, such that (eye-2*u*u')*v is proportional to w. - //(eye-2*u*u') is the orthogonal Householder reflexion matrix - // - // w default value is eye(v). In this case vector (eye-2*u*u')*v is the - // vector eye(v)*(+-norm(v)) - [lhs,rhs]=argn(0) - if rhs<2 then w=eye(v),end - a=-sqrt((v'*v)/(w'*w)) - u=v+a*w - u=u/norm(u) +function [u, hm] = householder(v,w) + + fname = "householder" + [lhs,rhs] = argn(0) + + // Example = demo + // -------------- + if rhs==0 then + disp("householder() example: Reflect an object using the Householder matrix") + [funs, path] = libraryinfo(whereis("householder")); + editor(path+"householder.sce") + exec(path+"householder.sce", -1) + return + end + + // CHECKING INPUT PARAMETERS + // ------------------------- + if rhs<2 then + w = eye(v) + end + if typeof(v(:))~="constant" + msg = _("%s: Wrong type for argument %d: Decimal or complex numbers expected.\n") + error(msprintf(msg, fname, 1)) + end + if typeof(w(:))~="constant" + msg = _("%s: Wrong type for argument %d: Decimal or complex numbers expected.\n") + error(msprintf(msg, fname, 2)) + end + if length(find(size(v)>1))>1 then + msg = _("%s: Wrong size for input argument #%d: Column vector expected.\n") + error(msprintf(msg, fname, 1)) + end + if length(find(size(w)>1))>1 then + msg = _("%s: Wrong size for input argument #%d: Column vector expected.\n") + error(msprintf(msg, fname, 2)) + end + if length(v)~=length(w) then + msg = _("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n") + error(msprintf(msg, fname, 1, 2)) + end + v = v(:) + w = w(:) + + // PROCESSING + // ---------- + if v==[] | w==[] then + u = [] + hm = [] + return + end + a = sqrt((v'*v) / (w'*w)) + + // Are v and w colinear? + colinear = %t + kn = find(w==0) + if kn~=[] then + colinear = colinear & and(v(kn)==0) + end + k = find(w~=0) + if colinear & k~=[] then + tmp = v(k)./w(k) + colinear = colinear & and(tmp==tmp(1)) + end + + // v and w coliear and reals + if colinear & isreal(v) & isreal(w) then + if tmp(1)<0 // v and w are opposite + u = v + else // v and w in the same direction + // => u is in the hyperplane orthogonal to v, and then u'*v==0 + u = zeros(v) + if kn~=[] + u(kn(1)) = 1 + else + if length(u)>1 + u(1:2) = [1 ; -w(1)/w(2)] + else + u = %nan + end + end + end + else + // v and w are not colinear or are complex + if or(v~=a*w) // v and w not colinear + u = v - a*w + else + u = v + a*w + end + end + try + u = u / norm(u) + end + if lhs>1 then + hm = eye() - 2*u*u' + end endfunction