1 /*
2 * Copyright (c) 2011, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 /* Based on
34 *
35 SUBROUTINE DGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO )
36 *
37 * -- LAPACK routine (version 3.3.1) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * -- April 2011 --
41 */
42
43 #ifndef FLENS_LAPACK_EIG_HD2_TCC
44 #define FLENS_LAPACK_EIG_HD2_TCC 1
45
46 #include <flens/blas/blas.h>
47 #include <flens/lapack/lapack.h>
48
49 namespace flens { namespace lapack {
50
51 //== generic lapack implementation =============================================
52
53 template <typename IndexType, typename MA, typename VTAU, typename VW>
54 void
55 hd2_generic(IndexType iLo,
56 IndexType iHi,
57 GeMatrix<MA> &A,
58 DenseVector<VTAU> &tau,
59 DenseVector<VW> &work)
60 {
61 using lapack::larf;
62 using lapack::larfg;
63 using std::max;
64 using std::min;
65
66 typedef typename GeMatrix<MA>::ElementType T;
67
68 const Underscore<IndexType> _;
69
70 const IndexType n = A.numCols();
71
72 for (IndexType i=iLo; i<iHi; ++i) {
73 //
74 // Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
75 //
76 larfg(iHi-i, A(i+1,i), A(_(min(i+2,n),iHi),i), tau(i));
77
78 const T Aii = A(i+1,i);
79 A(i+1,i) = 1;
80 //
81 // Apply H(i) to A(1:ihi,i+1:ihi) from the right
82 //
83 larf(Right, A(_(i+1,iHi),i), tau(i), A(_(1,iHi),_(i+1,iHi)), work);
84 //
85 // Apply H(i) to A(i+1:ihi,i+1:n) from the left
86 //
87 larf(Left, A(_(i+1,iHi),i), tau(i), A(_(i+1,iHi),_(i+1,n)), work);
88
89 A(i+1,i) = Aii;
90 }
91 }
92
93 //== interface for native lapack ===============================================
94
95 #ifdef CHECK_CXXLAPACK
96
97 template <typename IndexType, typename MA, typename VTAU, typename VWORK>
98 void
99 hd2_native(IndexType iLo, IndexType iHi, GeMatrix<MA> &A,
100 DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
101 {
102 // TODO: add this assertion to other ..._native implementations
103 assert(A.order()==ColMajor);
104 typedef typename GeMatrix<MA>::ElementType T;
105
106 const INTEGER N = A.numCols();
107 const INTEGER ILO = iLo;
108 const INTEGER IHI = iHi;
109 const INTEGER LDA = A.leadingDimension();
110
111 INTEGER INFO;
112
113 if (IsSame<T,DOUBLE>::value) {
114 LAPACK_IMPL(dgehd2)(&N,
115 &ILO,
116 &IHI,
117 A.data(),
118 &LDA,
119 tau.data(),
120 work.data(),
121 &INFO);
122 } else {
123 ASSERT(0);
124 }
125
126 ASSERT(INFO==0);
127 }
128
129 #endif // CHECK_CXXLAPACK
130
131 //== public interface ==========================================================
132
133 template <typename IndexType, typename MA, typename VTAU, typename VW>
134 void
135 hd2(IndexType iLo,
136 IndexType iHi,
137 GeMatrix<MA> &A,
138 DenseVector<VTAU> &tau,
139 DenseVector<VW> &work)
140 {
141 LAPACK_DEBUG_OUT("hd2");
142
143 using std::max;
144 //
145 // Test the input parameters
146 //
147 ASSERT(A.firstRow()==1);
148 ASSERT(A.firstCol()==1);
149 ASSERT(A.numRows()==A.numCols());
150 ASSERT(tau.firstIndex()<=iLo);
151 ASSERT(tau.lastIndex()>=iHi-1);
152 ASSERT(tau.inc()>0);
153 ASSERT(work.length()>=A.numRows());
154
155 ASSERT(1<=iLo);
156 ASSERT(iLo<=iHi);
157 ASSERT(iHi<=max(IndexType(1), A.numCols()));
158
159 # ifdef CHECK_CXXLAPACK
160 //
161 // Make copies of output arguments
162 //
163 typename GeMatrix<MA>::NoView _A = A;
164 typename DenseVector<VTAU>::NoView _tau = tau;
165 typename DenseVector<VW>::NoView _work = work;
166 # endif
167
168 //
169 // Call implementation
170 //
171 hd2_generic(iLo, iHi, A, tau, work);
172
173 # ifdef CHECK_CXXLAPACK
174 //
175 // Compare results
176 //
177 hd2_native(iLo, iHi, _A, _tau, _work);
178
179 if (! isIdentical(A, _A, "A", "A_")) {
180 std::cerr << "CXXLAPACK: A = " << A << std::endl;
181 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
182 ASSERT(0);
183 }
184
185 if (! isIdentical(tau, _tau, "tau", "tau_")) {
186 std::cerr << "CXXLAPACK: tau = " << tau << std::endl;
187 std::cerr << "F77LAPACK: _tau = " << _tau << std::endl;
188 ASSERT(0);
189 }
190
191 if (! isIdentical(work, _work, "work", "work_")) {
192 std::cerr << "CXXLAPACK: work = " << work << std::endl;
193 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
194 ASSERT(0);
195 }
196 # endif
197 }
198
199 //-- forwarding ----------------------------------------------------------------
200 template <typename IndexType, typename MA, typename VTAU, typename VW>
201 void
202 hd2(IndexType iLo,
203 IndexType iHi,
204 MA &&A,
205 VTAU &&tau,
206 VW &&work)
207 {
208 hd2(iLo, iHi, A, tau, work);
209 }
210
211 } } // namespace lapack, flens
212
213 #endif // FLENS_LAPACK_EIG_HD2_TCC
2 * Copyright (c) 2011, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 /* Based on
34 *
35 SUBROUTINE DGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO )
36 *
37 * -- LAPACK routine (version 3.3.1) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * -- April 2011 --
41 */
42
43 #ifndef FLENS_LAPACK_EIG_HD2_TCC
44 #define FLENS_LAPACK_EIG_HD2_TCC 1
45
46 #include <flens/blas/blas.h>
47 #include <flens/lapack/lapack.h>
48
49 namespace flens { namespace lapack {
50
51 //== generic lapack implementation =============================================
52
53 template <typename IndexType, typename MA, typename VTAU, typename VW>
54 void
55 hd2_generic(IndexType iLo,
56 IndexType iHi,
57 GeMatrix<MA> &A,
58 DenseVector<VTAU> &tau,
59 DenseVector<VW> &work)
60 {
61 using lapack::larf;
62 using lapack::larfg;
63 using std::max;
64 using std::min;
65
66 typedef typename GeMatrix<MA>::ElementType T;
67
68 const Underscore<IndexType> _;
69
70 const IndexType n = A.numCols();
71
72 for (IndexType i=iLo; i<iHi; ++i) {
73 //
74 // Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
75 //
76 larfg(iHi-i, A(i+1,i), A(_(min(i+2,n),iHi),i), tau(i));
77
78 const T Aii = A(i+1,i);
79 A(i+1,i) = 1;
80 //
81 // Apply H(i) to A(1:ihi,i+1:ihi) from the right
82 //
83 larf(Right, A(_(i+1,iHi),i), tau(i), A(_(1,iHi),_(i+1,iHi)), work);
84 //
85 // Apply H(i) to A(i+1:ihi,i+1:n) from the left
86 //
87 larf(Left, A(_(i+1,iHi),i), tau(i), A(_(i+1,iHi),_(i+1,n)), work);
88
89 A(i+1,i) = Aii;
90 }
91 }
92
93 //== interface for native lapack ===============================================
94
95 #ifdef CHECK_CXXLAPACK
96
97 template <typename IndexType, typename MA, typename VTAU, typename VWORK>
98 void
99 hd2_native(IndexType iLo, IndexType iHi, GeMatrix<MA> &A,
100 DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
101 {
102 // TODO: add this assertion to other ..._native implementations
103 assert(A.order()==ColMajor);
104 typedef typename GeMatrix<MA>::ElementType T;
105
106 const INTEGER N = A.numCols();
107 const INTEGER ILO = iLo;
108 const INTEGER IHI = iHi;
109 const INTEGER LDA = A.leadingDimension();
110
111 INTEGER INFO;
112
113 if (IsSame<T,DOUBLE>::value) {
114 LAPACK_IMPL(dgehd2)(&N,
115 &ILO,
116 &IHI,
117 A.data(),
118 &LDA,
119 tau.data(),
120 work.data(),
121 &INFO);
122 } else {
123 ASSERT(0);
124 }
125
126 ASSERT(INFO==0);
127 }
128
129 #endif // CHECK_CXXLAPACK
130
131 //== public interface ==========================================================
132
133 template <typename IndexType, typename MA, typename VTAU, typename VW>
134 void
135 hd2(IndexType iLo,
136 IndexType iHi,
137 GeMatrix<MA> &A,
138 DenseVector<VTAU> &tau,
139 DenseVector<VW> &work)
140 {
141 LAPACK_DEBUG_OUT("hd2");
142
143 using std::max;
144 //
145 // Test the input parameters
146 //
147 ASSERT(A.firstRow()==1);
148 ASSERT(A.firstCol()==1);
149 ASSERT(A.numRows()==A.numCols());
150 ASSERT(tau.firstIndex()<=iLo);
151 ASSERT(tau.lastIndex()>=iHi-1);
152 ASSERT(tau.inc()>0);
153 ASSERT(work.length()>=A.numRows());
154
155 ASSERT(1<=iLo);
156 ASSERT(iLo<=iHi);
157 ASSERT(iHi<=max(IndexType(1), A.numCols()));
158
159 # ifdef CHECK_CXXLAPACK
160 //
161 // Make copies of output arguments
162 //
163 typename GeMatrix<MA>::NoView _A = A;
164 typename DenseVector<VTAU>::NoView _tau = tau;
165 typename DenseVector<VW>::NoView _work = work;
166 # endif
167
168 //
169 // Call implementation
170 //
171 hd2_generic(iLo, iHi, A, tau, work);
172
173 # ifdef CHECK_CXXLAPACK
174 //
175 // Compare results
176 //
177 hd2_native(iLo, iHi, _A, _tau, _work);
178
179 if (! isIdentical(A, _A, "A", "A_")) {
180 std::cerr << "CXXLAPACK: A = " << A << std::endl;
181 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
182 ASSERT(0);
183 }
184
185 if (! isIdentical(tau, _tau, "tau", "tau_")) {
186 std::cerr << "CXXLAPACK: tau = " << tau << std::endl;
187 std::cerr << "F77LAPACK: _tau = " << _tau << std::endl;
188 ASSERT(0);
189 }
190
191 if (! isIdentical(work, _work, "work", "work_")) {
192 std::cerr << "CXXLAPACK: work = " << work << std::endl;
193 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
194 ASSERT(0);
195 }
196 # endif
197 }
198
199 //-- forwarding ----------------------------------------------------------------
200 template <typename IndexType, typename MA, typename VTAU, typename VW>
201 void
202 hd2(IndexType iLo,
203 IndexType iHi,
204 MA &&A,
205 VTAU &&tau,
206 VW &&work)
207 {
208 hd2(iLo, iHi, A, tau, work);
209 }
210
211 } } // namespace lapack, flens
212
213 #endif // FLENS_LAPACK_EIG_HD2_TCC