1 /*
2 * Copyright (c) 2012, 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 DGELQ2( M, N, 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_LQ_LQ2_TCC
44 #define FLENS_LAPACK_LQ_LQ2_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 MA, typename VTAU, typename VWORK>
54 void
55 lq2_generic(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
56 {
57 typedef typename GeMatrix<MA>::ElementType ElementType;
58 typedef typename GeMatrix<MA>::IndexType IndexType;
59
60 const Underscore<IndexType> _;
61
62 const IndexType m = A.numRows();
63 const IndexType n = A.numCols();
64 const IndexType k = std::min(m, n);
65
66 const ElementType One(1);
67
68 for (IndexType i=1; i<=k; ++i) {
69 //
70 // Generate elementary reflector H(i) to annihilate A(i,i+1:n)
71 //
72 larfg(n-i+1, A(i,i), A(i,_(std::min(i+1,n),n)), tau(i));
73
74 if (i<m) {
75 //
76 // Apply H(i) to A(i+1:m,i:n) from the right
77 //
78 const ElementType Aii = A(i,i);
79 A(i,i) = One;
80 auto _work = work(_(1, m-i));
81 larf(Right, A(i,_(i,n)), tau(i), A(_(i+1,m), _(i,n)), _work);
82 A(i,i) = Aii;
83 }
84 }
85 }
86
87 //== interface for native lapack ===============================================
88
89 #ifdef CHECK_CXXLAPACK
90
91 template <typename MA, typename VTAU, typename VWORK>
92 void
93 lq2_native(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
94 {
95 typedef typename GeMatrix<MA>::ElementType T;
96
97 const INTEGER M = A.numRows();
98 const INTEGER N = A.numCols();
99 const INTEGER LDA = A.leadingDimension();
100 INTEGER INFO;
101
102 if (IsSame<T, DOUBLE>::value) {
103 LAPACK_IMPL(dgelq2)(&M, &N, A.data(), &LDA,
104 tau.data(), work.data(),
105 &INFO);
106 assert(INFO==0);
107 } else {
108 ASSERT(0);
109 }
110 }
111
112 #endif // CHECK_CXXLAPACK
113
114 //== public interface ==========================================================
115
116 template <typename MA, typename VTAU, typename VWORK>
117 void
118 lq2(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
119 {
120 typedef typename GeMatrix<MA>::IndexType IndexType;
121
122 //
123 // Test the input parameters
124 //
125 # ifndef NDEBUG
126 ASSERT(A.firstRow()==1);
127 ASSERT(A.firstCol()==1);
128 ASSERT(tau.firstIndex()==1);
129 ASSERT(work.firstIndex()==1);
130
131 const IndexType m = A.numRows();
132 const IndexType n = A.numCols();
133 const IndexType k = std::min(m, n);
134
135 ASSERT(tau.length()>=k);
136 ASSERT(work.length()>=m);
137 # endif
138
139 # ifdef CHECK_CXXLAPACK
140 //
141 // Make copies of output arguments
142 //
143 typename GeMatrix<MA>::NoView A_org = A;
144 typename DenseVector<VTAU>::NoView tau_org = tau;
145 typename DenseVector<VTAU>::NoView work_org = work;
146 # endif
147
148 //
149 // Call implementation
150 //
151 lq2_generic(A, tau, work);
152
153 # ifdef CHECK_CXXLAPACK
154 //
155 // Restore output arguments
156 //
157 typename GeMatrix<MA>::NoView A_generic = A;
158 typename DenseVector<VTAU>::NoView tau_generic = tau;
159 typename DenseVector<VTAU>::NoView work_generic = work;
160
161 A = A_org;
162 tau = tau_org;
163 work = work_org;
164 //
165 // Compare results
166 //
167 lq2_native(A, tau, work);
168
169 bool failed = false;
170 if (! isIdentical(A_generic, A, "A_generic", "A")) {
171 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
172 std::cerr << "F77LAPACK: A = " << A << std::endl;
173 failed = true;
174 }
175
176 if (! isIdentical(tau_generic, tau, "tau_generic", "tau")) {
177 std::cerr << "CXXLAPACK: tau_generic = " << tau_generic << std::endl;
178 std::cerr << "F77LAPACK: tau = " << tau << std::endl;
179 failed = true;
180 }
181
182 if (! isIdentical(work_generic, work, "work_generic", "work")) {
183 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
184 std::cerr << "F77LAPACK: work = " << work << std::endl;
185 failed = true;
186 }
187
188 if (failed) {
189 ASSERT(0);
190 }
191 # endif
192 }
193
194 //-- forwarding ----------------------------------------------------------------
195 template <typename MA, typename VTAU, typename VWORK>
196 void
197 lq2(MA &&A, VTAU &&tau, VWORK &&work)
198 {
199 CHECKPOINT_ENTER;
200 lq2(A, tau, work);
201 CHECKPOINT_LEAVE;
202 }
203
204 } } // namespace lapack, flens
205
206 #endif // FLENS_LAPACK_LQ_LQ2_TCC
2 * Copyright (c) 2012, 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 DGELQ2( M, N, 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_LQ_LQ2_TCC
44 #define FLENS_LAPACK_LQ_LQ2_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 MA, typename VTAU, typename VWORK>
54 void
55 lq2_generic(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
56 {
57 typedef typename GeMatrix<MA>::ElementType ElementType;
58 typedef typename GeMatrix<MA>::IndexType IndexType;
59
60 const Underscore<IndexType> _;
61
62 const IndexType m = A.numRows();
63 const IndexType n = A.numCols();
64 const IndexType k = std::min(m, n);
65
66 const ElementType One(1);
67
68 for (IndexType i=1; i<=k; ++i) {
69 //
70 // Generate elementary reflector H(i) to annihilate A(i,i+1:n)
71 //
72 larfg(n-i+1, A(i,i), A(i,_(std::min(i+1,n),n)), tau(i));
73
74 if (i<m) {
75 //
76 // Apply H(i) to A(i+1:m,i:n) from the right
77 //
78 const ElementType Aii = A(i,i);
79 A(i,i) = One;
80 auto _work = work(_(1, m-i));
81 larf(Right, A(i,_(i,n)), tau(i), A(_(i+1,m), _(i,n)), _work);
82 A(i,i) = Aii;
83 }
84 }
85 }
86
87 //== interface for native lapack ===============================================
88
89 #ifdef CHECK_CXXLAPACK
90
91 template <typename MA, typename VTAU, typename VWORK>
92 void
93 lq2_native(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
94 {
95 typedef typename GeMatrix<MA>::ElementType T;
96
97 const INTEGER M = A.numRows();
98 const INTEGER N = A.numCols();
99 const INTEGER LDA = A.leadingDimension();
100 INTEGER INFO;
101
102 if (IsSame<T, DOUBLE>::value) {
103 LAPACK_IMPL(dgelq2)(&M, &N, A.data(), &LDA,
104 tau.data(), work.data(),
105 &INFO);
106 assert(INFO==0);
107 } else {
108 ASSERT(0);
109 }
110 }
111
112 #endif // CHECK_CXXLAPACK
113
114 //== public interface ==========================================================
115
116 template <typename MA, typename VTAU, typename VWORK>
117 void
118 lq2(GeMatrix<MA> &A, DenseVector<VTAU> &tau, DenseVector<VWORK> &work)
119 {
120 typedef typename GeMatrix<MA>::IndexType IndexType;
121
122 //
123 // Test the input parameters
124 //
125 # ifndef NDEBUG
126 ASSERT(A.firstRow()==1);
127 ASSERT(A.firstCol()==1);
128 ASSERT(tau.firstIndex()==1);
129 ASSERT(work.firstIndex()==1);
130
131 const IndexType m = A.numRows();
132 const IndexType n = A.numCols();
133 const IndexType k = std::min(m, n);
134
135 ASSERT(tau.length()>=k);
136 ASSERT(work.length()>=m);
137 # endif
138
139 # ifdef CHECK_CXXLAPACK
140 //
141 // Make copies of output arguments
142 //
143 typename GeMatrix<MA>::NoView A_org = A;
144 typename DenseVector<VTAU>::NoView tau_org = tau;
145 typename DenseVector<VTAU>::NoView work_org = work;
146 # endif
147
148 //
149 // Call implementation
150 //
151 lq2_generic(A, tau, work);
152
153 # ifdef CHECK_CXXLAPACK
154 //
155 // Restore output arguments
156 //
157 typename GeMatrix<MA>::NoView A_generic = A;
158 typename DenseVector<VTAU>::NoView tau_generic = tau;
159 typename DenseVector<VTAU>::NoView work_generic = work;
160
161 A = A_org;
162 tau = tau_org;
163 work = work_org;
164 //
165 // Compare results
166 //
167 lq2_native(A, tau, work);
168
169 bool failed = false;
170 if (! isIdentical(A_generic, A, "A_generic", "A")) {
171 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
172 std::cerr << "F77LAPACK: A = " << A << std::endl;
173 failed = true;
174 }
175
176 if (! isIdentical(tau_generic, tau, "tau_generic", "tau")) {
177 std::cerr << "CXXLAPACK: tau_generic = " << tau_generic << std::endl;
178 std::cerr << "F77LAPACK: tau = " << tau << std::endl;
179 failed = true;
180 }
181
182 if (! isIdentical(work_generic, work, "work_generic", "work")) {
183 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
184 std::cerr << "F77LAPACK: work = " << work << std::endl;
185 failed = true;
186 }
187
188 if (failed) {
189 ASSERT(0);
190 }
191 # endif
192 }
193
194 //-- forwarding ----------------------------------------------------------------
195 template <typename MA, typename VTAU, typename VWORK>
196 void
197 lq2(MA &&A, VTAU &&tau, VWORK &&work)
198 {
199 CHECKPOINT_ENTER;
200 lq2(A, tau, work);
201 CHECKPOINT_LEAVE;
202 }
203
204 } } // namespace lapack, flens
205
206 #endif // FLENS_LAPACK_LQ_LQ2_TCC