1 /*
2 Copyright (c) 2009-2014, Jack Poulson
3 All rights reserved.
4
5 This file is part of Elemental and is under the BSD 2-Clause License,
6 which can be found in the LICENSE file in the root directory, or at
7 http://opensource.org/licenses/BSD-2-Clause
8 */
9 #pragma once
10 #ifndef ELEM_LQ_HOUSEHOLDER_HPP
11 #define ELEM_LQ_HOUSEHOLDER_HPP
12
13 #include "./ApplyQ.hpp"
14 #include "./PanelHouseholder.hpp"
15
16 namespace elem {
17 namespace lq {
18
19 template<typename F>
20 inline void
Householder(Matrix<F> & A,Matrix<F> & t,Matrix<Base<F>> & d)21 Householder( Matrix<F>& A, Matrix<F>& t, Matrix<Base<F>>& d )
22 {
23 DEBUG_ONLY(CallStackEntry cse("lq::Householder"))
24 const Int m = A.Height();
25 const Int n = A.Width();
26 const Int minDim = Min(m,n);
27 t.Resize( minDim, 1 );
28 d.Resize( minDim, 1 );
29
30 const Int bsize = Blocksize();
31 for( Int k=0; k<minDim; k+=bsize )
32 {
33 const Int nb = Min(bsize,minDim-k);
34 auto ATopPan = ViewRange( A, k, k, k+nb, n );
35 auto ABottomPan = ViewRange( A, k+nb, k, m, n );
36 auto t1 = View( t, k, 0, nb, 1 );
37 auto d1 = View( d, k, 0, nb, 1 );
38
39 PanelHouseholder( ATopPan, t1, d1 );
40 ApplyQ( RIGHT, ADJOINT, ATopPan, t1, d1, ABottomPan );
41 }
42 }
43
44 template<typename F>
45 inline void
Householder(Matrix<F> & A)46 Householder( Matrix<F>& A )
47 {
48 DEBUG_ONLY(CallStackEntry cse("lq::Householder"))
49 Matrix<F> t;
50 Matrix<Base<F>> d;
51 Householder( A, t, d );
52 MakeTriangular( LOWER, A );
53 }
54
55 template<typename F>
56 inline void
Householder(DistMatrix<F> & A,DistMatrix<F,MD,STAR> & t,DistMatrix<Base<F>,MD,STAR> & d)57 Householder
58 ( DistMatrix<F>& A, DistMatrix<F,MD,STAR>& t, DistMatrix<Base<F>,MD,STAR>& d )
59 {
60 DEBUG_ONLY(
61 CallStackEntry cse("Householder");
62 if( A.Grid() != t.Grid() || t.Grid() != d.Grid() )
63 LogicError("{A,t,d} must be distributed over the same grid");
64 )
65 const Int m = A.Height();
66 const Int n = A.Width();
67 const Int minDim = Min(m,n);
68
69 t.SetRoot( A.DiagonalRoot() );
70 d.SetRoot( A.DiagonalRoot() );
71 t.AlignCols( A.DiagonalAlign() );
72 d.AlignCols( A.DiagonalAlign() );
73 t.Resize( minDim, 1 );
74 d.Resize( minDim, 1 );
75
76 const Int bsize = Blocksize();
77 for( Int k=0; k<minDim; k+=bsize )
78 {
79 const Int nb = Min(bsize,minDim-k);
80 auto ATopPan = ViewRange( A, k, k, k+nb, n );
81 auto ABottomPan = ViewRange( A, k+nb, k, m, n );
82 auto t1 = View( t, k, 0, nb, 1 );
83 auto d1 = View( d, k, 0, nb, 1 );
84
85 PanelHouseholder( ATopPan, t1, d1 );
86 ApplyQ( RIGHT, ADJOINT, ATopPan, t1, d1, ABottomPan );
87 }
88 }
89
90 template<typename F>
91 inline void
Householder(DistMatrix<F> & A)92 Householder( DistMatrix<F>& A )
93 {
94 DEBUG_ONLY(CallStackEntry cse("Householder"))
95 DistMatrix<F,MD,STAR> t(A.Grid());
96 DistMatrix<Base<F>,MD,STAR> d(A.Grid());
97 Householder( A, t, d );
98 MakeTriangular( LOWER, A );
99 }
100
101 } // namespace lq
102 } // namespace elem
103
104 #endif // ifndef ELEM_LQ_HOUSEHOLDER_HPP
105