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 SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO )
35 *
36 * -- LAPACK routine (version 3.3.1) --
37 * -- LAPACK is a software package provided by Univ. of Tennessee, --
38 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
39 * -- April 2011 --
40 */
41
42 #ifndef FLENS_LAPACK_QR_QR2_TCC
43 #define FLENS_LAPACK_QR_QR2_TCC 1
44
45 #include <flens/blas/blas.h>
46 #include <flens/lapack/lapack.h>
47
48 namespace flens { namespace lapack {
49
50 //== generic lapack implementation =============================================
51
52 template <typename MA, typename VTAU, typename VWORK>
53 void
54 qr2_generic(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
55 {
56 typedef typename GeMatrix<MA>::IndexType IndexType;
57 typedef typename GeMatrix<MA>::ElementType T;
58
59 const Underscore<IndexType> _;
60
61 const IndexType m = A.numRows();
62 const IndexType n = A.numCols();
63 const IndexType k = std::min(m, n);
64
65 for (IndexType i=1; i<=k; ++i) {
66 //
67 // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
68 //
69 larfg(m-i+1, A(i,i), A(_(std::min(i+1,m),m), i), tau(i));
70
71 if (i<n) {
72 //
73 // Apply H(i) to A(i:m,i+1:n) from the left
74 //
75 const T Aii = A(i,i);
76 A(i,i) = T(1);
77 auto _work = work(_(1, n-i));
78 larf(Left, A(_(i,m),i), tau(i), A(_(i,m), _(i+1,n)), _work);
79 A(i,i) = Aii;
80 }
81 }
82 }
83
84 //== interface for native lapack ===============================================
85
86 #ifdef CHECK_CXXLAPACK
87
88 template <typename MA, typename VTAU, typename VWORK>
89 void
90 qr2_native(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
91 {
92 typedef typename GeMatrix<MA>::ElementType T;
93
94 const INTEGER M = A.numRows();
95 const INTEGER N = A.numCols();
96 const INTEGER LDA = A.leadingDimension();
97 INTEGER INFO;
98
99 if (IsSame<T, DOUBLE>::value) {
100 LAPACK_IMPL(dgeqr2)(&M, &N, A.data(), &LDA,
101 tau.data(), work.data(),
102 &INFO);
103 assert(INFO==0);
104 } else {
105 ASSERT(0);
106 }
107 }
108
109 #endif // CHECK_CXXLAPACK
110
111 //== public interface ==========================================================
112
113 template <typename MA, typename VTAU, typename VWORK>
114 void
115 qr2(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
116 {
117 typedef typename GeMatrix<MA>::IndexType IndexType;
118
119 //
120 // Test the input parameters
121 //
122 # ifndef NDEBUG
123 ASSERT(A.firstRow()==1);
124 ASSERT(A.firstCol()==1);
125 ASSERT(tau.firstIndex()==1);
126 ASSERT(work.firstIndex()==1);
127
128 const IndexType m = A.numRows();
129 const IndexType n = A.numCols();
130 const IndexType k = std::min(m, n);
131
132 ASSERT(tau.length()==k);
133 ASSERT(work.length()==n);
134 # endif
135
136 # ifdef CHECK_CXXLAPACK
137 //
138 // Make copies of output arguments
139 //
140 typename GeMatrix<MA>::NoView A_org = A;
141 typename DenseVector<VTAU>::NoView tau_org = tau;
142 typename DenseVector<VTAU>::NoView work_org = work;
143 # endif
144
145 //
146 // Call implementation
147 //
148 qr2_generic(A, tau, work);
149
150 # ifdef CHECK_CXXLAPACK
151 //
152 // Restore output arguments
153 //
154 typename GeMatrix<MA>::NoView A_generic = A;
155 typename DenseVector<VTAU>::NoView tau_generic = tau;
156 typename DenseVector<VTAU>::NoView work_generic = work;
157
158 A = A_org;
159 tau = tau_org;
160 work = work_org;
161 //
162 // Compare results
163 //
164 qr2_native(A, tau, work);
165
166 bool failed = false;
167 if (! isIdentical(A_generic, A, "A_generic", "A")) {
168 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
169 std::cerr << "F77LAPACK: A = " << A << std::endl;
170 failed = true;
171 }
172
173 if (! isIdentical(tau_generic, tau, "tau_generic", "tau")) {
174 std::cerr << "CXXLAPACK: tau_generic = " << tau_generic << std::endl;
175 std::cerr << "F77LAPACK: tau = " << tau << std::endl;
176 failed = true;
177 }
178
179 if (! isIdentical(work_generic, work, "work_generic", "work")) {
180 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
181 std::cerr << "F77LAPACK: work = " << work << std::endl;
182 failed = true;
183 }
184
185 if (failed) {
186 ASSERT(0);
187 }
188 # endif
189 }
190
191 //-- forwarding ----------------------------------------------------------------
192 template <typename MA, typename VTAU, typename VWORK>
193 void
194 qr2(MA &&A, VTAU &&tau, VWORK &&work)
195 {
196 CHECKPOINT_ENTER;
197 qr2(A, tau, work);
198 CHECKPOINT_LEAVE;
199 }
200
201 } } // namespace lapack, flens
202
203 #endif // FLENS_LAPACK_QR_QR2_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 SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO )
35 *
36 * -- LAPACK routine (version 3.3.1) --
37 * -- LAPACK is a software package provided by Univ. of Tennessee, --
38 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
39 * -- April 2011 --
40 */
41
42 #ifndef FLENS_LAPACK_QR_QR2_TCC
43 #define FLENS_LAPACK_QR_QR2_TCC 1
44
45 #include <flens/blas/blas.h>
46 #include <flens/lapack/lapack.h>
47
48 namespace flens { namespace lapack {
49
50 //== generic lapack implementation =============================================
51
52 template <typename MA, typename VTAU, typename VWORK>
53 void
54 qr2_generic(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
55 {
56 typedef typename GeMatrix<MA>::IndexType IndexType;
57 typedef typename GeMatrix<MA>::ElementType T;
58
59 const Underscore<IndexType> _;
60
61 const IndexType m = A.numRows();
62 const IndexType n = A.numCols();
63 const IndexType k = std::min(m, n);
64
65 for (IndexType i=1; i<=k; ++i) {
66 //
67 // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
68 //
69 larfg(m-i+1, A(i,i), A(_(std::min(i+1,m),m), i), tau(i));
70
71 if (i<n) {
72 //
73 // Apply H(i) to A(i:m,i+1:n) from the left
74 //
75 const T Aii = A(i,i);
76 A(i,i) = T(1);
77 auto _work = work(_(1, n-i));
78 larf(Left, A(_(i,m),i), tau(i), A(_(i,m), _(i+1,n)), _work);
79 A(i,i) = Aii;
80 }
81 }
82 }
83
84 //== interface for native lapack ===============================================
85
86 #ifdef CHECK_CXXLAPACK
87
88 template <typename MA, typename VTAU, typename VWORK>
89 void
90 qr2_native(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
91 {
92 typedef typename GeMatrix<MA>::ElementType T;
93
94 const INTEGER M = A.numRows();
95 const INTEGER N = A.numCols();
96 const INTEGER LDA = A.leadingDimension();
97 INTEGER INFO;
98
99 if (IsSame<T, DOUBLE>::value) {
100 LAPACK_IMPL(dgeqr2)(&M, &N, A.data(), &LDA,
101 tau.data(), work.data(),
102 &INFO);
103 assert(INFO==0);
104 } else {
105 ASSERT(0);
106 }
107 }
108
109 #endif // CHECK_CXXLAPACK
110
111 //== public interface ==========================================================
112
113 template <typename MA, typename VTAU, typename VWORK>
114 void
115 qr2(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
116 {
117 typedef typename GeMatrix<MA>::IndexType IndexType;
118
119 //
120 // Test the input parameters
121 //
122 # ifndef NDEBUG
123 ASSERT(A.firstRow()==1);
124 ASSERT(A.firstCol()==1);
125 ASSERT(tau.firstIndex()==1);
126 ASSERT(work.firstIndex()==1);
127
128 const IndexType m = A.numRows();
129 const IndexType n = A.numCols();
130 const IndexType k = std::min(m, n);
131
132 ASSERT(tau.length()==k);
133 ASSERT(work.length()==n);
134 # endif
135
136 # ifdef CHECK_CXXLAPACK
137 //
138 // Make copies of output arguments
139 //
140 typename GeMatrix<MA>::NoView A_org = A;
141 typename DenseVector<VTAU>::NoView tau_org = tau;
142 typename DenseVector<VTAU>::NoView work_org = work;
143 # endif
144
145 //
146 // Call implementation
147 //
148 qr2_generic(A, tau, work);
149
150 # ifdef CHECK_CXXLAPACK
151 //
152 // Restore output arguments
153 //
154 typename GeMatrix<MA>::NoView A_generic = A;
155 typename DenseVector<VTAU>::NoView tau_generic = tau;
156 typename DenseVector<VTAU>::NoView work_generic = work;
157
158 A = A_org;
159 tau = tau_org;
160 work = work_org;
161 //
162 // Compare results
163 //
164 qr2_native(A, tau, work);
165
166 bool failed = false;
167 if (! isIdentical(A_generic, A, "A_generic", "A")) {
168 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
169 std::cerr << "F77LAPACK: A = " << A << std::endl;
170 failed = true;
171 }
172
173 if (! isIdentical(tau_generic, tau, "tau_generic", "tau")) {
174 std::cerr << "CXXLAPACK: tau_generic = " << tau_generic << std::endl;
175 std::cerr << "F77LAPACK: tau = " << tau << std::endl;
176 failed = true;
177 }
178
179 if (! isIdentical(work_generic, work, "work_generic", "work")) {
180 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
181 std::cerr << "F77LAPACK: work = " << work << std::endl;
182 failed = true;
183 }
184
185 if (failed) {
186 ASSERT(0);
187 }
188 # endif
189 }
190
191 //-- forwarding ----------------------------------------------------------------
192 template <typename MA, typename VTAU, typename VWORK>
193 void
194 qr2(MA &&A, VTAU &&tau, VWORK &&work)
195 {
196 CHECKPOINT_ENTER;
197 qr2(A, tau, work);
198 CHECKPOINT_LEAVE;
199 }
200
201 } } // namespace lapack, flens
202
203 #endif // FLENS_LAPACK_QR_QR2_TCC