1 /*
2 Copyright (C) 1992-2008 The University of Tennessee
3 All rights reserved.
4
5 Copyright (c) 2009-2014, Jack Poulson
6 All rights reserved.
7
8 This file is loosely based upon the LAPACK routines dlarfg.f and zlarfg.f.
9
10 This file is part of Elemental and is under the BSD 2-Clause License,
11 which can be found in the LICENSE file in the root directory, or at
12 http://opensource.org/licenses/BSD-2-Clause
13 */
14 #pragma once
15 #ifndef ELEM_REFLECTOR_HPP
16 #define ELEM_REFLECTOR_HPP
17
18 #include ELEM_NRM2_INC
19 #include ELEM_SCALE_INC
20
21 #include "./Reflector/Col.hpp"
22 #include "./Reflector/Row.hpp"
23
24 namespace elem {
25
26 //
27 // The LAPACK convention defines tau such that
28 //
29 // H = I - tau [1; v] [1, v'],
30 //
31 // but adjoint(H) [chi; x] = [beta; 0].
32 //
33 // Elemental simply uses H [chi; x] = [beta; 0].
34 //
35 // On exit, chi is overwritten with beta, and x is overwritten with v.
36 //
37 // Another major difference from LAPACK is in the treatment of the special case
38 // of x=0, where LAPACK would put H := I, which is not a valid Householder
39 // reflector. We instead use the valid Householder reflector:
40 // H [chi; 0] = [-chi; 0],
41 // which is accomplished by setting tau=2, and v=0.
42 //
43
44 // TODO: Switch to 1/tau to be simplify discussions of UT transforms
45
46 template<typename F>
47 inline F
LeftReflector(Matrix<F> & chi,Matrix<F> & x)48 LeftReflector( Matrix<F>& chi, Matrix<F>& x )
49 {
50 DEBUG_ONLY(
51 CallStackEntry cse("LeftReflector");
52 if( chi.Height() != 1 || chi.Width() != 1 )
53 LogicError("chi must be a scalar");
54 if( x.Height() != 1 && x.Width() != 1 )
55 LogicError("x must be a vector");
56 )
57 typedef Base<F> Real;
58
59 Real norm = Nrm2( x );
60 F alpha = chi.Get(0,0);
61
62 if( norm == Real(0) && ImagPart(alpha) == Real(0) )
63 {
64 chi.Set(0,0,-chi.Get(0,0));
65 return F(2);
66 }
67
68 Real beta;
69 if( RealPart(alpha) <= 0 )
70 beta = lapack::SafeNorm( alpha, norm );
71 else
72 beta = -lapack::SafeNorm( alpha, norm );
73
74 // Rescale if the vector is too small
75 const Real safeMin = lapack::MachineSafeMin<Real>();
76 const Real epsilon = lapack::MachineEpsilon<Real>();
77 const Real safeInv = safeMin/epsilon;
78 Int count = 0;
79 if( Abs(beta) < safeInv )
80 {
81 Real invOfSafeInv = Real(1)/safeInv;
82 do
83 {
84 ++count;
85 Scale( invOfSafeInv, x );
86 alpha *= invOfSafeInv;
87 beta *= invOfSafeInv;
88 } while( Abs(beta) < safeInv );
89
90 norm = Nrm2( x );
91 if( RealPart(alpha) <= 0 )
92 beta = lapack::SafeNorm( alpha, norm );
93 else
94 beta = -lapack::SafeNorm( alpha, norm );
95 }
96
97 F tau = (beta-Conj(alpha)) / beta;
98 Scale( Real(1)/(alpha-beta), x );
99
100 // Undo the scaling
101 for( Int j=0; j<count; ++j )
102 beta *= safeInv;
103
104 chi.Set(0,0,beta);
105 return tau;
106 }
107
108 template<typename F>
109 inline F
LeftReflector(F & chi,Matrix<F> & x)110 LeftReflector( F& chi, Matrix<F>& x )
111 {
112 DEBUG_ONLY(
113 CallStackEntry cse("LeftReflector");
114 if( x.Height() != 1 && x.Width() != 1 )
115 LogicError("x must be a vector");
116 )
117 typedef Base<F> Real;
118
119 Real norm = Nrm2( x );
120 F alpha = chi;
121
122 if( norm == Real(0) && ImagPart(alpha) == Real(0) )
123 {
124 chi = -chi;
125 return F(2);
126 }
127
128 Real beta;
129 if( RealPart(alpha) <= 0 )
130 beta = lapack::SafeNorm( alpha, norm );
131 else
132 beta = -lapack::SafeNorm( alpha, norm );
133
134 // Rescale if the vector is too small
135 const Real safeMin = lapack::MachineSafeMin<Real>();
136 const Real epsilon = lapack::MachineEpsilon<Real>();
137 const Real safeInv = safeMin/epsilon;
138 Int count = 0;
139 if( Abs(beta) < safeInv )
140 {
141 Real invOfSafeInv = Real(1)/safeInv;
142 do
143 {
144 ++count;
145 Scale( invOfSafeInv, x );
146 alpha *= invOfSafeInv;
147 beta *= invOfSafeInv;
148 } while( Abs(beta) < safeInv );
149
150 norm = Nrm2( x );
151 if( RealPart(alpha) <= 0 )
152 beta = lapack::SafeNorm( alpha, norm );
153 else
154 beta = -lapack::SafeNorm( alpha, norm );
155 }
156
157 F tau = (beta-Conj(alpha)) / beta;
158 Scale( Real(1)/(alpha-beta), x );
159
160 // Undo the scaling
161 for( Int j=0; j<count; ++j )
162 beta *= safeInv;
163
164 chi = beta;
165 return tau;
166 }
167
168
169 template<typename F,Dist U,Dist V>
170 inline F
LeftReflector(DistMatrix<F,U,V> & chi,DistMatrix<F,U,V> & x)171 LeftReflector( DistMatrix<F,U,V>& chi, DistMatrix<F,U,V>& x )
172 {
173 DEBUG_ONLY(
174 CallStackEntry cse("LeftReflector");
175 if( chi.Grid() != x.Grid() )
176 LogicError("chi and x must be distributed over the same grid");
177 if( chi.Height() != 1 || chi.Width() != 1 )
178 LogicError("chi must be a scalar");
179 if( x.Width() != 1 )
180 LogicError("x must be a column vector");
181 )
182 F tau;
183 if( x.RowRank() == x.RowAlign() )
184 tau = reflector::Col( chi, x );
185 mpi::Broadcast( tau, x.RowAlign(), x.RowComm() );
186 return tau;
187 }
188
189 template<typename F,Dist U,Dist V>
190 inline F
LeftReflector(F & chi,DistMatrix<F,U,V> & x)191 LeftReflector( F& chi, DistMatrix<F,U,V>& x )
192 {
193 DEBUG_ONLY(
194 CallStackEntry cse("LeftReflector");
195 if( x.Width() != 1 )
196 LogicError("x must be a column vector");
197 )
198 F tau;
199 if( x.RowRank() == x.RowAlign() )
200 tau = reflector::Col( chi, x );
201 mpi::Broadcast( tau, x.RowAlign(), x.RowComm() );
202 return tau;
203 }
204
205 //
206 // Defines tau and v such that
207 //
208 // H = I - tau [1; v] [1, v'],
209 //
210 // and [chi x] H = [beta 0]
211 //
212
213 template<typename F>
214 inline F
RightReflector(Matrix<F> & chi,Matrix<F> & x)215 RightReflector( Matrix<F>& chi, Matrix<F>& x )
216 {
217 DEBUG_ONLY(
218 CallStackEntry cse("RightReflector");
219 if( chi.Height() != 1 || chi.Width() != 1 )
220 LogicError("chi must be a scalar");
221 if( x.Height() != 1 && x.Width() != 1 )
222 LogicError("x must be a vector");
223 )
224 const F tau = LeftReflector( chi, x );
225 // There is no need to conjugate chi, it should be real now
226 Conjugate( x );
227 return tau;
228 }
229
230 template<typename F>
231 inline F
RightReflector(F & chi,Matrix<F> & x)232 RightReflector( F& chi, Matrix<F>& x )
233 {
234 DEBUG_ONLY(
235 CallStackEntry cse("RightReflector");
236 if( x.Height() != 1 && x.Width() != 1 )
237 LogicError("x must be a vector");
238 )
239 const F tau = LeftReflector( chi, x );
240 // There is no need to conjugate chi, it should be real now
241 Conjugate( x );
242 return tau;
243 }
244
245 template<typename F,Dist U,Dist V>
246 inline F
RightReflector(DistMatrix<F,U,V> & chi,DistMatrix<F,U,V> & x)247 RightReflector( DistMatrix<F,U,V>& chi, DistMatrix<F,U,V>& x )
248 {
249 DEBUG_ONLY(
250 CallStackEntry cse("RightReflector");
251 if( chi.Grid() != x.Grid() )
252 LogicError("chi and x must be distributed over the same grid");
253 if( chi.Height() != 1 || chi.Width() != 1 )
254 LogicError("chi must be a scalar");
255 if( x.Height() != 1 )
256 LogicError("x must be a row vector");
257 )
258 F tau;
259 if( x.ColRank() == x.ColAlign() )
260 tau = reflector::Row( chi, x );
261 mpi::Broadcast( tau, x.ColAlign(), x.ColComm() );
262 return tau;
263 }
264
265 template<typename F,Dist U,Dist V>
266 inline F
RightReflector(F & chi,DistMatrix<F,U,V> & x)267 RightReflector( F& chi, DistMatrix<F,U,V>& x )
268 {
269 DEBUG_ONLY(
270 CallStackEntry cse("RightReflector");
271 if( x.Height() != 1 )
272 LogicError("x must be a row vector");
273 )
274 F tau;
275 if( x.ColRank() == x.ColAlign() )
276 tau = reflector::Row( chi, x );
277 mpi::Broadcast( tau, x.ColAlign(), x.ColComm() );
278 return tau;
279 }
280
281 } // namespace elem
282
283 #endif // ifndef ELEM_REFLECTOR_HPP
284