* Bug #13810 fixed - householder(v, k*v) returned column of %nan. Input parameters...
[scilab.git] / scilab / modules / linear_algebra / macros / householder.sci
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
14 function [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
109 endfunction