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