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 DORGHR( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
36 *
37 * -- LAPACK routine (version 3.2) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * November 2006
41 */
42
43 #ifndef FLENS_LAPACK_EIG_ORGHR_TCC
44 #define FLENS_LAPACK_EIG_ORGHR_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>
54 IndexType
55 orghr_generic_wsq(IndexType iLo,
56 IndexType iHi,
57 const GeMatrix<MA> &A,
58 const DenseVector<VTAU> &)
59 {
60 using std::max;
61
62 typedef typename GeMatrix<MA>::ElementType T;
63
64 const Underscore<IndexType> _;
65 const IndexType n = A.numRows();
66 const IndexType nh = iHi - iLo;
67 const IndexType nb = ilaenv<T>(1, "ORGQR", "", nh, nh, nh);
68
69 return max(IndexType(1), nh)*nb;
70 }
71
72 template <typename IndexType, typename MA, typename VTAU, typename VW>
73 void
74 orghr_generic(IndexType iLo,
75 IndexType iHi,
76 GeMatrix<MA> &A,
77 const DenseVector<VTAU> &tau,
78 DenseVector<VW> &work)
79 {
80 using std::max;
81
82 typedef typename GeMatrix<MA>::ElementType T;
83
84 const Underscore<IndexType> _;
85 const IndexType n = A.numRows();
86 const IndexType nh = iHi - iLo;
87 const IndexType nb = ilaenv<T>(1, "ORGQR", "", nh, nh, nh);
88 const IndexType lWorkOpt = max(IndexType(1), nh)*nb;
89
90 //
91 // Apply workspace query if needed
92 //
93 if (work.length()==0) {
94 work.resize(lWorkOpt);
95 }
96 //
97 // Quick return if possible
98 //
99 if (n==0) {
100 work(1) = 1;
101 }
102 //
103 // Shift the vectors which define the elementary reflectors one
104 // column to the right, and set the first ilo and the last n-ihi
105 // rows and columns to those of the unit matrix
106 //
107 for (IndexType j=iHi; j>iLo; --j) {
108 for (IndexType i=1; i<=j-1; ++i) {
109 A(i,j) = T(0);
110 }
111 for (IndexType i=j+1; i<=iHi; ++i) {
112 A(i,j) = A(i,j-1);
113 }
114 for (IndexType i=iHi+1; i<=n; ++i) {
115 A(i,j) = T(0);
116 }
117 }
118 for (IndexType j=1; j<=iLo; ++j) {
119 for (IndexType i=1; i<=n; ++i) {
120 A(i,j) = T(0);
121 }
122 A(j,j) = 1;
123 }
124 for (IndexType j=iHi+1; j<=n; ++j) {
125 for (IndexType i=1; i<=n; ++i) {
126 A(i,j) = T(0);
127 }
128 A(j,j) = T(1);
129 }
130
131 if (nh>0) {
132 //
133 // Generate Q(ilo+1:ihi,ilo+1:ihi)
134 //
135 orgqr(nh, A(_(iLo+1,iHi),_(iLo+1,iHi)), tau(_(iLo,iHi-1)), work);
136 }
137 work(1) = lWorkOpt;
138 }
139
140 //== interface for native lapack ===============================================
141
142 #ifdef CHECK_CXXLAPACK
143
144 template <typename IndexType, typename MA, typename VTAU>
145 IndexType
146 orghr_native_wsq(IndexType iLo,
147 IndexType iHi,
148 const GeMatrix<MA> &A,
149 const DenseVector<VTAU> &tau)
150 {
151 typedef typename GeMatrix<MA>::ElementType T;
152
153 const INTEGER N = A.numRows();
154 const INTEGER ILO = iLo;
155 const INTEGER IHI = iHi;
156 const INTEGER LDA = A.leadingDimension();
157 T WORK;
158 T DUMMY;
159 const INTEGER LWORK = -1;
160 INTEGER INFO;
161
162 if (IsSame<T, DOUBLE>::value) {
163 LAPACK_IMPL(dorghr)(&N,
164 &ILO,
165 &IHI,
166 &DUMMY,
167 &LDA,
168 tau.data(),
169 &WORK,
170 &LWORK,
171 &INFO);
172 } else {
173 ASSERT(0);
174 }
175 ASSERT(INFO==0);
176 return WORK;
177 }
178
179
180 template <typename IndexType, typename MA, typename VTAU, typename VW>
181 void
182 orghr_native(IndexType iLo,
183 IndexType iHi,
184 GeMatrix<MA> &A,
185 const DenseVector<VTAU> &tau,
186 DenseVector<VW> &work)
187 {
188 typedef typename GeMatrix<MA>::ElementType T;
189
190 const INTEGER N = A.numRows();
191 const INTEGER ILO = iLo;
192 const INTEGER IHI = iHi;
193 const INTEGER LDA = A.leadingDimension();
194 const INTEGER LWORK = work.length();
195 INTEGER INFO;
196
197 if (IsSame<T, DOUBLE>::value) {
198 LAPACK_IMPL(dorghr)(&N,
199 &ILO,
200 &IHI,
201 A.data(),
202 &LDA,
203 tau.data(),
204 work.data(),
205 &LWORK,
206 &INFO);
207 } else {
208 ASSERT(0);
209 }
210 ASSERT(INFO==0);
211 }
212
213 #endif // CHECK_CXXLAPACK
214
215 //== public interface ==========================================================
216
217 template <typename IndexType, typename MA, typename VTAU>
218 IndexType
219 orghr_wsq(IndexType iLo,
220 IndexType iHi,
221 const GeMatrix<MA> &A,
222 const DenseVector<VTAU> &tau)
223 {
224 //
225 // Test the input parameters
226 //
227 # ifndef NDEBUG
228 ASSERT(A.numRows()==A.numCols());
229
230 const IndexType n = A.numCols();
231 if (n==0) {
232 ASSERT(iLo==1);
233 ASSERT(iHi==0);
234 } else {
235 ASSERT(1<=iLo);
236 ASSERT(iLo<=iHi);
237 ASSERT(iHi<=n);
238 }
239 ASSERT(tau.length()==(n-1));
240 # endif
241
242 //
243 // Call implementation
244 //
245 IndexType ws = orghr_generic_wsq(iLo, iHi, A, tau);
246
247 # ifdef CHECK_CXXLAPACK
248 //
249 // Compare results
250 //
251 IndexType _ws = orghr_native_wsq(iLo, iHi, A, tau);
252
253 if (ws!=_ws) {
254 std::cerr << "CXXLAPACK: ws = " << ws << std::endl;
255 std::cerr << "F77LAPACK: _ws = " << _ws << std::endl;
256 ASSERT(0);
257 }
258 # endif
259
260 return ws;
261 }
262
263 template <typename IndexType, typename MA, typename VTAU, typename VW>
264 void
265 orghr(IndexType iLo,
266 IndexType iHi,
267 GeMatrix<MA> &A,
268 const DenseVector<VTAU> &tau,
269 DenseVector<VW> &work)
270 {
271 LAPACK_DEBUG_OUT("orghr");
272
273 //
274 // Test the input parameters
275 //
276 # ifndef NDEBUG
277 ASSERT(A.numRows()==A.numCols());
278
279 const IndexType n = A.numCols();
280 if (n==0) {
281 ASSERT(iLo==1);
282 ASSERT(iHi==0);
283 } else {
284 ASSERT(1<=iLo);
285 ASSERT(iLo<=iHi);
286 ASSERT(iHi<=n);
287 }
288 ASSERT(tau.length()==(n-1));
289 ASSERT(work.length()>=iHi-iLo);
290 # endif
291
292 //
293 // Make copies of output arguments
294 //
295 # ifdef CHECK_CXXLAPACK
296 typename GeMatrix<MA>::NoView _A = A;
297 typename DenseVector<VW>::NoView _work = work;
298 # endif
299
300 //
301 // Call implementation
302 //
303 orghr_generic(iLo, iHi, A, tau, work);
304
305 # ifdef CHECK_CXXLAPACK
306 //
307 // Compare results
308 //
309 if (_work.length()==0) {
310 _work.resize(work.length());
311 }
312 orghr_native(iLo, iHi, _A, tau, _work);
313
314 bool failed = false;
315 if (! isIdentical(A, _A, " A", "A_")) {
316 std::cerr << "CXXLAPACK: A = " << A << std::endl;
317 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
318 failed = true;
319 }
320
321 if (! isIdentical(work, _work, " work", "_work")) {
322 std::cerr << "CXXLAPACK: work = " << work << std::endl;
323 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
324 failed = true;
325 }
326
327 if (failed) {
328 std::cerr << "error in: hrd.tcc" << std::endl;
329 ASSERT(0);
330 } else {
331 // std::cerr << "passed: hrd.tcc" << std::endl;
332 }
333 # endif
334 }
335
336 //-- forwarding ----------------------------------------------------------------
337 template <typename IndexType, typename MA, typename VTAU>
338 IndexType
339 orghr_wsq(IndexType iLo,
340 IndexType iHi,
341 const MA &&A,
342 const VTAU &tau)
343 {
344 return orghr_wsq(iLo, iHi, A, tau);
345 }
346
347 template <typename IndexType, typename MA, typename VTAU, typename VW>
348 void
349 orghr(IndexType iLo,
350 IndexType iHi,
351 MA &&A,
352 const VTAU &tau,
353 VW &&work)
354 {
355 orghr(iLo, iHi, A, tau, work);
356 }
357
358 } } // namespace lapack, flens
359
360 #endif // FLENS_LAPACK_EIG_ORGHR_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 DORGHR( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
36 *
37 * -- LAPACK routine (version 3.2) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * November 2006
41 */
42
43 #ifndef FLENS_LAPACK_EIG_ORGHR_TCC
44 #define FLENS_LAPACK_EIG_ORGHR_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>
54 IndexType
55 orghr_generic_wsq(IndexType iLo,
56 IndexType iHi,
57 const GeMatrix<MA> &A,
58 const DenseVector<VTAU> &)
59 {
60 using std::max;
61
62 typedef typename GeMatrix<MA>::ElementType T;
63
64 const Underscore<IndexType> _;
65 const IndexType n = A.numRows();
66 const IndexType nh = iHi - iLo;
67 const IndexType nb = ilaenv<T>(1, "ORGQR", "", nh, nh, nh);
68
69 return max(IndexType(1), nh)*nb;
70 }
71
72 template <typename IndexType, typename MA, typename VTAU, typename VW>
73 void
74 orghr_generic(IndexType iLo,
75 IndexType iHi,
76 GeMatrix<MA> &A,
77 const DenseVector<VTAU> &tau,
78 DenseVector<VW> &work)
79 {
80 using std::max;
81
82 typedef typename GeMatrix<MA>::ElementType T;
83
84 const Underscore<IndexType> _;
85 const IndexType n = A.numRows();
86 const IndexType nh = iHi - iLo;
87 const IndexType nb = ilaenv<T>(1, "ORGQR", "", nh, nh, nh);
88 const IndexType lWorkOpt = max(IndexType(1), nh)*nb;
89
90 //
91 // Apply workspace query if needed
92 //
93 if (work.length()==0) {
94 work.resize(lWorkOpt);
95 }
96 //
97 // Quick return if possible
98 //
99 if (n==0) {
100 work(1) = 1;
101 }
102 //
103 // Shift the vectors which define the elementary reflectors one
104 // column to the right, and set the first ilo and the last n-ihi
105 // rows and columns to those of the unit matrix
106 //
107 for (IndexType j=iHi; j>iLo; --j) {
108 for (IndexType i=1; i<=j-1; ++i) {
109 A(i,j) = T(0);
110 }
111 for (IndexType i=j+1; i<=iHi; ++i) {
112 A(i,j) = A(i,j-1);
113 }
114 for (IndexType i=iHi+1; i<=n; ++i) {
115 A(i,j) = T(0);
116 }
117 }
118 for (IndexType j=1; j<=iLo; ++j) {
119 for (IndexType i=1; i<=n; ++i) {
120 A(i,j) = T(0);
121 }
122 A(j,j) = 1;
123 }
124 for (IndexType j=iHi+1; j<=n; ++j) {
125 for (IndexType i=1; i<=n; ++i) {
126 A(i,j) = T(0);
127 }
128 A(j,j) = T(1);
129 }
130
131 if (nh>0) {
132 //
133 // Generate Q(ilo+1:ihi,ilo+1:ihi)
134 //
135 orgqr(nh, A(_(iLo+1,iHi),_(iLo+1,iHi)), tau(_(iLo,iHi-1)), work);
136 }
137 work(1) = lWorkOpt;
138 }
139
140 //== interface for native lapack ===============================================
141
142 #ifdef CHECK_CXXLAPACK
143
144 template <typename IndexType, typename MA, typename VTAU>
145 IndexType
146 orghr_native_wsq(IndexType iLo,
147 IndexType iHi,
148 const GeMatrix<MA> &A,
149 const DenseVector<VTAU> &tau)
150 {
151 typedef typename GeMatrix<MA>::ElementType T;
152
153 const INTEGER N = A.numRows();
154 const INTEGER ILO = iLo;
155 const INTEGER IHI = iHi;
156 const INTEGER LDA = A.leadingDimension();
157 T WORK;
158 T DUMMY;
159 const INTEGER LWORK = -1;
160 INTEGER INFO;
161
162 if (IsSame<T, DOUBLE>::value) {
163 LAPACK_IMPL(dorghr)(&N,
164 &ILO,
165 &IHI,
166 &DUMMY,
167 &LDA,
168 tau.data(),
169 &WORK,
170 &LWORK,
171 &INFO);
172 } else {
173 ASSERT(0);
174 }
175 ASSERT(INFO==0);
176 return WORK;
177 }
178
179
180 template <typename IndexType, typename MA, typename VTAU, typename VW>
181 void
182 orghr_native(IndexType iLo,
183 IndexType iHi,
184 GeMatrix<MA> &A,
185 const DenseVector<VTAU> &tau,
186 DenseVector<VW> &work)
187 {
188 typedef typename GeMatrix<MA>::ElementType T;
189
190 const INTEGER N = A.numRows();
191 const INTEGER ILO = iLo;
192 const INTEGER IHI = iHi;
193 const INTEGER LDA = A.leadingDimension();
194 const INTEGER LWORK = work.length();
195 INTEGER INFO;
196
197 if (IsSame<T, DOUBLE>::value) {
198 LAPACK_IMPL(dorghr)(&N,
199 &ILO,
200 &IHI,
201 A.data(),
202 &LDA,
203 tau.data(),
204 work.data(),
205 &LWORK,
206 &INFO);
207 } else {
208 ASSERT(0);
209 }
210 ASSERT(INFO==0);
211 }
212
213 #endif // CHECK_CXXLAPACK
214
215 //== public interface ==========================================================
216
217 template <typename IndexType, typename MA, typename VTAU>
218 IndexType
219 orghr_wsq(IndexType iLo,
220 IndexType iHi,
221 const GeMatrix<MA> &A,
222 const DenseVector<VTAU> &tau)
223 {
224 //
225 // Test the input parameters
226 //
227 # ifndef NDEBUG
228 ASSERT(A.numRows()==A.numCols());
229
230 const IndexType n = A.numCols();
231 if (n==0) {
232 ASSERT(iLo==1);
233 ASSERT(iHi==0);
234 } else {
235 ASSERT(1<=iLo);
236 ASSERT(iLo<=iHi);
237 ASSERT(iHi<=n);
238 }
239 ASSERT(tau.length()==(n-1));
240 # endif
241
242 //
243 // Call implementation
244 //
245 IndexType ws = orghr_generic_wsq(iLo, iHi, A, tau);
246
247 # ifdef CHECK_CXXLAPACK
248 //
249 // Compare results
250 //
251 IndexType _ws = orghr_native_wsq(iLo, iHi, A, tau);
252
253 if (ws!=_ws) {
254 std::cerr << "CXXLAPACK: ws = " << ws << std::endl;
255 std::cerr << "F77LAPACK: _ws = " << _ws << std::endl;
256 ASSERT(0);
257 }
258 # endif
259
260 return ws;
261 }
262
263 template <typename IndexType, typename MA, typename VTAU, typename VW>
264 void
265 orghr(IndexType iLo,
266 IndexType iHi,
267 GeMatrix<MA> &A,
268 const DenseVector<VTAU> &tau,
269 DenseVector<VW> &work)
270 {
271 LAPACK_DEBUG_OUT("orghr");
272
273 //
274 // Test the input parameters
275 //
276 # ifndef NDEBUG
277 ASSERT(A.numRows()==A.numCols());
278
279 const IndexType n = A.numCols();
280 if (n==0) {
281 ASSERT(iLo==1);
282 ASSERT(iHi==0);
283 } else {
284 ASSERT(1<=iLo);
285 ASSERT(iLo<=iHi);
286 ASSERT(iHi<=n);
287 }
288 ASSERT(tau.length()==(n-1));
289 ASSERT(work.length()>=iHi-iLo);
290 # endif
291
292 //
293 // Make copies of output arguments
294 //
295 # ifdef CHECK_CXXLAPACK
296 typename GeMatrix<MA>::NoView _A = A;
297 typename DenseVector<VW>::NoView _work = work;
298 # endif
299
300 //
301 // Call implementation
302 //
303 orghr_generic(iLo, iHi, A, tau, work);
304
305 # ifdef CHECK_CXXLAPACK
306 //
307 // Compare results
308 //
309 if (_work.length()==0) {
310 _work.resize(work.length());
311 }
312 orghr_native(iLo, iHi, _A, tau, _work);
313
314 bool failed = false;
315 if (! isIdentical(A, _A, " A", "A_")) {
316 std::cerr << "CXXLAPACK: A = " << A << std::endl;
317 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
318 failed = true;
319 }
320
321 if (! isIdentical(work, _work, " work", "_work")) {
322 std::cerr << "CXXLAPACK: work = " << work << std::endl;
323 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
324 failed = true;
325 }
326
327 if (failed) {
328 std::cerr << "error in: hrd.tcc" << std::endl;
329 ASSERT(0);
330 } else {
331 // std::cerr << "passed: hrd.tcc" << std::endl;
332 }
333 # endif
334 }
335
336 //-- forwarding ----------------------------------------------------------------
337 template <typename IndexType, typename MA, typename VTAU>
338 IndexType
339 orghr_wsq(IndexType iLo,
340 IndexType iHi,
341 const MA &&A,
342 const VTAU &tau)
343 {
344 return orghr_wsq(iLo, iHi, A, tau);
345 }
346
347 template <typename IndexType, typename MA, typename VTAU, typename VW>
348 void
349 orghr(IndexType iLo,
350 IndexType iHi,
351 MA &&A,
352 const VTAU &tau,
353 VW &&work)
354 {
355 orghr(iLo, iHi, A, tau, work);
356 }
357
358 } } // namespace lapack, flens
359
360 #endif // FLENS_LAPACK_EIG_ORGHR_TCC