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 DORG2R( M, N, K, A, LDA, TAU, WORK, 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_QR_ORG2R_TCC
44 #define FLENS_LAPACK_QR_ORG2R_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 VWORK>
54 void
55 org2r_generic(IndexType k,
56 GeMatrix<MA> &A,
57 const DenseVector<VTAU> &tau,
58 DenseVector<VWORK> &work)
59 {
60 typedef typename GeMatrix<MA>::ElementType T;
61
62 const Underscore<IndexType> _;
63
64 const IndexType m = A.numRows();
65 const IndexType n = A.numCols();
66
67 const T Zero(0), One(1);
68 //
69 // Quick return if possible
70 //
71 if (n<=0) {
72 return;
73 }
74 //
75 // Initialise columns k+1:n to columns of the unit matrix
76 //
77 for (IndexType j=k+1; j<=n; ++j) {
78 A(_(1,m),j) = Zero;
79 A(j,j) = One;
80 }
81
82 for (IndexType i=k; i>=1; --i) {
83 //
84 // Apply H(i) to A(i:m,i+1:n) from the left
85 //
86 if (i<n) {
87 A(i,i) = One;
88 larf(Left, A(_(i,m),i), tau(i), A(_(i,m),_(i+1,n)), work);
89 }
90 if (i<m) {
91 blas::scal(-tau(i), A(_(i+1,m),i));
92 }
93 A(i,i) = One - tau(i);
94 //
95 // Set A(1:i-1,i) to zero
96 //
97 A(_(1,i-1),i) = Zero;
98 }
99 }
100
101 //== interface for native lapack ===============================================
102
103 #ifdef CHECK_CXXLAPACK
104
105 template <typename IndexType, typename MA, typename VTAU, typename VWORK>
106 void
107 org2r_native(IndexType k,
108 GeMatrix<MA> &A,
109 const DenseVector<VTAU> &tau,
110 DenseVector<VWORK> &work)
111 {
112 typedef typename GeMatrix<MA>::ElementType T;
113
114 const INTEGER M = A.numRows();
115 const INTEGER N = A.numCols();
116 const INTEGER K = k;
117 const INTEGER LDA = A.leadingDimension();
118 INTEGER INFO;
119
120 if (IsSame<T,DOUBLE>::value) {
121 LAPACK_IMPL(dorg2r)(&M,
122 &N,
123 &K,
124 A.data(),
125 &LDA,
126 tau.data(),
127 work.data(),
128 &INFO);
129 } else {
130 ASSERT(0);
131 }
132 ASSERT(INFO==0);
133 }
134
135 #endif // CHECK_CXXLAPACK
136
137 //== public interface ==========================================================
138
139 template <typename IndexType, typename MA, typename VTAU, typename VWORK>
140 void
141 org2r(IndexType k,
142 GeMatrix<MA> &A,
143 const DenseVector<VTAU> &tau,
144 DenseVector<VWORK> &work)
145 {
146 //
147 // Test the input parameters
148 //
149 # ifndef NDEBUG
150 ASSERT(A.firstRow()==IndexType(1));
151 ASSERT(A.firstCol()==IndexType(1));
152 ASSERT(tau.firstIndex()==IndexType(1));
153 ASSERT(work.length()>=A.numCols());
154
155 const IndexType m = A.numRows();
156 const IndexType n = A.numCols();
157
158 ASSERT(n<=m);
159 ASSERT(k<=n);
160 ASSERT(0<=k);
161 # endif
162
163 //
164 // Make copies of output arguments
165 //
166 # ifdef CHECK_CXXLAPACK
167 typename GeMatrix<MA>::NoView _A = A;
168 typename DenseVector<VWORK>::NoView _work = work;
169 # endif
170
171 //
172 // Call implementation
173 //
174 org2r_generic(k, A, tau, work);
175
176 //
177 // Compare results
178 //
179 # ifdef CHECK_CXXLAPACK
180 org2r_native(k, _A, tau, _work);
181
182 bool failed = false;
183 if (! isIdentical(A, _A, " A", "A_")) {
184 std::cerr << "CXXLAPACK: A = " << A << std::endl;
185 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
186 failed = true;
187 }
188
189 if (! isIdentical(work, _work, " work", "_work")) {
190 std::cerr << "CXXLAPACK: work = " << work << std::endl;
191 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
192 failed = true;
193 }
194
195 if (failed) {
196 std::cerr << "error in: org2r.tcc" << std::endl;
197 ASSERT(0);
198 } else {
199 // std::cerr << "passed: org2r.tcc" << std::endl;
200 }
201 # endif
202 }
203
204 //-- forwarding ----------------------------------------------------------------
205 template <typename IndexType, typename MA, typename VTAU, typename VWORK>
206 void
207 org2r(IndexType k, MA &&A, const VTAU &tau, VWORK &&work)
208 {
209 org2r(k, A, tau, work);
210 }
211
212 } } // namespace lapack, flens
213
214 #endif // FLENS_LAPACK_QR_ORG2R_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 DORG2R( M, N, K, A, LDA, TAU, WORK, 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_QR_ORG2R_TCC
44 #define FLENS_LAPACK_QR_ORG2R_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 VWORK>
54 void
55 org2r_generic(IndexType k,
56 GeMatrix<MA> &A,
57 const DenseVector<VTAU> &tau,
58 DenseVector<VWORK> &work)
59 {
60 typedef typename GeMatrix<MA>::ElementType T;
61
62 const Underscore<IndexType> _;
63
64 const IndexType m = A.numRows();
65 const IndexType n = A.numCols();
66
67 const T Zero(0), One(1);
68 //
69 // Quick return if possible
70 //
71 if (n<=0) {
72 return;
73 }
74 //
75 // Initialise columns k+1:n to columns of the unit matrix
76 //
77 for (IndexType j=k+1; j<=n; ++j) {
78 A(_(1,m),j) = Zero;
79 A(j,j) = One;
80 }
81
82 for (IndexType i=k; i>=1; --i) {
83 //
84 // Apply H(i) to A(i:m,i+1:n) from the left
85 //
86 if (i<n) {
87 A(i,i) = One;
88 larf(Left, A(_(i,m),i), tau(i), A(_(i,m),_(i+1,n)), work);
89 }
90 if (i<m) {
91 blas::scal(-tau(i), A(_(i+1,m),i));
92 }
93 A(i,i) = One - tau(i);
94 //
95 // Set A(1:i-1,i) to zero
96 //
97 A(_(1,i-1),i) = Zero;
98 }
99 }
100
101 //== interface for native lapack ===============================================
102
103 #ifdef CHECK_CXXLAPACK
104
105 template <typename IndexType, typename MA, typename VTAU, typename VWORK>
106 void
107 org2r_native(IndexType k,
108 GeMatrix<MA> &A,
109 const DenseVector<VTAU> &tau,
110 DenseVector<VWORK> &work)
111 {
112 typedef typename GeMatrix<MA>::ElementType T;
113
114 const INTEGER M = A.numRows();
115 const INTEGER N = A.numCols();
116 const INTEGER K = k;
117 const INTEGER LDA = A.leadingDimension();
118 INTEGER INFO;
119
120 if (IsSame<T,DOUBLE>::value) {
121 LAPACK_IMPL(dorg2r)(&M,
122 &N,
123 &K,
124 A.data(),
125 &LDA,
126 tau.data(),
127 work.data(),
128 &INFO);
129 } else {
130 ASSERT(0);
131 }
132 ASSERT(INFO==0);
133 }
134
135 #endif // CHECK_CXXLAPACK
136
137 //== public interface ==========================================================
138
139 template <typename IndexType, typename MA, typename VTAU, typename VWORK>
140 void
141 org2r(IndexType k,
142 GeMatrix<MA> &A,
143 const DenseVector<VTAU> &tau,
144 DenseVector<VWORK> &work)
145 {
146 //
147 // Test the input parameters
148 //
149 # ifndef NDEBUG
150 ASSERT(A.firstRow()==IndexType(1));
151 ASSERT(A.firstCol()==IndexType(1));
152 ASSERT(tau.firstIndex()==IndexType(1));
153 ASSERT(work.length()>=A.numCols());
154
155 const IndexType m = A.numRows();
156 const IndexType n = A.numCols();
157
158 ASSERT(n<=m);
159 ASSERT(k<=n);
160 ASSERT(0<=k);
161 # endif
162
163 //
164 // Make copies of output arguments
165 //
166 # ifdef CHECK_CXXLAPACK
167 typename GeMatrix<MA>::NoView _A = A;
168 typename DenseVector<VWORK>::NoView _work = work;
169 # endif
170
171 //
172 // Call implementation
173 //
174 org2r_generic(k, A, tau, work);
175
176 //
177 // Compare results
178 //
179 # ifdef CHECK_CXXLAPACK
180 org2r_native(k, _A, tau, _work);
181
182 bool failed = false;
183 if (! isIdentical(A, _A, " A", "A_")) {
184 std::cerr << "CXXLAPACK: A = " << A << std::endl;
185 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
186 failed = true;
187 }
188
189 if (! isIdentical(work, _work, " work", "_work")) {
190 std::cerr << "CXXLAPACK: work = " << work << std::endl;
191 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
192 failed = true;
193 }
194
195 if (failed) {
196 std::cerr << "error in: org2r.tcc" << std::endl;
197 ASSERT(0);
198 } else {
199 // std::cerr << "passed: org2r.tcc" << std::endl;
200 }
201 # endif
202 }
203
204 //-- forwarding ----------------------------------------------------------------
205 template <typename IndexType, typename MA, typename VTAU, typename VWORK>
206 void
207 org2r(IndexType k, MA &&A, const VTAU &tau, VWORK &&work)
208 {
209 org2r(k, A, tau, work);
210 }
211
212 } } // namespace lapack, flens
213
214 #endif // FLENS_LAPACK_QR_ORG2R_TCC