* Bug #13810 fixed - householder(v, k*v) returned column of %nan. Input parameters... 49/17849/6
Samuel GOUGEON [Thu, 3 Mar 2016 15:44:16 +0000 (16:44 +0100)]
Change-Id: Ief87ef84ecba005b2bd3def5f4b5c940d29ccb4e

scilab/CHANGES
scilab/modules/helptools/images/householder_1.png [new file with mode: 0644]
scilab/modules/linear_algebra/help/en_US/factorization/householder.xml
scilab/modules/linear_algebra/help/fr_FR/factorization/householder.xml
scilab/modules/linear_algebra/help/ja_JP/factorization/householder.xml
scilab/modules/linear_algebra/help/pt_BR/factorization/householder.xml
scilab/modules/linear_algebra/macros/householder.sce [new file with mode: 0644]
scilab/modules/linear_algebra/macros/householder.sci
scilab/modules/linear_algebra/tests/nonreg_tests/bug_13810.dia.ref [new file with mode: 0644]
scilab/modules/linear_algebra/tests/nonreg_tests/bug_13810.tst [new file with mode: 0644]

index aebd6f4..4545c20 100644 (file)
@@ -302,6 +302,8 @@ In 6.0.0:
 
 * Bug #9456 fixed  - bench_run did not work on a path or in a toolbox
 
+* Bug #13810  fixed - householder(v, k*v) returned column of %nan. Input parameters were not checked. The Householder matrix could not be returned. Help pages were inaccurate and without examples. There was no householder() demo.
+
 * Bug #13869 fixed - bench_run with option nb_run=10 did not override the NB RUN tags
 
 * Bug #14035 fixed - ndgrid did not manage all homogeneous data type (booleans, integers, polynomials, rationals, strings, [])
diff --git a/scilab/modules/helptools/images/householder_1.png b/scilab/modules/helptools/images/householder_1.png
new file mode 100644 (file)
index 0000000..7111f56
Binary files /dev/null and b/scilab/modules/helptools/images/householder_1.png differ
index 131bd99..5a9fde0 100644 (file)
@@ -2,7 +2,8 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
- * 
+ * Copyright (C) 2015 - Samuel GOUGEON
+ *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
  *
  * This file is hereby licensed under the terms of the GNU GPL v2.0,
@@ -11,7 +12,6 @@
  * and continues to be available under such terms.
  * For more information, see the COPYING file which you should have received
  * along with this program.
- *
  -->
 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en" xml:id="householder">
     <refnamediv>
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
-        <synopsis>u=householder(v [,w])</synopsis>
+        <synopsis>         householder         // demo</synopsis>
+        <synopsis> u     = householder(v [,w])</synopsis>
+        <synopsis>[u, H] = householder(v [,w])</synopsis>
     </refsynopsisdiv>
-    <refsection>
+    <refsection role="arguments">
         <title>Arguments</title>
         <variablelist>
             <varlistentry>
                 <term>w</term>
                 <listitem>
                     <para>
-                        real or complex column vector with same size as <literal>v</literal>. Default value is <literal>eye(v)</literal>
+                        real or complex column vector with same size as <literal>v</literal>.
+                        Default value is <literal>eye(v)</literal> ((Ox) axis).
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>u</term>
                 <listitem>
-                    <para>real or complex column vector</para>
+                    <para>
+                        unit vector lying in the <literal>(v,w)</literal> plane and orthogonal
+                        to the bisectrix of <literal>(v,w)</literal>.
+                        Column of size(v) of real or complex numbers.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>H</term>
+                <listitem>
+                    <para>
+                        Orthogonal Householder reflexion matrix: <literal>H= eye() - 2*u*u'</literal>.
+                        <varname>H</varname> is such that <literal>inv(H)==H</literal>,
+                        <literal>H'==H</literal>, and <literal>det(H)==-1</literal>.
+                    </para>
+                    <para>
+                        If <varname>v</varname> and <varname>w</varname> are real,
+                        <literal>H*v</literal> is proportional to <varname>w</varname>.
+                    </para>
                 </listitem>
             </varlistentry>
         </variablelist>
     </refsection>
-    <refsection>
+    <refsection role="description">
         <title>Description</title>
         <para>
-            given 2 column vectors <literal>v</literal>, <literal> w</literal> of same size, <literal>householder(v,w)</literal> returns a unitary 
-            column vector <literal>u</literal>, such that <literal> (eye()-2*u*u')*v</literal> is proportional to <literal>w</literal>.
-            <literal>(eye()-2*u*u')</literal> is the orthogonal Householder reflexion matrix .
+            <literal>householder(..)</literal> computes the unit vector <varname>u</varname>
+            lying in the <literal>(v,w)</literal> plane and orthogonal to the bisectrix of
+            <literal>(v,w)</literal>.
+        </para>
+        <para>
+            If <varname>v</varname> and <varname>w</varname> are proportional:
+            <itemizedlist>
+                <listitem>
+                    <para>
+                        If they are opposite, <literal>u= v/|v|</literal> is returned.
+                    </para>
+                </listitem>
+                <listitem>
+                    If they are real and have the same direction, <varname>u</varname> is set
+                    in the (xOy) plane with a priori <literal>u(1)>0</literal>, and orthogonal to
+                    <varname>v</varname> (<literal>u'*v==0</literal>). However,
+                    <itemizedlist>
+                        <listitem>
+                            If they are along (Ox), <literal>u = (Oy+)</literal> is returned instead.
+                        </listitem>
+                        <listitem>
+                            If <varname>v</varname> and <varname>w</varname> are scalars with same
+                            signs, the orthogonal sub-space is restricted to <literal>{0}</literal>
+                            that can't be normalized:
+                            <varname>u</varname> and <varname>H</varname> are then set to
+                            <literal>%nan</literal>.
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </itemizedlist>
         </para>
         <para>
-            <literal>w</literal> default value is <literal> eye(v)</literal>. In this case vector <literal> (eye()-2*u*u')*v</literal> is the 
-            vector  <literal> eye(v)*norm(v)</literal>.
+            If the related reflexion matrix <varname>H</varname> is computed, for any point A
+            of column coordinates <literal>a</literal>, <literal>H*a</literal> are the coordinates of
+            the reflected image of A on P.
         </para>
+        <note>
+            If <varname>v</varname> or/and <varname>w</varname> are in row, they are priorly
+            transposed into columns.
+        </note>
+        <warning>
+            If <varname>v</varname> or/and <varname>w</varname> are <literal>[]</literal>,
+            <literal>[]</literal> is returned for <varname>u</varname> and <varname>H</varname>.
+        </warning>
+    </refsection>
+    <refsection role="examples">
+        <title>Examples</title>
+        <programlisting role="example"><![CDATA[
+a = [ rand(1,1) 0  0 ]';
+[ra hm] = householder(a);
+[a ra hm*a ]
+norm(ra)
+
+b = rand(3,1);
+[rb, hm] = householder(b);
+[b rb eye(b) clean(hm*b) ]
+norm(rb)
+
+[rb2b, hm] = householder(b, 2*b);
+[b rb2b clean(hm*b ./ b) ]  // last column must be uniform
+norm(rb2b)                  // must be 1
+
+c = rand(3,1);
+[rbc, hm] = householder(b,c);
+norm(rbc)          // must be 1
+hm*b ./c           // must be uniform
+
+d = b + %i*c;
+e = rand(3,1) + %i*rand(3,1);
+[rde, hm] = householder(d,e);
+norm(rbc)               // must be 1
+clean(inv(hm) - hm)     // must be zeros(3,3)
+clean(hm' - hm)         // must be zeros(3,3)
+clean(det(hm))          // must be -1
+ ]]></programlisting>
+        
+        <para>APPLICATION : Reflected image of an object</para>
+        <programlisting role="example"><![CDATA[
+// (OA) = [0 0 1] is reflected in O into (OB) = [ 1 1 0.3 ]:
+[n, H] = householder([0 0 1]', [ 1 1 0.3 ]');
+// "n" is the unit vector orthogonal to the reflecting plane
+
+// Emitting object (feature from shell demo):
+u = linspace(0,2*%pi,40);
+v = linspace(0,2*%pi,20);
+Xe = (cos(u).*u)'*(1+cos(v)/2)+10;
+Ye = (u/2)'*sin(v);
+Ze = (sin(u).*u)'*(1+cos(v)/2);
+
+// Reflected object:
+Pe = [ Xe(:)' ; Ye(:)' ; Ze(:)'];
+Pr = H*Pe;
+Xr = matrix(Pr(1,:),40,-1);
+Yr = matrix(Pr(2,:),40,-1);
+Zr = matrix(Pr(3,:),40,-1);
+
+// Reflecting plane containing O: n(1).x + n(2).y + n(3).z = 0
+//   Sampling space:
+x = linspace(min([Xe(:);Xr(:)]), max([Xe(:);Xr(:)]),20);
+y = linspace(min([Ye(:);Yr(:)]), max([Ye(:);Yr(:)]),20);
+[X, Y] = meshgrid(x,y);
+//   Generating the mirror:
+deff("z = mirror(x,y,n)","z = -n(1)/n(3)*x - n(2)/n(3)*y")
+Zm = mirror(X,Y,n);
+
+// Plotting:
+clf
+drawlater
+f = gcf();
+f.color_map = [ 0.8 0.8 0.8 ; jetcolormap(100)];
+surf(Xe,Ye,Ze)
+surf(X,Y,Zm)
+surf(Xr,Yr,Zr)
+a = gca();
+a.isoview = "on";
+a.rotation_angles = [74 123];
+a.children.color_flag = 0;
+a.children.color_mode = 0;
+a.children(1).foreground = color("red");
+a.children(2).foreground = 1;
+a.children(3).foreground = color("green");
+drawnow
+ ]]></programlisting>
+        <scilab:image>
+            householder();
+        </scilab:image>
     </refsection>
     <refsection role="see also">
         <title>See Also</title>
             </member>
         </simplelist>
     </refsection>
+    <refsection role="history">
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.0</revnumber>
+                <revdescription>
+                    <para>Householder reflexion matrix added as second output parameter.
+                        Demo householder() added. Help page reviewed.
+                    </para>
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
index fffe07c..6158e2f 100644 (file)
@@ -2,7 +2,8 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
- * 
+ * Copyright (C) 2015 - Samuel GOUGEON
+ *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
  *
  * This file is hereby licensed under the terms of the GNU GPL v2.0,
@@ -11,7 +12,6 @@
  * and continues to be available under such terms.
  * For more information, see the COPYING file which you should have received
  * along with this program.
- *
  -->
 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="fr" xml:id="householder">
     <refnamediv>
     </refnamediv>
     <refsynopsisdiv>
         <title>Séquence d'appel</title>
-        <synopsis>u=householder(v [,w])</synopsis>
+        <synopsis>         householder         // demo</synopsis>
+        <synopsis> u     = householder(v [,w])</synopsis>
+        <synopsis>[u, H] = householder(v [,w])</synopsis>
     </refsynopsisdiv>
-    <refsection>
+    <refsection role="arguments">
         <title>Paramètres</title>
         <variablelist>
             <varlistentry>
                 <term>w  </term>
                 <listitem>
                     <para>
-                        vecteur colonne réel ou complexe de même taille que <literal>v</literal> (la valeur par défaut est <literal>eye(v)</literal>)
+                        vecteur colonne réel ou complexe de même taille que <literal>v</literal>
+                        La valeur par défaut est <literal>eye(v)</literal> (axe (Ox)).
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>u  </term>
                 <listitem>
-                    <para>vecteur colonne réel ou complexe
+                    <para>
+                        vecteur unitaire résidant dans le plan <literal>(v,w)</literal> et orthogonal
+                        à la bissectrice de <literal>(v,w)</literal>.
+                        Colonne de nombres réels ou complexes, de taille size(v).
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>H</term>
+                <listitem>
+                    <para>
+                        Matrice orthogonale de réflexion de Householder : <literal>H= eye() - 2*u*u'</literal>.
+                        <varname>H</varname> est telle que <literal>inv(H)==H</literal>,
+                        <literal>H'==H</literal>, et <literal>det(H)==-1</literal>.
+                    </para>
+                    <para>
+                        Si <varname>v</varname> et <varname>w</varname> sont des vecteurs réels,
+                        <literal>H*v</literal> est colinéaire à <varname>w</varname>.
                     </para>
                 </listitem>
             </varlistentry>
         </variablelist>
     </refsection>
     <refsection>
-        <title>Description</title>
+        <title role="description">Description</title>
         <para>
-            Etant donnés deux vecteurs colonnes <literal>v</literal> et <literal>w</literal> de même taille, <literal>householder(v,w)</literal> renvoie un vecteur normé <literal>u</literal>, tel que 
-            <literal>(eye()-2*u*u')*v</literal> est colinéaire à <literal>w</literal>.
-            <literal>(eye()-2*u*u')</literal> est la matrice de la transformation de Householder correspondante.
+            <literal>householder(..)</literal> calcule le vecteur unitaire <varname>u</varname>
+            résidant dans le plan <literal>(v,w)</literal> et orthogonal à la bissectrice de
+            <literal>(v,w)</literal>.
         </para>
         <para>
-            La valeur par défaut de <literal>w</literal> est <literal> eye(v)</literal>. Dans ce cas le vecteur <literal> (eye()-2*u*u')*v</literal> est égal à <literal> eye(v)*norm(v)</literal>.
+            Si <varname>v</varname> et <varname>w</varname> sont proportionels :
+            <itemizedlist>
+                <listitem>
+                    S'ils sont opposés, <literal>u= v/|v|</literal> est retourné.
+                </listitem>
+                <listitem>
+                    S'ils sont à coordonnées réelles et sont dans la même direction,
+                    <varname>u</varname> est choisi dans le
+                    plan (xOy) avec a priori <literal>u(1)>0</literal>, et orthogonal à
+                    <varname>v</varname> (<literal>u'*v==0</literal>). Cependant,
+                    <itemizedlist>
+                        <listitem>
+                            S'ils sont selon (Ox), <literal>u = (Oy+)</literal> est retourné.
+                        </listitem>
+                        <listitem>
+                            S'ils sont scalaires (et de mêmes signes),
+                            le sous-espace orthogonal est réduit à <literal>{0}</literal> qui n'est
+                            pas normalisable :
+                            <varname>u</varname> et <varname>H</varname> sont mis à <literal>%nan</literal>.
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </itemizedlist>
         </para>
+        
+        <para>
+            Si la matrice de réflexion <varname>H</varname> correspondante est calculée,
+            pour tout point A de coordonnées <literal>a</literal> en colonne, <literal>H*a</literal>
+            sont les coordonnées de l'image de A réfléchie sur P.
+        </para>
+        <note>
+            Si <varname>v</varname> ou/et <varname>w</varname> sont des vecteurs ligne,
+            ils sont préalablement transposés en colonnes.
+        </note>
+        <warning>
+            Si <varname>v</varname> ou/et <varname>w</varname> sont <literal>[]</literal>,
+            <varname>u</varname> et <varname>H</varname> valent alors <literal>[]</literal>.
+        </warning>
+    </refsection>
+    <refsection role="examples">
+        <title>Exemples</title>
+        <programlisting role="example"><![CDATA[
+a = [ rand(1,1) 0  0 ]';
+[ra hm] = householder(a);
+[a ra hm*a ]
+norm(ra)
+
+b = rand(3,1);
+[rb, hm] = householder(b);
+[b rb eye(b) clean(hm*b) ]
+norm(rb)
+
+[rb2b, hm] = householder(b, 2*b);
+[b rb2b clean(hm*b ./ b) ]  // last column must be uniform
+norm(rb2b)                  // must be 1
+
+c = rand(3,1);
+[rbc, hm] = householder(b,c);
+norm(rbc)          // must be 1
+hm*b ./c           // must be uniform
+
+d = b + %i*c;
+e = rand(3,1) + %i*rand(3,1);
+[rde, hm] = householder(d,e);
+norm(rbc)               // must be 1
+clean(inv(hm) - hm)     // must be zeros(3,3)
+clean(hm' - hm)         // must be zeros(3,3)
+clean(det(hm))          // must be -1
+ ]]></programlisting>
+        
+        <para>APPLICATION : Image réfléchie d'un objet</para>
+        <programlisting role="example"><![CDATA[
+// (OA) = [0 0 1] se réfléchit en O en (OB) = [ 1 1 0.3 ]:
+[n, H] = householder([0 0 1]', [ 1 1 0.3 ]');
+// "n" est la normale au plan réflecteur
+
+// Objet source :
+u = linspace(0,2*%pi,40);
+v = linspace(0,2*%pi,20);
+Xe = (cos(u).*u)'*(1+cos(v)/2)+10;
+Ye = (u/2)'*sin(v);
+Ze = (sin(u).*u)'*(1+cos(v)/2);
+
+// Image réfléchie de l'object :
+P = [ Xe(:)' ; Ye(:)' ; Ze(:)'];
+Pr = H*P;
+Xr = matrix(Pr(1,:),40,-1);
+Yr = matrix(Pr(2,:),40,-1);
+Zr = matrix(Pr(3,:),40,-1);
+
+// Plan réflecteur contnant O : n(1).x + n(2).y + n(3).z = 0:
+//  Grille d'espace :
+x = linspace(min([Xe(:);Xr(:)]), max([Xe(:);Xr(:)]),20);
+y = linspace(min([Ye(:);Yr(:)]), max([Ye(:);Yr(:)]),20);
+[X, Y] = meshgrid(x,y);
+//  Création du plan:
+deff("z = mirror(x,y,n)","z = -n(1)/n(3)*x - n(2)/n(3)*y")
+Zm = mirror(X,Y,n);
+
+// Illustration graphique :
+clf
+drawlater
+f = gcf();
+f.color_map = [ 0.8 0.8 0.8 ; jetcolormap(100)];
+surf(Xe,Ye,Ze)
+surf(X,Y,Zm)
+surf(Xr,Yr,Zr)
+a = gca();
+a.isoview = "on";
+a.rotation_angles = [74 123];
+a.children.color_flag = 0;
+a.children.color_mode = 0;
+a.children(1).foreground = color("red");
+a.children(2).foreground = 1;
+a.children(3).foreground = color("green");
+drawnow
+ ]]></programlisting>
+        <scilab:image>
+            householder();
+        </scilab:image>
     </refsection>
     <refsection role="see also">
         <title>Voir aussi</title>
             </member>
         </simplelist>
     </refsection>
+    <refsection role="history">
+        <title>Historique</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.0</revnumber>
+                <revdescription>
+                    <para>Matrice de réflexion de Householder fournie en second paramètre de sortie.
+                        Demo householder() ajoutée. Page d'aide revue.
+                    </para>
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
index 6f8a22e..94c17ed 100644 (file)
@@ -1,9 +1,9 @@
 <?xml version="1.0" encoding="UTF-8"?>
-
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
- * 
+ * Copyright (C) 2015 - Samuel GOUGEON
+ *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
  *
  * This file is hereby licensed under the terms of the GNU GPL v2.0,
  * and continues to be available under such terms.
  * For more information, see the COPYING file which you should have received
  * along with this program.
- *
  -->
 
 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="ja" xml:id="householder">
     
     <refnamediv>
-        
         <refname>householder</refname>
-        
         <refpurpose>ハウスホルダー直交鏡映行列</refpurpose>
-        
     </refnamediv>
-    
     <refsynopsisdiv>
-        
         <title>呼び出し手順</title>
-        
-        <synopsis>u=householder(v [,w])</synopsis>
-        
+        <synopsis>         householder         // demo</synopsis>
+        <synopsis> u     = householder(v [,w])</synopsis>
+        <synopsis>[u, H] = householder(v [,w])</synopsis>
     </refsynopsisdiv>
-    
-    <refsection>
-        
+    <refsection role="arguments">
         <title>引数</title>
-        
         <variablelist>
-            
             <varlistentry>
-                
                 <term>v</term>
-                
                 <listitem>
-                    
                     <para>実数または複素数の列ベクトル</para>
-                    
                 </listitem>
-                
             </varlistentry>
-            
             <varlistentry>
-                
                 <term>w</term>
-                
                 <listitem>
-                    
                     <para>
-                        
                         <literal>v</literal>と同じ大きさの実数または複素数の列ベクトル.
-                        
-                        デフォルト値は<literal>eye(v)</literal>
-                        
+                        デフォルト値は<literal>eye(v)</literal>((Ox) axis).
                     </para>
-                    
                 </listitem>
-                
             </varlistentry>
-            
             <varlistentry>
-                
                 <term>u</term>
-                
                 <listitem>
-                    
-                    <para>実数または複素数の列ベクトル</para>
-                    
+                    <para>
+                        unit vector lying in the <literal>(v,w)</literal> plane and orthogonal
+                        to the bisectrix of <literal>(v,w)</literal>.
+                        Column of size(v) of real or complex numbers.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>H</term>
+                <listitem>
+                    <para>
+                        Orthogonal Householder reflexion matrix: <literal>H= eye() - 2*u*u'</literal>.
+                        <varname>H</varname> is such that <literal>inv(H)==H</literal>,
+                        <literal>H'==H</literal>, and <literal>det(H)==-1</literal>.
+                    </para>
+                    <para>
+                        If <varname>v</varname> and <varname>w</varname> are real,
+                        <literal>H*v</literal> is proportional to <varname>w</varname>.
+                    </para>
                 </listitem>
-                
             </varlistentry>
-            
         </variablelist>
-        
     </refsection>
-    
-    <refsection>
-        
+    <refsection role="description">
         <title>説明</title>
-        
         <para>
-            
-            同じ大きさの列ベクトル
-            
-            <literal>v</literal>, <literal> w</literal> を指定すると, 
-            
-            <literal>householder(v,w)</literal> は,
-            
-            <literal> (eye()-2*u*u')*v</literal>が<literal>w</literal>に比例するような
-            
-            ユニタリ列ベクトル<literal>u</literal>を返します.
-            
-            <literal>(eye()-2*u*u')</literal> はハウスホルダー直交鏡映行列です.
-            
+            <literal>householder(..)</literal> computes the unit vector <varname>u</varname>
+            lying in the <literal>(v,w)</literal> plane and orthogonal to the bisectrix of
+            <literal>(v,w)</literal>.
         </para>
-        
         <para>
-            
-            <literal>w</literal> のデフォルト値は <literal> eye(v)</literal>です. 
-            
-            この場合,ベクトル<literal> (eye()-2*u*u')*v</literal> はベクトル 
-            
-            <literal> eye(v)*norm(v)</literal>です.
-            
+            If <varname>v</varname> and <varname>w</varname> are proportional:
+            <itemizedlist>
+                <listitem>
+                    <para>
+                        If they are opposite, <literal>u= v/|v|</literal> is returned.
+                    </para>
+                </listitem>
+                <listitem>
+                    If they are real and have the same direction, <varname>u</varname> is set
+                    in the (xOy) plane with a priori <literal>u(1)>0</literal>, and orthogonal to
+                    <varname>v</varname> (<literal>u'*v==0</literal>). However,
+                    <itemizedlist>
+                        <listitem>
+                            If they are along (Ox), <literal>u = (Oy+)</literal> is returned instead.
+                        </listitem>
+                        <listitem>
+                            If <varname>v</varname> and <varname>w</varname> are scalars with same
+                            signs, the orthogonal sub-space is restricted to <literal>{0}</literal>
+                            that can't be normalized:
+                            <varname>u</varname> and <varname>H</varname> are then set to
+                            <literal>%nan</literal>.
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </itemizedlist>
         </para>
         
+        <para>
+            If the related reflexion matrix <varname>H</varname> is computed, for any point A
+            of column coordinates <literal>a</literal>, <literal>H*a</literal> are the coordinates of
+            the reflected image of A on P.
+        </para>
+        <note>
+            If <varname>v</varname> or/and <varname>w</varname> are in row, they are priorly
+            transposed into columns.
+        </note>
+        <warning>
+            If <varname>v</varname> or/and <varname>w</varname> are <literal>[]</literal>,
+            <literal>[]</literal> is returned for <varname>u</varname> and <varname>H</varname>.
+        </warning>
+    </refsection>
+    <refsection role="examples">
+        <title>例</title>
+        <programlisting role="example"><![CDATA[
+a = [ rand(1,1) 0  0 ]';
+[ra hm] = householder(a);
+[a ra hm*a ]
+norm(ra)
+
+b = rand(3,1);
+[rb, hm] = householder(b);
+[b rb eye(b) clean(hm*b) ]
+norm(rb)
+
+[rb2b, hm] = householder(b, 2*b);
+[b rb2b clean(hm*b ./ b) ]  // last column must be uniform
+norm(rb2b)                  // must be 1
+
+c = rand(3,1);
+[rbc, hm] = householder(b,c);
+norm(rbc)          // must be 1
+hm*b ./c           // must be uniform
+
+d = b + %i*c;
+e = rand(3,1) + %i*rand(3,1);
+[rde, hm] = householder(d,e);
+norm(rbc)               // must be 1
+clean(inv(hm) - hm)     // must be zeros(3,3)
+clean(hm' - hm)         // must be zeros(3,3)
+clean(det(hm))          // must be -1
+ ]]></programlisting>
+        
+        <para>APPLICATION : Reflected image of an object</para>
+        <programlisting role="example"><![CDATA[
+// (OA) = [0 0 1] is reflected in O into (OB) = [ 1 1 0.3 ]:
+[n, H] = householder([0 0 1]', [ 1 1 0.3 ]');
+// "n" is the unit vector orthogonal to the reflecting plane
+
+// Emitting object (feature from shell demo):
+u = linspace(0,2*%pi,40);
+v = linspace(0,2*%pi,20);
+Xe = (cos(u).*u)'*(1+cos(v)/2)+10;
+Ye = (u/2)'*sin(v);
+Ze = (sin(u).*u)'*(1+cos(v)/2);
+
+// Reflected object:
+Pe = [ Xe(:)' ; Ye(:)' ; Ze(:)'];
+Pr = H*Pe;
+Xr = matrix(Pr(1,:),40,-1);
+Yr = matrix(Pr(2,:),40,-1);
+Zr = matrix(Pr(3,:),40,-1);
+
+// Reflecting plane containing O: n(1).x + n(2).y + n(3).z = 0
+//   Sampling space:
+x = linspace(min([Xe(:);Xr(:)]), max([Xe(:);Xr(:)]),20);
+y = linspace(min([Ye(:);Yr(:)]), max([Ye(:);Yr(:)]),20);
+[X, Y] = meshgrid(x,y);
+//   Generating the mirror:
+deff("z = mirror(x,y,n)","z = -n(1)/n(3)*x - n(2)/n(3)*y")
+Zm = mirror(X,Y,n);
+
+// Plotting:
+clf
+drawlater
+f = gcf();
+f.color_map = [ 0.8 0.8 0.8 ; jetcolormap(100)];
+surf(Xe,Ye,Ze)
+surf(X,Y,Zm)
+surf(Xr,Yr,Zr)
+a = gca();
+a.isoview = "on";
+a.rotation_angles = [74 123];
+a.children.color_flag = 0;
+a.children.color_mode = 0;
+a.children(1).foreground = color("red");
+a.children(2).foreground = 1;
+a.children(3).foreground = color("green");
+drawnow
+ ]]></programlisting>
+        <scilab:image>
+            householder();
+        </scilab:image>
     </refsection>
     
     <refsection role="see also">
-        
         <title>参照</title>
-        
         <simplelist type="inline">
-            
             <member>
-                
                 <link linkend="qr">qr</link>
-                
             </member>
-            
             <member>
-                
                 <link linkend="givens">givens</link>
-                
             </member>
-            
         </simplelist>
         
     </refsection>
+    <refsection role="history">
+        <title>履歴</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.0</revnumber>
+                <revdescription>
+                    <para>Householder reflexion matrix added as second output parameter.
+                        Demo householder() added. Help page reviewed.
+                    </para>
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
     
 </refentry>
-
index da47d6d..318eef7 100644 (file)
@@ -2,7 +2,8 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
- * 
+ * Copyright (C) 2015 - Samuel GOUGEON
+ *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
  *
  * This file is hereby licensed under the terms of the GNU GPL v2.0,
  * and continues to be available under such terms.
  * For more information, see the COPYING file which you should have received
  * along with this program.
- *
  -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="householder" xml:lang="en">
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="householder" xml:lang="pt">
     <refnamediv>
         <refname>householder</refname>
         <refpurpose>matriz de reflexão ortogonal de Householder</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Seqüência de Chamamento</title>
-        <synopsis>u=householder(v [,w])</synopsis>
+        <synopsis>         householder         // demo</synopsis>
+        <synopsis> u     = householder(v [,w])</synopsis>
+        <synopsis>[u, H] = householder(v [,w])</synopsis>
     </refsynopsisdiv>
-    <refsection>
+    <refsection role="arguments">
         <title>Parâmetros</title>
         <variablelist>
             <varlistentry>
                 <term>w</term>
                 <listitem>
                     <para>vetor coluna de reais ou complexos com o mesmo tamanho que
-                        <literal>v</literal>. Valor padrão é
-                        <literal>eye(v)</literal>
+                        <varname>v</varname>. Valor padrão é <literal>eye(v)</literal> ((Ox) axis)
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>u</term>
                 <listitem>
-                    <para>vetor coluna de reais ou complexos</para>
+                    <para>
+                        unit vector lying in the <literal>(v,w)</literal> plane and orthogonal
+                        to the bisectrix of <literal>(v,w)</literal>.
+                        Column of size(v) of real or complex numbers.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>H</term>
+                <listitem>
+                    <para>
+                        Orthogonal Householder reflexion matrix: <literal>H= eye() - 2*u*u'</literal>.
+                        <varname>H</varname> is such that <literal>inv(H)==H</literal>,
+                        <literal>H'==H</literal>, and <literal>det(H)==-1</literal>.
+                    </para>
+                    <para>
+                        If <varname>v</varname> and <varname>w</varname> are real,
+                        <literal>H*v</literal> is proportional to <varname>w</varname>.
+                    </para>
                 </listitem>
             </varlistentry>
         </variablelist>
     </refsection>
-    <refsection>
+    <refsection role="description">
         <title>Descrição</title>
         <para>
-            Dados dois vetores coluna <literal>v</literal>, <literal>
-                w
-            </literal>
-            de mesmo tamanho, <literal>householder(v,w)</literal> retorna
-            um vetor coluna unitário <literal>u</literal>, tal que<literal>
-                (eye()-2*u*u')*v
-            </literal>
-            éproporcional a <literal>w</literal>.
-            <literal>(eye()-2*u*u')</literal> é a matriz de reflexão ortogonal de
-            Householder .
+            <literal>householder(..)</literal> computes the unit vector <varname>u</varname>
+            lying in the <literal>(v,w)</literal> plane and orthogonal to the bisectrix of
+            <literal>(v,w)</literal>.
         </para>
         <para>
-            O valor padrão de<literal> w</literal> é<literal> eye(v)</literal>.
-            Neste caso, o vetor <literal> (eye()-2*u*u')*v</literal> é o
-            vetor<literal> eye(v)*norm(v)</literal>.
+            If <varname>v</varname> and <varname>w</varname> are proportional:
+            <itemizedlist>
+                <listitem>
+                    <para>
+                        If they are opposite, <literal>u= v/|v|</literal> is returned.
+                    </para>
+                </listitem>
+                <listitem>
+                    If they are real and have the same direction, <varname>u</varname> is set
+                    in the (xOy) plane with a priori <literal>u(1)>0</literal>, and orthogonal to
+                    <varname>v</varname> (<literal>u'*v==0</literal>). However,
+                    <itemizedlist>
+                        <listitem>
+                            If they are along (Ox), <literal>u = (Oy+)</literal> is returned instead.
+                        </listitem>
+                        <listitem>
+                            If <varname>v</varname> and <varname>w</varname> are scalars with same
+                            signs, the orthogonal sub-space is restricted to <literal>{0}</literal>
+                            that can't be normalized:
+                            <varname>u</varname> and <varname>H</varname> are then set to
+                            <literal>%nan</literal>.
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </itemizedlist>
         </para>
+        
+        <para>
+            If the related reflexion matrix <varname>H</varname> is computed, for any point A
+            of column coordinates <literal>a</literal>, <literal>H*a</literal> are the coordinates of
+            the reflected image of A on P.
+        </para>
+        <note>
+            If <varname>v</varname> or/and <varname>w</varname> are in row, they are priorly
+            transposed into columns.
+        </note>
+        <warning>
+            If <varname>v</varname> or/and <varname>w</varname> are <literal>[]</literal>,
+            <literal>[]</literal> is returned for <varname>u</varname> and <varname>H</varname>.
+        </warning>
+    </refsection>
+    <refsection role="examples">
+        <title>Exemplos</title>
+        <programlisting role="example"><![CDATA[
+a = [ rand(1,1) 0  0 ]';
+[ra hm] = householder(a);
+[a ra hm*a ]
+norm(ra)
+
+b = rand(3,1);
+[rb, hm] = householder(b);
+[b rb eye(b) clean(hm*b) ]
+norm(rb)
+
+[rb2b, hm] = householder(b, 2*b);
+[b rb2b clean(hm*b ./ b) ]  // last column must be uniform
+norm(rb2b)                  // must be 1
+
+c = rand(3,1);
+[rbc, hm] = householder(b,c);
+norm(rbc)          // must be 1
+hm*b ./c           // must be uniform
+
+d = b + %i*c;
+e = rand(3,1) + %i*rand(3,1);
+[rde, hm] = householder(d,e);
+norm(rbc)               // must be 1
+clean(inv(hm) - hm)     // must be zeros(3,3)
+clean(hm' - hm)         // must be zeros(3,3)
+clean(det(hm))          // must be -1
+ ]]></programlisting>
+        
+        <para>APPLICATION : Reflected image of an object</para>
+        <programlisting role="example"><![CDATA[
+// (OA) = [0 0 1] is reflected in O into (OB) = [ 1 1 0.3 ]:
+[n, H] = householder([0 0 1]', [ 1 1 0.3 ]');
+// "n" is the unit vector orthogonal to the reflecting plane
+
+// Emitting object (feature from shell demo):
+u = linspace(0,2*%pi,40);
+v = linspace(0,2*%pi,20);
+Xe = (cos(u).*u)'*(1+cos(v)/2)+10;
+Ye = (u/2)'*sin(v);
+Ze = (sin(u).*u)'*(1+cos(v)/2);
+
+// Reflected object:
+Pe = [ Xe(:)' ; Ye(:)' ; Ze(:)'];
+Pr = H*Pe;
+Xr = matrix(Pr(1,:),40,-1);
+Yr = matrix(Pr(2,:),40,-1);
+Zr = matrix(Pr(3,:),40,-1);
+
+// Reflecting plane containing O: n(1).x + n(2).y + n(3).z = 0
+//   Sampling space:
+x = linspace(min([Xe(:);Xr(:)]), max([Xe(:);Xr(:)]),20);
+y = linspace(min([Ye(:);Yr(:)]), max([Ye(:);Yr(:)]),20);
+[X, Y] = meshgrid(x,y);
+//   Generating the mirror:
+deff("z = mirror(x,y,n)","z = -n(1)/n(3)*x - n(2)/n(3)*y")
+Zm = mirror(X,Y,n);
+
+// Plotting:
+clf
+drawlater
+f = gcf();
+f.color_map = [ 0.8 0.8 0.8 ; jetcolormap(100)];
+surf(Xe,Ye,Ze)
+surf(X,Y,Zm)
+surf(Xr,Yr,Zr)
+a = gca();
+a.isoview = "on";
+a.rotation_angles = [74 123];
+a.children.color_flag = 0;
+a.children.color_mode = 0;
+a.children(1).foreground = color("red");
+a.children(2).foreground = 1;
+a.children(3).foreground = color("green");
+drawnow
+ ]]></programlisting>
+        <scilab:image>
+            householder();
+        </scilab:image>
     </refsection>
     <refsection role="see also">
         <title>Ver Também</title>
             </member>
         </simplelist>
     </refsection>
+    <refsection role="history">
+        <title>Histórico</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.0</revnumber>
+                <revdescription>
+                    <para>Householder reflexion matrix added as second output parameter.
+                        Demo householder() added. Help page reviewed.
+                    </para>
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
diff --git a/scilab/modules/linear_algebra/macros/householder.sce b/scilab/modules/linear_algebra/macros/householder.sce
new file mode 100644 (file)
index 0000000..4d9becb
--- /dev/null
@@ -0,0 +1,60 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2016 - Samuel GOUGEON
+//
+// Copyright (C) 2012 - 2016 - Scilab Enterprises
+//
+// This file is hereby licensed under the terms of the GNU GPL v2.0,
+// pursuant to article 5.3.4 of the CeCILL v.2.1.
+// This file was originally licensed under the terms of the CeCILL v2.1,
+// and continues to be available under such terms.
+// For more information, see the COPYING file which you should have received
+// along with this program.
+
+// Example for householder() : Reflection of an object using Householder matrix
+// ----------------------------------------------------------------------------
+// (OA) = [0 0 1] is reflected in O into (OB) = [ 1 1 0.3 ]:
+[n, H] = householder([0 0 1]', [ 1 1 0.3 ]');
+// "n" is the unit vector orthogonal to the reflecting plane
+
+// Emitting object (feature from shell demo):
+u = linspace(0,2*%pi,40);
+v = linspace(0,2*%pi,20);
+Xe = (cos(u).*u)'*(1+cos(v)/2)+10;
+Ye = (u/2)'*sin(v);
+Ze = (sin(u).*u)'*(1+cos(v)/2);
+
+// Reflected object:
+Pe = [ Xe(:)' ; Ye(:)' ; Ze(:)'];
+Pr = H*Pe;
+Xr = matrix(Pr(1,:),40,-1);
+Yr = matrix(Pr(2,:),40,-1);
+Zr = matrix(Pr(3,:),40,-1);
+
+// Reflecting plane containing O: n(1).x + n(2).y + n(3).z = 0
+//   Sampling space:
+x = linspace(min([Xe(:);Xr(:)]), max([Xe(:);Xr(:)]),20);
+y = linspace(min([Ye(:);Yr(:)]), max([Ye(:);Yr(:)]),20);
+[X, Y] = meshgrid(x,y);
+//   Generating the mirror:
+deff("z = mirror(x,y,n)","z = -n(1)/n(3)*x - n(2)/n(3)*y")
+Zm = mirror(X,Y,n);
+
+// Plotting:
+clf
+drawlater
+f = gcf();
+f.color_map = [ 0.8 0.8 0.8 ; jetcolormap(100)];
+surf(Xe,Ye,Ze)
+surf(X,Y,Zm)
+surf(Xr,Yr,Zr)
+a = gca();
+a.isoview = "on";
+a.rotation_angles = [74 123];
+a.children.color_flag = 0;
+a.children.color_mode = 0;
+a.children(1).foreground = color("red");
+a.children(2).foreground = 1;
+a.children(3).foreground = color("green");
+drawnow
+u = []
+hm = []
index 9d0f544..d480506 100644 (file)
@@ -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
 //
 // 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
diff --git a/scilab/modules/linear_algebra/tests/nonreg_tests/bug_13810.dia.ref b/scilab/modules/linear_algebra/tests/nonreg_tests/bug_13810.dia.ref
new file mode 100644 (file)
index 0000000..481ecbb
--- /dev/null
@@ -0,0 +1,23 @@
+// ============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2015 - Samuel GOUGEON
+//
+//  This file is distributed under the same license as the Scilab package.
+// ============================================================================
+// <-- Non-regression test for bug 13810 -->
+//
+// <-- Bugzilla URL -->
+//     http://bugzilla.scilab.org/13810
+//
+// <-- Short Description -->
+// householder(v, k*v) returned [%nan %nan %nan]'
+assert_checkequal(householder([%e ; 0 ; 0]), [0 ; 1 ; 0]);
+v = [ 2 3 1]';
+ref = [ 0.8320503  -0.5547002  0 ]';
+assert_checkalmostequal(householder(v,v), ref, 1e-6);
+assert_checkalmostequal(householder(v,2*v), ref, 1e-6);
+ref = [  0.5345225  0.8017837  0.2672612 ]';
+assert_checkalmostequal(householder(v,-2*v), ref, 1e-6);
+v = [ 2 -%i 1+%i].';
+ref = [0.7559289 -0.3779645*%i  0.3779645 + 0.3779645*%i].';
+assert_checkalmostequal(householder(v,v), ref, 1e-6);
diff --git a/scilab/modules/linear_algebra/tests/nonreg_tests/bug_13810.tst b/scilab/modules/linear_algebra/tests/nonreg_tests/bug_13810.tst
new file mode 100644 (file)
index 0000000..481ecbb
--- /dev/null
@@ -0,0 +1,23 @@
+// ============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2015 - Samuel GOUGEON
+//
+//  This file is distributed under the same license as the Scilab package.
+// ============================================================================
+// <-- Non-regression test for bug 13810 -->
+//
+// <-- Bugzilla URL -->
+//     http://bugzilla.scilab.org/13810
+//
+// <-- Short Description -->
+// householder(v, k*v) returned [%nan %nan %nan]'
+assert_checkequal(householder([%e ; 0 ; 0]), [0 ; 1 ; 0]);
+v = [ 2 3 1]';
+ref = [ 0.8320503  -0.5547002  0 ]';
+assert_checkalmostequal(householder(v,v), ref, 1e-6);
+assert_checkalmostequal(householder(v,2*v), ref, 1e-6);
+ref = [  0.5345225  0.8017837  0.2672612 ]';
+assert_checkalmostequal(householder(v,-2*v), ref, 1e-6);
+v = [ 2 -%i 1+%i].';
+ref = [0.7559289 -0.3779645*%i  0.3779645 + 0.3779645*%i].';
+assert_checkalmostequal(householder(v,v), ref, 1e-6);