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 DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
36 IMPLICIT NONE
37 *
38 * -- LAPACK auxiliary routine (version 3.3.1) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * -- April 2011 --
42 */
43
44 #ifndef FLENS_LAPACK_AUX_LARFT_TCC
45 #define FLENS_LAPACK_AUX_LARFT_TCC 1
46
47 #include <flens/blas/blas.h>
48
49 namespace flens { namespace lapack {
50
51 //== generic lapack implementation =============================================
52
53 template <typename IndexType, typename MV, typename VTAU, typename MT>
54 void
55 larft_generic(Direction direction, StoreVectors storeVectors, IndexType n,
56 GeMatrix<MV> &V, const DenseVector<VTAU> &tau, TrMatrix<MT> &Tr)
57 {
58 using std::max;
59 using std::min;
60
61 const Underscore<IndexType> _;
62
63 typedef typename GeMatrix<MV>::ElementType T;
64
65 const T Zero(0);
66
67 //
68 // Quick return if possible
69 //
70 if (n==0) {
71 return;
72 }
73
74 const IndexType k = Tr.dim();
75
76 // Lehn: as long as we do not have col-views for TrMatrix we get them
77 // via a GeMatrixView
78 auto _Tr = Tr.general();
79
80 if (direction==Forward) {
81 IndexType lastV = -1;
82 IndexType prevLastV = n;
83 for (IndexType i=1; i<=k; ++i) {
84 prevLastV = max(i, prevLastV);
85 if (tau(i)==Zero) {
86 //
87 // H(i) = I
88 //
89 for (IndexType j=1; j<=i; ++j) {
90 Tr(j,i) = Zero;
91 }
92 } else {
93 //
94 // general case
95 //
96 T Vii = V(i,i);
97 V(i,i) = 1;
98 if (storeVectors==ColumnWise) {
99 // Skip any trailing zeros.
100 for (lastV=n; lastV>=i+1; --lastV) {
101 if (V(lastV,i)!=Zero) {
102 break;
103 }
104 }
105 IndexType j = min(lastV, prevLastV);
106 //
107 // T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
108 //
109 blas::mv(Trans, -tau(i),
110 V(_(i,j),_(1,i-1)), V(_(i,j),i),
111 Zero,
112 _Tr(_(1,i-1),i));
113 } else { /* storeVectors==RowWise */
114 // Skip any trailing zeros.
115 for (lastV=n; lastV>=i+1; --lastV) {
116 if (V(i,lastV)!=Zero) {
117 break;
118 }
119 }
120 IndexType j = min(lastV, prevLastV);
121 //
122 // T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
123 //
124 blas::mv(NoTrans, -tau(i),
125 V(_(1,i-1),_(i,j)), V(i,_(i,j)),
126 Zero,
127 _Tr(_(1,i-1),i));
128 }
129 V(i,i) = Vii;
130 //
131 // T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
132 //
133 blas::mv(NoTrans,
134 _Tr(_(1,i-1),_(1,i-1)).upper(),
135 _Tr(_(1,i-1),i));
136 Tr(i,i) = tau(i);
137 if (i>1) {
138 prevLastV = max(prevLastV, lastV);
139 } else {
140 prevLastV = lastV;
141 }
142 }
143 }
144 } else {
145 // Lehn: I will implement it as soon as someone needs this case
146 ASSERT(0);
147 }
148 }
149
150 //== interface for native lapack ===============================================
151
152 #ifdef CHECK_CXXLAPACK
153
154 template <typename IndexType, typename MV, typename VTAU, typename MT>
155 void
156 larft_native(Direction direction, StoreVectors storeVectors, IndexType n,
157 GeMatrix<MV> &V, const DenseVector<VTAU> &tau, TrMatrix<MT> &Tr)
158 {
159 typedef typename TrMatrix<MT>::ElementType T;
160
161 const char DIRECT = (direction==Forward) ? 'F' : 'B';
162 const char STOREV = (storeVectors==ColumnWise) ? 'C' : 'R';
163 const INTEGER N = n;
164 const INTEGER K = Tr.dim();
165 const INTEGER LDV = V.leadingDimension();
166 const INTEGER LDT = Tr.leadingDimension();
167
168 if (IsSame<T, DOUBLE>::value) {
169 LAPACK_IMPL(dlarft)(&DIRECT,
170 &STOREV,
171 &N,
172 &K,
173 V.data(),
174 &LDV,
175 tau.data(),
176 Tr.data(),
177 &LDT);
178 } else {
179 ASSERT(0);
180 }
181 }
182
183 #endif // CHECK_CXXLAPACK
184
185 //== public interface ==========================================================
186
187 template <typename IndexType, typename MV, typename VTAU, typename MT>
188 void
189 larft(Direction direction, StoreVectors storeVectors, IndexType n,
190 GeMatrix<MV> &V, const DenseVector<VTAU> &tau, TrMatrix<MT> &Tr)
191 {
192 LAPACK_DEBUG_OUT("larft");
193
194 //
195 // Test the input parameters
196 //
197 ASSERT(V.firstRow()==1);
198 ASSERT(V.firstCol()==1);
199 ASSERT(tau.firstIndex()==1);
200 ASSERT(Tr.firstRow()==1);
201 ASSERT(Tr.firstCol()==1);
202
203 # ifdef CHECK_CXXLAPACK
204 //
205 // Make copies of output arguments
206 //
207 typename GeMatrix<MV>::NoView V_org = V;
208 // copy the full storage!
209 typename GeMatrix<MT>::NoView Tr_org = Tr.general();
210 # endif
211
212 //
213 // Call implementation
214 //
215 larft_generic(direction, storeVectors, n, V, tau, Tr);
216
217 # ifdef CHECK_CXXLAPACK
218 //
219 // Restore output arguments
220 //
221 typename GeMatrix<MV>::NoView V_generic = V;
222 typename GeMatrix<MT>::NoView Tr_generic = Tr.general();
223
224 V = V_org;
225 Tr.general() = Tr_org;
226 //
227 // Compare results
228 //
229 larft_native(direction, storeVectors, n, V, tau, Tr);
230
231 bool failed = false;
232 if (! isIdentical(V_generic, V, "V_generic", "V")) {
233 std::cerr << "CXXLAPACK: V_generic = " << V_generic << std::endl;
234 std::cerr << "F77LAPACK: V = " << V << std::endl;
235 failed = true;
236 }
237
238 if (! isIdentical(Tr_generic, Tr.general(), "Tr_generic", "_Tr"))
239 {
240 std::cerr << "CXXLAPACK: Tr_generic = " << Tr_generic << std::endl;
241 std::cerr << "F77LAPACK: Tr = " << Tr.general() << std::endl;
242 failed = true;
243 }
244
245 if (failed) {
246 ASSERT(0);
247 }
248 # endif
249 }
250
251 //-- forwarding ----------------------------------------------------------------
252 template <typename IndexType, typename MV, typename VTAU, typename MT>
253 void
254 larft(Direction direction, StoreVectors storeVectors, IndexType n,
255 MV &&V, const VTAU &tau, MT &&Tr)
256 {
257 larft(direction, storeVectors, n, V, tau, Tr);
258 }
259
260 } } // namespace lapack, flens
261
262 #endif // FLENS_LAPACK_AUX_LARFT_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 DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
36 IMPLICIT NONE
37 *
38 * -- LAPACK auxiliary routine (version 3.3.1) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * -- April 2011 --
42 */
43
44 #ifndef FLENS_LAPACK_AUX_LARFT_TCC
45 #define FLENS_LAPACK_AUX_LARFT_TCC 1
46
47 #include <flens/blas/blas.h>
48
49 namespace flens { namespace lapack {
50
51 //== generic lapack implementation =============================================
52
53 template <typename IndexType, typename MV, typename VTAU, typename MT>
54 void
55 larft_generic(Direction direction, StoreVectors storeVectors, IndexType n,
56 GeMatrix<MV> &V, const DenseVector<VTAU> &tau, TrMatrix<MT> &Tr)
57 {
58 using std::max;
59 using std::min;
60
61 const Underscore<IndexType> _;
62
63 typedef typename GeMatrix<MV>::ElementType T;
64
65 const T Zero(0);
66
67 //
68 // Quick return if possible
69 //
70 if (n==0) {
71 return;
72 }
73
74 const IndexType k = Tr.dim();
75
76 // Lehn: as long as we do not have col-views for TrMatrix we get them
77 // via a GeMatrixView
78 auto _Tr = Tr.general();
79
80 if (direction==Forward) {
81 IndexType lastV = -1;
82 IndexType prevLastV = n;
83 for (IndexType i=1; i<=k; ++i) {
84 prevLastV = max(i, prevLastV);
85 if (tau(i)==Zero) {
86 //
87 // H(i) = I
88 //
89 for (IndexType j=1; j<=i; ++j) {
90 Tr(j,i) = Zero;
91 }
92 } else {
93 //
94 // general case
95 //
96 T Vii = V(i,i);
97 V(i,i) = 1;
98 if (storeVectors==ColumnWise) {
99 // Skip any trailing zeros.
100 for (lastV=n; lastV>=i+1; --lastV) {
101 if (V(lastV,i)!=Zero) {
102 break;
103 }
104 }
105 IndexType j = min(lastV, prevLastV);
106 //
107 // T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
108 //
109 blas::mv(Trans, -tau(i),
110 V(_(i,j),_(1,i-1)), V(_(i,j),i),
111 Zero,
112 _Tr(_(1,i-1),i));
113 } else { /* storeVectors==RowWise */
114 // Skip any trailing zeros.
115 for (lastV=n; lastV>=i+1; --lastV) {
116 if (V(i,lastV)!=Zero) {
117 break;
118 }
119 }
120 IndexType j = min(lastV, prevLastV);
121 //
122 // T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
123 //
124 blas::mv(NoTrans, -tau(i),
125 V(_(1,i-1),_(i,j)), V(i,_(i,j)),
126 Zero,
127 _Tr(_(1,i-1),i));
128 }
129 V(i,i) = Vii;
130 //
131 // T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
132 //
133 blas::mv(NoTrans,
134 _Tr(_(1,i-1),_(1,i-1)).upper(),
135 _Tr(_(1,i-1),i));
136 Tr(i,i) = tau(i);
137 if (i>1) {
138 prevLastV = max(prevLastV, lastV);
139 } else {
140 prevLastV = lastV;
141 }
142 }
143 }
144 } else {
145 // Lehn: I will implement it as soon as someone needs this case
146 ASSERT(0);
147 }
148 }
149
150 //== interface for native lapack ===============================================
151
152 #ifdef CHECK_CXXLAPACK
153
154 template <typename IndexType, typename MV, typename VTAU, typename MT>
155 void
156 larft_native(Direction direction, StoreVectors storeVectors, IndexType n,
157 GeMatrix<MV> &V, const DenseVector<VTAU> &tau, TrMatrix<MT> &Tr)
158 {
159 typedef typename TrMatrix<MT>::ElementType T;
160
161 const char DIRECT = (direction==Forward) ? 'F' : 'B';
162 const char STOREV = (storeVectors==ColumnWise) ? 'C' : 'R';
163 const INTEGER N = n;
164 const INTEGER K = Tr.dim();
165 const INTEGER LDV = V.leadingDimension();
166 const INTEGER LDT = Tr.leadingDimension();
167
168 if (IsSame<T, DOUBLE>::value) {
169 LAPACK_IMPL(dlarft)(&DIRECT,
170 &STOREV,
171 &N,
172 &K,
173 V.data(),
174 &LDV,
175 tau.data(),
176 Tr.data(),
177 &LDT);
178 } else {
179 ASSERT(0);
180 }
181 }
182
183 #endif // CHECK_CXXLAPACK
184
185 //== public interface ==========================================================
186
187 template <typename IndexType, typename MV, typename VTAU, typename MT>
188 void
189 larft(Direction direction, StoreVectors storeVectors, IndexType n,
190 GeMatrix<MV> &V, const DenseVector<VTAU> &tau, TrMatrix<MT> &Tr)
191 {
192 LAPACK_DEBUG_OUT("larft");
193
194 //
195 // Test the input parameters
196 //
197 ASSERT(V.firstRow()==1);
198 ASSERT(V.firstCol()==1);
199 ASSERT(tau.firstIndex()==1);
200 ASSERT(Tr.firstRow()==1);
201 ASSERT(Tr.firstCol()==1);
202
203 # ifdef CHECK_CXXLAPACK
204 //
205 // Make copies of output arguments
206 //
207 typename GeMatrix<MV>::NoView V_org = V;
208 // copy the full storage!
209 typename GeMatrix<MT>::NoView Tr_org = Tr.general();
210 # endif
211
212 //
213 // Call implementation
214 //
215 larft_generic(direction, storeVectors, n, V, tau, Tr);
216
217 # ifdef CHECK_CXXLAPACK
218 //
219 // Restore output arguments
220 //
221 typename GeMatrix<MV>::NoView V_generic = V;
222 typename GeMatrix<MT>::NoView Tr_generic = Tr.general();
223
224 V = V_org;
225 Tr.general() = Tr_org;
226 //
227 // Compare results
228 //
229 larft_native(direction, storeVectors, n, V, tau, Tr);
230
231 bool failed = false;
232 if (! isIdentical(V_generic, V, "V_generic", "V")) {
233 std::cerr << "CXXLAPACK: V_generic = " << V_generic << std::endl;
234 std::cerr << "F77LAPACK: V = " << V << std::endl;
235 failed = true;
236 }
237
238 if (! isIdentical(Tr_generic, Tr.general(), "Tr_generic", "_Tr"))
239 {
240 std::cerr << "CXXLAPACK: Tr_generic = " << Tr_generic << std::endl;
241 std::cerr << "F77LAPACK: Tr = " << Tr.general() << std::endl;
242 failed = true;
243 }
244
245 if (failed) {
246 ASSERT(0);
247 }
248 # endif
249 }
250
251 //-- forwarding ----------------------------------------------------------------
252 template <typename IndexType, typename MV, typename VTAU, typename MT>
253 void
254 larft(Direction direction, StoreVectors storeVectors, IndexType n,
255 MV &&V, const VTAU &tau, MT &&Tr)
256 {
257 larft(direction, storeVectors, n, V, tau, Tr);
258 }
259
260 } } // namespace lapack, flens
261
262 #endif // FLENS_LAPACK_AUX_LARFT_TCC