1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 2// Copyright (C) 2008 - INRIA 3// Copyright (C) 2015 - Samuel GOUGEON : http://bugzilla.scilab.org/13810 4// 5// Copyright (C) 2012 - 2016 - Scilab Enterprises 6// 7// This file is hereby licensed under the terms of the GNU GPL v2.0, 8// pursuant to article 5.3.4 of the CeCILL v.2.1. 9// This file was originally licensed under the terms of the CeCILL v2.1, 10// and continues to be available under such terms. 11// For more information, see the COPYING file which you should have received 12// along with this program. 13 14function [u, hm] = householder(v,w) 15 16 fname = "householder" 17 [lhs,rhs] = argn(0) 18 19 // Example = demo 20 // -------------- 21 if rhs==0 then 22 disp("householder() example: Reflect an object using the Householder matrix") 23 [funs, path] = libraryinfo(whereis("householder")); 24 editor(path+"householder.sce") 25 exec(path+"householder.sce", -1) 26 return 27 end 28 29 // CHECKING INPUT PARAMETERS 30 // ------------------------- 31 if rhs<2 then 32 w = eye(v) 33 end 34 if typeof(v(:))~="constant" 35 msg = _("%s: Wrong type for argument %d: Decimal or complex numbers expected.\n") 36 error(msprintf(msg, fname, 1)) 37 end 38 if typeof(w(:))~="constant" 39 msg = _("%s: Wrong type for argument %d: Decimal or complex numbers expected.\n") 40 error(msprintf(msg, fname, 2)) 41 end 42 if length(find(size(v)>1))>1 then 43 msg = _("%s: Wrong size for input argument #%d: Column vector expected.\n") 44 error(msprintf(msg, fname, 1)) 45 end 46 if length(find(size(w)>1))>1 then 47 msg = _("%s: Wrong size for input argument #%d: Column vector expected.\n") 48 error(msprintf(msg, fname, 2)) 49 end 50 if length(v)~=length(w) then 51 msg = _("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n") 52 error(msprintf(msg, fname, 1, 2)) 53 end 54 v = v(:) 55 w = w(:) 56 57 // PROCESSING 58 // ---------- 59 if v==[] | w==[] then 60 u = [] 61 hm = [] 62 return 63 end 64 a = sqrt((v'*v) / (w'*w)) 65 66 // Are v and w colinear? 67 colinear = %t 68 kn = find(w==0) 69 if kn~=[] then 70 colinear = colinear & and(v(kn)==0) 71 end 72 k = find(w~=0) 73 if colinear & k~=[] then 74 tmp = v(k)./w(k) 75 colinear = colinear & and(tmp==tmp(1)) 76 end 77 78 // v and w coliear and reals 79 if colinear & isreal(v) & isreal(w) then 80 if tmp(1)<0 // v and w are opposite 81 u = v 82 else // v and w in the same direction 83 // => u is in the hyperplane orthogonal to v, and then u'*v==0 84 u = zeros(v) 85 if kn~=[] 86 u(kn(1)) = 1 87 else 88 if length(u)>1 89 u(1:2) = [1 ; -w(1)/w(2)] 90 else 91 u = %nan 92 end 93 end 94 end 95 else 96 // v and w are not colinear or are complex 97 if or(v~=a*w) // v and w not colinear 98 u = v - a*w 99 else 100 u = v + a*w 101 end 102 end 103 try 104 u = u / norm(u) 105 end 106 if lhs>1 then 107 hm = eye() - 2*u*u' 108 end 109endfunction 110