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 DORMHR( SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C,
36 $ LDC, WORK, LWORK, INFO )
37 *
38 * -- LAPACK routine (version 3.2) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * November 2006
42 */
43
44 #ifndef FLENS_LAPACK_EIG_ORMHR_TCC
45 #define FLENS_LAPACK_EIG_ORMHR_TCC 1
46
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename IndexType, typename MC>
55 IndexType
56 ormhr_generic_wsq(Side side,
57 Transpose trans,
58 IndexType iLo,
59 IndexType iHi,
60 const GeMatrix<MC> &C)
61 {
62 using std::max;
63
64 typedef typename GeMatrix<MC>::ElementType T;
65
66 const IndexType m = C.numRows();
67 const IndexType n = C.numCols();
68 const IndexType nh = iHi - iLo;
69
70 //
71 // nw is the minimum dimension of WORK
72 //
73 const IndexType nw = (side==Left) ? n : m;
74
75 // TODO: implement a better way for setting up the opt string
76 char opt[3];
77 opt[0] = getF77LapackChar(side);
78 opt[1] = getF77LapackChar(trans);
79 opt[2] = 0;
80
81 IndexType nb;
82 if (side==Left) {
83 nb = ilaenv<T>(1, "ORMQR", opt, nh, n, nh);
84 } else {
85 nb = ilaenv<T>(1, "ORMQR", opt, m, nh, nh);
86 }
87 return max(IndexType(1), nw)*nb;
88 }
89
90 template <typename IndexType, typename MA, typename VTAU,
91 typename MC, typename VWORK>
92 void
93 ormhr_generic(Side side,
94 Transpose trans,
95 IndexType iLo,
96 IndexType iHi,
97 GeMatrix<MA> &A,
98 const DenseVector<VTAU> &tau,
99 GeMatrix<MC> &C,
100 DenseVector<VWORK> &work)
101 {
102 using std::max;
103
104 typedef typename GeMatrix<MC>::ElementType T;
105
106 const Underscore<IndexType> _;
107
108 const IndexType m = C.numRows();
109 const IndexType n = C.numCols();
110 const IndexType nh = iHi - iLo;
111
112 //
113 // nw is the minimum dimension of WORK
114 //
115 const IndexType nw = (side==Left) ? n : m;
116
117 // TODO: implement a better way for setting up the opt string
118 char opt[3];
119 opt[0] = getF77LapackChar(side);
120 opt[1] = getF77LapackChar(trans);
121 opt[2] = 0;
122
123 IndexType nb;
124 if (side==Left) {
125 nb = ilaenv<T>(1, "ORMQR", opt, nh, n, nh);
126 } else {
127 nb = ilaenv<T>(1, "ORMQR", opt, m, nh, nh);
128 }
129 IndexType lWorkOpt = max(IndexType(1), nw)*nb;
130
131 //
132 // Apply worksize query
133 //
134 if (work.length()==0) {
135 work.resize(lWorkOpt);
136 }
137
138 //
139 // Quick return if possible
140 //
141 if ((m==0) || (n==0) || (nh==0)) {
142 work(1) = 1;
143 return;
144 }
145
146 const auto _tau = tau(_(iLo,iHi-1));
147 auto _A = A(_(iLo+1,iHi),_(iLo,iHi-1));
148
149 if (side==Left) {
150
151 auto _C = C(_(iLo+1,iHi),_);
152
153 ormqr(Left, trans, _A, _tau, _C, work);
154
155 } else {
156
157 auto _C = C(_,_(iLo+1,iHi));
158
159 ormqr(Right, trans, _A, _tau, _C, work);
160 }
161
162 work(1) = lWorkOpt;
163 }
164
165 //== interface for native lapack ===============================================
166
167 #ifdef CHECK_CXXLAPACK
168
169 template <typename IndexType, typename MC>
170 IndexType
171 ormhr_native_wsq(Side side,
172 Transpose trans,
173 IndexType iLo,
174 IndexType iHi,
175 const GeMatrix<MC> &C)
176 {
177 using std::max;
178
179 typedef typename GeMatrix<MC>::ElementType T;
180
181 IndexType m = C.numRows();
182 IndexType n = C.numCols();
183
184 const char SIDE = getF77LapackChar(side);
185 const char TRANS = getF77LapackChar(trans);
186 const INTEGER M = m;
187 const INTEGER N = n;
188 const INTEGER ILO = iLo;
189 const INTEGER IHI = iHi;
190 const INTEGER LDA = (SIDE=='L') ? max(IndexType(1), m)
191 : max(IndexType(1), n);
192 const INTEGER LDC = C.leadingDimension();
193 T WORK;
194 T DUMMY;
195 const INTEGER LWORK = -1;
196 INTEGER INFO;
197
198 if (IsSame<T,DOUBLE>::value) {
199 LAPACK_IMPL(dormhr)(&SIDE,
200 &TRANS,
201 &M,
202 &N,
203 &ILO,
204 &IHI,
205 &DUMMY,
206 &LDA,
207 &DUMMY,
208 &DUMMY,
209 &LDC,
210 &WORK,
211 &LWORK,
212 &INFO);
213 } else {
214 ASSERT(0);
215 }
216 ASSERT(INFO==0);
217 return WORK;
218 }
219
220 template <typename IndexType, typename MA, typename VTAU,
221 typename MC, typename VWORK>
222 void
223 ormhr_native(Side side,
224 Transpose trans,
225 IndexType iLo,
226 IndexType iHi,
227 GeMatrix<MA> &A,
228 const DenseVector<VTAU> &tau,
229 GeMatrix<MC> &C,
230 DenseVector<VWORK> &work)
231 {
232 typedef typename GeMatrix<MC>::ElementType T;
233
234 const char SIDE = getF77LapackChar(side);
235 const char TRANS = getF77LapackChar(trans);
236 const INTEGER M = C.numRows();
237 const INTEGER N = C.numCols();
238 const INTEGER ILO = iLo;
239 const INTEGER IHI = iHi;
240 const INTEGER LDA = A.leadingDimension();
241 const INTEGER LDC = C.leadingDimension();
242 const INTEGER LWORK = work.length();
243 INTEGER INFO;
244
245 if (IsSame<T,DOUBLE>::value) {
246 LAPACK_IMPL(dormhr)(&SIDE,
247 &TRANS,
248 &M,
249 &N,
250 &ILO,
251 &IHI,
252 A.data(),
253 &LDA,
254 tau.data(),
255 C.data(),
256 &LDC,
257 work.data(),
258 &LWORK,
259 &INFO);
260 } else {
261 ASSERT(0);
262 }
263 ASSERT(INFO==0);
264 }
265
266 #endif // CHECK_CXXLAPACK
267
268 //== public interface ==========================================================
269
270 template <typename IndexType, typename MC>
271 IndexType
272 ormhr_wsq(Side side,
273 Transpose trans,
274 IndexType iLo,
275 IndexType iHi,
276 const GeMatrix<MC> &C)
277 {
278 LAPACK_DEBUG_OUT("ormhr_wsq");
279
280 //
281 // Test the input parameters
282 //
283 # ifndef NDEBUG
284 const IndexType m = C.numRows();
285 const IndexType n = C.numCols();
286
287 if (side==Left) {
288 if (m==0) {
289 ASSERT(iLo==1);
290 ASSERT(iHi==0);
291 } else {
292 ASSERT(1<=iLo);
293 ASSERT(iLo<=iHi);
294 ASSERT(iHi<=m);
295 }
296 }
297 if (side==Right) {
298 if (n==0) {
299 ASSERT(iLo==1);
300 ASSERT(iHi==0);
301 } else {
302 ASSERT(1<=iLo);
303 ASSERT(iLo<=iHi);
304 ASSERT(iHi<=n);
305 }
306 }
307 # endif
308
309 //
310 // Call implementation
311 //
312 IndexType ws = ormhr_generic_wsq(side, trans, iLo, iHi, C);
313
314 //
315 // Compare results
316 //
317 # ifdef CHECK_CXXLAPACK
318 IndexType _ws = ormhr_native_wsq(side, trans, iLo, iHi, C);
319
320 if (ws!=_ws) {
321 std::cerr << "CXXLAPACK: ws = " << ws << std::endl;
322 std::cerr << "F77LAPACK: _ws = " << _ws << std::endl;
323 ASSERT(0);
324 }
325 # endif
326 return ws;
327 }
328
329 template <typename IndexType, typename MA, typename VTAU,
330 typename MC, typename VWORK>
331 void
332 ormhr(Side side,
333 Transpose trans,
334 IndexType iLo,
335 IndexType iHi,
336 GeMatrix<MA> &A,
337 const DenseVector<VTAU> &tau,
338 GeMatrix<MC> &C,
339 DenseVector<VWORK> &work)
340 {
341 LAPACK_DEBUG_OUT("ormhr");
342
343 using std::max;
344
345 //
346 // Test the input parameters
347 //
348 # ifndef NDEBUG
349 const IndexType m = C.numRows();
350 const IndexType n = C.numCols();
351
352 ASSERT(A.firstRow()==1);
353 ASSERT(A.firstCol()==1);
354 ASSERT(tau.firstIndex()==1);
355 if (side==Left) {
356 if (m==0) {
357 ASSERT(iLo==1);
358 ASSERT(iHi==0);
359 } else {
360 ASSERT(1<=iLo);
361 ASSERT(iLo<=iHi);
362 ASSERT(iHi<=m);
363 }
364 }
365 if (side==Right) {
366 if (n==0) {
367 ASSERT(iLo==1);
368 ASSERT(iHi==0);
369 } else {
370 ASSERT(1<=iLo);
371 ASSERT(iLo<=iHi);
372 ASSERT(iHi<=n);
373 }
374 }
375 ASSERT((side==Left && A.numCols()==m)
376 || (side==Right && A.numCols()==n));
377 ASSERT((side==Left && tau.length()==(m-1))
378 || (side==Right && tau.length()==(n-1)));
379 ASSERT((side==Left && work.length()>=max(IndexType(1),n))
380 || (side==Right && work.length()>=max(IndexType(1),m)));
381 # endif
382
383 //
384 // Make copies of output arguments
385 //
386 # ifdef CHECK_CXXLAPACK
387 typename GeMatrix<MA>::NoView A_org = A;
388 typename GeMatrix<MC>::NoView C_org = C;
389 typename DenseVector<VWORK>::NoView work_org = work;
390 # endif
391
392 //
393 // Call implementation
394 //
395 ormhr_generic(side, trans, iLo, iHi, A, tau, C, work);
396
397 # ifdef CHECK_CXXLAPACK
398 //
399 // Make copies of results computed by the generic implementation
400 //
401 typename GeMatrix<MA>::NoView A_generic = A;
402 typename GeMatrix<MC>::NoView C_generic = C;
403 typename DenseVector<VWORK>::NoView work_generic = work;
404
405 //
406 // restore output arguments
407 //
408 A = A_org;
409 C = C_org;
410 work = work_org;
411
412 //
413 // Compare generic results with results from the native implementation
414 //
415 ormhr_native(side, trans, iLo, iHi, A, tau, C, work);
416
417 bool failed = false;
418 if (! isIdentical(A_generic, A, "A_generic", "A")) {
419 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
420 std::cerr << "F77LAPACK: A = " << A << std::endl;
421 failed = true;
422 }
423 if (! isIdentical(C_generic, C, "C_generic", "C")) {
424 std::cerr << "CXXLAPACK: C_generic = " << C_generic << std::endl;
425 std::cerr << "F77LAPACK: C = " << C << std::endl;
426 failed = true;
427 }
428 if (! isIdentical(work_generic, work, "work_generic", "work")) {
429 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
430 std::cerr << "F77LAPACK: work = " << work << std::endl;
431 failed = true;
432 }
433
434 if (failed) {
435 std::cerr << "error in: ormhr.tcc" << std::endl;
436 ASSERT(0);
437 } else {
438 // std::cerr << "passed: ormhr.tcc" << std::endl;
439 }
440 # endif
441 }
442
443 //-- forwarding ----------------------------------------------------------------
444 template <typename IndexType, typename MA, typename VTAU, typename MC>
445 IndexType
446 ormhr_wsq(Side side,
447 Transpose trans,
448 IndexType iLo,
449 IndexType iHi,
450 const MC &&C)
451 {
452 CHECKPOINT_ENTER;
453 const IndexType info = ormhr_wsq(side, trans, iLo, iHi, C);
454 CHECKPOINT_LEAVE;
455
456 return info;
457 }
458
459 template <typename IndexType, typename MA, typename VTAU,
460 typename MC, typename VWORK>
461 void
462 ormhr(Side side,
463 Transpose trans,
464 IndexType iLo,
465 IndexType iHi,
466 MA &&A,
467 const VTAU &tau,
468 MC &&C,
469 VWORK &&work)
470 {
471 CHECKPOINT_ENTER;
472 ormhr(side, trans, iLo, iHi, A, tau, C, work);
473 CHECKPOINT_LEAVE;
474 }
475
476 } } // namespace lapack, flens
477
478 #endif // FLENS_LAPACK_EIG_ORMHR_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 DORMHR( SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C,
36 $ LDC, WORK, LWORK, INFO )
37 *
38 * -- LAPACK routine (version 3.2) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * November 2006
42 */
43
44 #ifndef FLENS_LAPACK_EIG_ORMHR_TCC
45 #define FLENS_LAPACK_EIG_ORMHR_TCC 1
46
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename IndexType, typename MC>
55 IndexType
56 ormhr_generic_wsq(Side side,
57 Transpose trans,
58 IndexType iLo,
59 IndexType iHi,
60 const GeMatrix<MC> &C)
61 {
62 using std::max;
63
64 typedef typename GeMatrix<MC>::ElementType T;
65
66 const IndexType m = C.numRows();
67 const IndexType n = C.numCols();
68 const IndexType nh = iHi - iLo;
69
70 //
71 // nw is the minimum dimension of WORK
72 //
73 const IndexType nw = (side==Left) ? n : m;
74
75 // TODO: implement a better way for setting up the opt string
76 char opt[3];
77 opt[0] = getF77LapackChar(side);
78 opt[1] = getF77LapackChar(trans);
79 opt[2] = 0;
80
81 IndexType nb;
82 if (side==Left) {
83 nb = ilaenv<T>(1, "ORMQR", opt, nh, n, nh);
84 } else {
85 nb = ilaenv<T>(1, "ORMQR", opt, m, nh, nh);
86 }
87 return max(IndexType(1), nw)*nb;
88 }
89
90 template <typename IndexType, typename MA, typename VTAU,
91 typename MC, typename VWORK>
92 void
93 ormhr_generic(Side side,
94 Transpose trans,
95 IndexType iLo,
96 IndexType iHi,
97 GeMatrix<MA> &A,
98 const DenseVector<VTAU> &tau,
99 GeMatrix<MC> &C,
100 DenseVector<VWORK> &work)
101 {
102 using std::max;
103
104 typedef typename GeMatrix<MC>::ElementType T;
105
106 const Underscore<IndexType> _;
107
108 const IndexType m = C.numRows();
109 const IndexType n = C.numCols();
110 const IndexType nh = iHi - iLo;
111
112 //
113 // nw is the minimum dimension of WORK
114 //
115 const IndexType nw = (side==Left) ? n : m;
116
117 // TODO: implement a better way for setting up the opt string
118 char opt[3];
119 opt[0] = getF77LapackChar(side);
120 opt[1] = getF77LapackChar(trans);
121 opt[2] = 0;
122
123 IndexType nb;
124 if (side==Left) {
125 nb = ilaenv<T>(1, "ORMQR", opt, nh, n, nh);
126 } else {
127 nb = ilaenv<T>(1, "ORMQR", opt, m, nh, nh);
128 }
129 IndexType lWorkOpt = max(IndexType(1), nw)*nb;
130
131 //
132 // Apply worksize query
133 //
134 if (work.length()==0) {
135 work.resize(lWorkOpt);
136 }
137
138 //
139 // Quick return if possible
140 //
141 if ((m==0) || (n==0) || (nh==0)) {
142 work(1) = 1;
143 return;
144 }
145
146 const auto _tau = tau(_(iLo,iHi-1));
147 auto _A = A(_(iLo+1,iHi),_(iLo,iHi-1));
148
149 if (side==Left) {
150
151 auto _C = C(_(iLo+1,iHi),_);
152
153 ormqr(Left, trans, _A, _tau, _C, work);
154
155 } else {
156
157 auto _C = C(_,_(iLo+1,iHi));
158
159 ormqr(Right, trans, _A, _tau, _C, work);
160 }
161
162 work(1) = lWorkOpt;
163 }
164
165 //== interface for native lapack ===============================================
166
167 #ifdef CHECK_CXXLAPACK
168
169 template <typename IndexType, typename MC>
170 IndexType
171 ormhr_native_wsq(Side side,
172 Transpose trans,
173 IndexType iLo,
174 IndexType iHi,
175 const GeMatrix<MC> &C)
176 {
177 using std::max;
178
179 typedef typename GeMatrix<MC>::ElementType T;
180
181 IndexType m = C.numRows();
182 IndexType n = C.numCols();
183
184 const char SIDE = getF77LapackChar(side);
185 const char TRANS = getF77LapackChar(trans);
186 const INTEGER M = m;
187 const INTEGER N = n;
188 const INTEGER ILO = iLo;
189 const INTEGER IHI = iHi;
190 const INTEGER LDA = (SIDE=='L') ? max(IndexType(1), m)
191 : max(IndexType(1), n);
192 const INTEGER LDC = C.leadingDimension();
193 T WORK;
194 T DUMMY;
195 const INTEGER LWORK = -1;
196 INTEGER INFO;
197
198 if (IsSame<T,DOUBLE>::value) {
199 LAPACK_IMPL(dormhr)(&SIDE,
200 &TRANS,
201 &M,
202 &N,
203 &ILO,
204 &IHI,
205 &DUMMY,
206 &LDA,
207 &DUMMY,
208 &DUMMY,
209 &LDC,
210 &WORK,
211 &LWORK,
212 &INFO);
213 } else {
214 ASSERT(0);
215 }
216 ASSERT(INFO==0);
217 return WORK;
218 }
219
220 template <typename IndexType, typename MA, typename VTAU,
221 typename MC, typename VWORK>
222 void
223 ormhr_native(Side side,
224 Transpose trans,
225 IndexType iLo,
226 IndexType iHi,
227 GeMatrix<MA> &A,
228 const DenseVector<VTAU> &tau,
229 GeMatrix<MC> &C,
230 DenseVector<VWORK> &work)
231 {
232 typedef typename GeMatrix<MC>::ElementType T;
233
234 const char SIDE = getF77LapackChar(side);
235 const char TRANS = getF77LapackChar(trans);
236 const INTEGER M = C.numRows();
237 const INTEGER N = C.numCols();
238 const INTEGER ILO = iLo;
239 const INTEGER IHI = iHi;
240 const INTEGER LDA = A.leadingDimension();
241 const INTEGER LDC = C.leadingDimension();
242 const INTEGER LWORK = work.length();
243 INTEGER INFO;
244
245 if (IsSame<T,DOUBLE>::value) {
246 LAPACK_IMPL(dormhr)(&SIDE,
247 &TRANS,
248 &M,
249 &N,
250 &ILO,
251 &IHI,
252 A.data(),
253 &LDA,
254 tau.data(),
255 C.data(),
256 &LDC,
257 work.data(),
258 &LWORK,
259 &INFO);
260 } else {
261 ASSERT(0);
262 }
263 ASSERT(INFO==0);
264 }
265
266 #endif // CHECK_CXXLAPACK
267
268 //== public interface ==========================================================
269
270 template <typename IndexType, typename MC>
271 IndexType
272 ormhr_wsq(Side side,
273 Transpose trans,
274 IndexType iLo,
275 IndexType iHi,
276 const GeMatrix<MC> &C)
277 {
278 LAPACK_DEBUG_OUT("ormhr_wsq");
279
280 //
281 // Test the input parameters
282 //
283 # ifndef NDEBUG
284 const IndexType m = C.numRows();
285 const IndexType n = C.numCols();
286
287 if (side==Left) {
288 if (m==0) {
289 ASSERT(iLo==1);
290 ASSERT(iHi==0);
291 } else {
292 ASSERT(1<=iLo);
293 ASSERT(iLo<=iHi);
294 ASSERT(iHi<=m);
295 }
296 }
297 if (side==Right) {
298 if (n==0) {
299 ASSERT(iLo==1);
300 ASSERT(iHi==0);
301 } else {
302 ASSERT(1<=iLo);
303 ASSERT(iLo<=iHi);
304 ASSERT(iHi<=n);
305 }
306 }
307 # endif
308
309 //
310 // Call implementation
311 //
312 IndexType ws = ormhr_generic_wsq(side, trans, iLo, iHi, C);
313
314 //
315 // Compare results
316 //
317 # ifdef CHECK_CXXLAPACK
318 IndexType _ws = ormhr_native_wsq(side, trans, iLo, iHi, C);
319
320 if (ws!=_ws) {
321 std::cerr << "CXXLAPACK: ws = " << ws << std::endl;
322 std::cerr << "F77LAPACK: _ws = " << _ws << std::endl;
323 ASSERT(0);
324 }
325 # endif
326 return ws;
327 }
328
329 template <typename IndexType, typename MA, typename VTAU,
330 typename MC, typename VWORK>
331 void
332 ormhr(Side side,
333 Transpose trans,
334 IndexType iLo,
335 IndexType iHi,
336 GeMatrix<MA> &A,
337 const DenseVector<VTAU> &tau,
338 GeMatrix<MC> &C,
339 DenseVector<VWORK> &work)
340 {
341 LAPACK_DEBUG_OUT("ormhr");
342
343 using std::max;
344
345 //
346 // Test the input parameters
347 //
348 # ifndef NDEBUG
349 const IndexType m = C.numRows();
350 const IndexType n = C.numCols();
351
352 ASSERT(A.firstRow()==1);
353 ASSERT(A.firstCol()==1);
354 ASSERT(tau.firstIndex()==1);
355 if (side==Left) {
356 if (m==0) {
357 ASSERT(iLo==1);
358 ASSERT(iHi==0);
359 } else {
360 ASSERT(1<=iLo);
361 ASSERT(iLo<=iHi);
362 ASSERT(iHi<=m);
363 }
364 }
365 if (side==Right) {
366 if (n==0) {
367 ASSERT(iLo==1);
368 ASSERT(iHi==0);
369 } else {
370 ASSERT(1<=iLo);
371 ASSERT(iLo<=iHi);
372 ASSERT(iHi<=n);
373 }
374 }
375 ASSERT((side==Left && A.numCols()==m)
376 || (side==Right && A.numCols()==n));
377 ASSERT((side==Left && tau.length()==(m-1))
378 || (side==Right && tau.length()==(n-1)));
379 ASSERT((side==Left && work.length()>=max(IndexType(1),n))
380 || (side==Right && work.length()>=max(IndexType(1),m)));
381 # endif
382
383 //
384 // Make copies of output arguments
385 //
386 # ifdef CHECK_CXXLAPACK
387 typename GeMatrix<MA>::NoView A_org = A;
388 typename GeMatrix<MC>::NoView C_org = C;
389 typename DenseVector<VWORK>::NoView work_org = work;
390 # endif
391
392 //
393 // Call implementation
394 //
395 ormhr_generic(side, trans, iLo, iHi, A, tau, C, work);
396
397 # ifdef CHECK_CXXLAPACK
398 //
399 // Make copies of results computed by the generic implementation
400 //
401 typename GeMatrix<MA>::NoView A_generic = A;
402 typename GeMatrix<MC>::NoView C_generic = C;
403 typename DenseVector<VWORK>::NoView work_generic = work;
404
405 //
406 // restore output arguments
407 //
408 A = A_org;
409 C = C_org;
410 work = work_org;
411
412 //
413 // Compare generic results with results from the native implementation
414 //
415 ormhr_native(side, trans, iLo, iHi, A, tau, C, work);
416
417 bool failed = false;
418 if (! isIdentical(A_generic, A, "A_generic", "A")) {
419 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
420 std::cerr << "F77LAPACK: A = " << A << std::endl;
421 failed = true;
422 }
423 if (! isIdentical(C_generic, C, "C_generic", "C")) {
424 std::cerr << "CXXLAPACK: C_generic = " << C_generic << std::endl;
425 std::cerr << "F77LAPACK: C = " << C << std::endl;
426 failed = true;
427 }
428 if (! isIdentical(work_generic, work, "work_generic", "work")) {
429 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
430 std::cerr << "F77LAPACK: work = " << work << std::endl;
431 failed = true;
432 }
433
434 if (failed) {
435 std::cerr << "error in: ormhr.tcc" << std::endl;
436 ASSERT(0);
437 } else {
438 // std::cerr << "passed: ormhr.tcc" << std::endl;
439 }
440 # endif
441 }
442
443 //-- forwarding ----------------------------------------------------------------
444 template <typename IndexType, typename MA, typename VTAU, typename MC>
445 IndexType
446 ormhr_wsq(Side side,
447 Transpose trans,
448 IndexType iLo,
449 IndexType iHi,
450 const MC &&C)
451 {
452 CHECKPOINT_ENTER;
453 const IndexType info = ormhr_wsq(side, trans, iLo, iHi, C);
454 CHECKPOINT_LEAVE;
455
456 return info;
457 }
458
459 template <typename IndexType, typename MA, typename VTAU,
460 typename MC, typename VWORK>
461 void
462 ormhr(Side side,
463 Transpose trans,
464 IndexType iLo,
465 IndexType iHi,
466 MA &&A,
467 const VTAU &tau,
468 MC &&C,
469 VWORK &&work)
470 {
471 CHECKPOINT_ENTER;
472 ormhr(side, trans, iLo, iHi, A, tau, C, work);
473 CHECKPOINT_LEAVE;
474 }
475
476 } } // namespace lapack, flens
477
478 #endif // FLENS_LAPACK_EIG_ORMHR_TCC