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 DGEBAK( JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV,
36 $ 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_BAK_TCC
45 #define FLENS_LAPACK_EIG_BAK_TCC 1
46
47 #include <flens/aux/aux.h>
48 #include <flens/blas/blas.h>
49 #include <flens/lapack/lapack.h>
50
51 namespace flens { namespace lapack {
52
53 //== generic lapack implementation =============================================
54
55 template <typename IndexType, typename VSCALE, typename MV>
56 void
57 bak_generic(BALANCE::Balance job,
58 Side side,
59 IndexType iLo,
60 IndexType iHi,
61 const DenseVector<VSCALE> &scale,
62 GeMatrix<MV> &V)
63 {
64 using namespace BALANCE;
65
66 typedef typename GeMatrix<MV>::ElementType T;
67 const T One(1);
68
69 const Underscore<IndexType> _;
70
71 const IndexType n = V.numRows();
72 const IndexType m = V.numCols();
73 //
74 // Quick return if possible
75 //
76 if (n==0) {
77 return;
78 }
79 if (m==0) {
80 return;
81 }
82 if (job==None) {
83 return;
84 }
85
86 if (iLo!=iHi) {
87 //
88 // Backward balance
89 //
90 if (job==ScaleOnly || job==Both) {
91
92 if (side==Right) {
93 for (IndexType i=iLo; i<=iHi; ++i) {
94 V(i,_) *= scale(i);
95 }
96 }
97
98 if (side==Left) {
99 for (IndexType i=iLo; i<=iHi; ++i) {
100 V(i,_) *= One / scale(i);
101 }
102 }
103
104 }
105 }
106 //
107 // Backward permutation
108 //
109 // For I = ILO-1 step -1 until 1,
110 // IHI+1 step 1 until N do --
111 //
112 if (job==PermuteOnly || job==Both) {
113 if (side==Right) {
114 for (IndexType ii=1; ii<=n; ++ii) {
115 IndexType i = ii;
116 if (i>=iLo && i<=iHi) {
117 continue;
118 }
119 if (i<iLo) {
120 i = iLo - ii;
121 }
122 const IndexType k = explicit_cast<T,IndexType>(scale(i));
123 if (k==i) {
124 continue;
125 }
126 blas::swap(V(i,_), V(k,_));
127 }
128 }
129
130 if (side==Left) {
131 for (IndexType ii=1; ii<=n; ++ii) {
132 IndexType i = ii;
133 if (i>=iLo && i<=iHi) {
134 continue;
135 }
136 if (i<iLo) {
137 i = iLo - ii;
138 }
139 const IndexType k = explicit_cast<T,IndexType>(scale(i));
140 if (k==i) {
141 continue;
142 }
143 blas::swap(V(i,_), V(k,_));
144 }
145 }
146 }
147 }
148
149 //== interface for native lapack ===============================================
150
151 #ifdef CHECK_CXXLAPACK
152
153 template <typename IndexType, typename VSCALE, typename MV>
154 void
155 bak_native(BALANCE::Balance job,
156 Side side,
157 IndexType iLo,
158 IndexType iHi,
159 const DenseVector<VSCALE> &scale,
160 GeMatrix<MV> &V)
161 {
162 typedef typename GeMatrix<MV>::ElementType T;
163
164 const char JOB = getF77LapackChar<BALANCE::Balance>(job);
165 const char SIDE = getF77LapackChar<Side>(side);
166 const INTEGER N = V.numRows();
167 const INTEGER ILO = iLo;
168 const INTEGER IHI = iHi;
169 const INTEGER M = V.numCols();
170 const INTEGER LDV = V.leadingDimension();
171 INTEGER INFO;
172
173 if (IsSame<T, DOUBLE>::value) {
174 LAPACK_IMPL(dgebak)(&JOB,
175 &SIDE,
176 &N,
177 &ILO,
178 &IHI,
179 scale.data(),
180 &M,
181 V.data(),
182 &LDV,
183 &INFO);
184 } else {
185 ASSERT(0);
186 }
187 ASSERT(INFO==0);
188 }
189
190 #endif
191
192 //== public interface ==========================================================
193
194 template <typename IndexType, typename VSCALE, typename MV>
195 void
196 bak(BALANCE::Balance job,
197 Side side,
198 IndexType iLo,
199 IndexType iHi,
200 const DenseVector<VSCALE> &scale,
201 GeMatrix<MV> &V)
202 {
203 LAPACK_DEBUG_OUT("bak");
204
205 # ifndef NDEBUG
206 //
207 // Test the input parameters
208 //
209 const IndexType n = V.numRows();
210
211 if (n>0) {
212 ASSERT(1<=iLo);
213 ASSERT(iLo<=iHi);
214 ASSERT(iHi<=n);
215 } else {
216 ASSERT(iLo==1);
217 ASSERT(iHi==0);
218 }
219 # endif
220
221 # ifdef CHECK_CXXLAPACK
222 //
223 // Make copies of output arguments
224 //
225 typename GeMatrix<MV>::NoView _V = V;
226 # endif
227
228 //
229 // Call implementation
230 //
231 bak_generic(job, side, iLo, iHi, scale, V);
232
233 # ifdef CHECK_CXXLAPACK
234 //
235 // Compare results
236 //
237 bak_native(job, side, iLo, iHi, scale, _V);
238
239 if (! isIdentical(V, _V, " V", "_V")) {
240 std::cerr << "CXXLAPACK: V = " << V << std::endl;
241 std::cerr << "F77LAPACK: _V = " << _V << std::endl;
242 ASSERT(0);
243 }
244 # endif
245 }
246
247 //-- forwarding ----------------------------------------------------------------
248 template <typename IndexType, typename SCALE, typename MV>
249 void
250 bak(BALANCE::Balance job,
251 Side side,
252 IndexType iLo,
253 IndexType iHi,
254 const SCALE &scale,
255 MV &&V)
256 {
257 CHECKPOINT_ENTER;
258 bak(job, side, iLo, iHi, scale, V);
259 CHECKPOINT_LEAVE;
260 }
261
262 } } // namespace lapack, flens
263
264 #endif // FLENS_LAPACK_EIG_BAK_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 DGEBAK( JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV,
36 $ 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_BAK_TCC
45 #define FLENS_LAPACK_EIG_BAK_TCC 1
46
47 #include <flens/aux/aux.h>
48 #include <flens/blas/blas.h>
49 #include <flens/lapack/lapack.h>
50
51 namespace flens { namespace lapack {
52
53 //== generic lapack implementation =============================================
54
55 template <typename IndexType, typename VSCALE, typename MV>
56 void
57 bak_generic(BALANCE::Balance job,
58 Side side,
59 IndexType iLo,
60 IndexType iHi,
61 const DenseVector<VSCALE> &scale,
62 GeMatrix<MV> &V)
63 {
64 using namespace BALANCE;
65
66 typedef typename GeMatrix<MV>::ElementType T;
67 const T One(1);
68
69 const Underscore<IndexType> _;
70
71 const IndexType n = V.numRows();
72 const IndexType m = V.numCols();
73 //
74 // Quick return if possible
75 //
76 if (n==0) {
77 return;
78 }
79 if (m==0) {
80 return;
81 }
82 if (job==None) {
83 return;
84 }
85
86 if (iLo!=iHi) {
87 //
88 // Backward balance
89 //
90 if (job==ScaleOnly || job==Both) {
91
92 if (side==Right) {
93 for (IndexType i=iLo; i<=iHi; ++i) {
94 V(i,_) *= scale(i);
95 }
96 }
97
98 if (side==Left) {
99 for (IndexType i=iLo; i<=iHi; ++i) {
100 V(i,_) *= One / scale(i);
101 }
102 }
103
104 }
105 }
106 //
107 // Backward permutation
108 //
109 // For I = ILO-1 step -1 until 1,
110 // IHI+1 step 1 until N do --
111 //
112 if (job==PermuteOnly || job==Both) {
113 if (side==Right) {
114 for (IndexType ii=1; ii<=n; ++ii) {
115 IndexType i = ii;
116 if (i>=iLo && i<=iHi) {
117 continue;
118 }
119 if (i<iLo) {
120 i = iLo - ii;
121 }
122 const IndexType k = explicit_cast<T,IndexType>(scale(i));
123 if (k==i) {
124 continue;
125 }
126 blas::swap(V(i,_), V(k,_));
127 }
128 }
129
130 if (side==Left) {
131 for (IndexType ii=1; ii<=n; ++ii) {
132 IndexType i = ii;
133 if (i>=iLo && i<=iHi) {
134 continue;
135 }
136 if (i<iLo) {
137 i = iLo - ii;
138 }
139 const IndexType k = explicit_cast<T,IndexType>(scale(i));
140 if (k==i) {
141 continue;
142 }
143 blas::swap(V(i,_), V(k,_));
144 }
145 }
146 }
147 }
148
149 //== interface for native lapack ===============================================
150
151 #ifdef CHECK_CXXLAPACK
152
153 template <typename IndexType, typename VSCALE, typename MV>
154 void
155 bak_native(BALANCE::Balance job,
156 Side side,
157 IndexType iLo,
158 IndexType iHi,
159 const DenseVector<VSCALE> &scale,
160 GeMatrix<MV> &V)
161 {
162 typedef typename GeMatrix<MV>::ElementType T;
163
164 const char JOB = getF77LapackChar<BALANCE::Balance>(job);
165 const char SIDE = getF77LapackChar<Side>(side);
166 const INTEGER N = V.numRows();
167 const INTEGER ILO = iLo;
168 const INTEGER IHI = iHi;
169 const INTEGER M = V.numCols();
170 const INTEGER LDV = V.leadingDimension();
171 INTEGER INFO;
172
173 if (IsSame<T, DOUBLE>::value) {
174 LAPACK_IMPL(dgebak)(&JOB,
175 &SIDE,
176 &N,
177 &ILO,
178 &IHI,
179 scale.data(),
180 &M,
181 V.data(),
182 &LDV,
183 &INFO);
184 } else {
185 ASSERT(0);
186 }
187 ASSERT(INFO==0);
188 }
189
190 #endif
191
192 //== public interface ==========================================================
193
194 template <typename IndexType, typename VSCALE, typename MV>
195 void
196 bak(BALANCE::Balance job,
197 Side side,
198 IndexType iLo,
199 IndexType iHi,
200 const DenseVector<VSCALE> &scale,
201 GeMatrix<MV> &V)
202 {
203 LAPACK_DEBUG_OUT("bak");
204
205 # ifndef NDEBUG
206 //
207 // Test the input parameters
208 //
209 const IndexType n = V.numRows();
210
211 if (n>0) {
212 ASSERT(1<=iLo);
213 ASSERT(iLo<=iHi);
214 ASSERT(iHi<=n);
215 } else {
216 ASSERT(iLo==1);
217 ASSERT(iHi==0);
218 }
219 # endif
220
221 # ifdef CHECK_CXXLAPACK
222 //
223 // Make copies of output arguments
224 //
225 typename GeMatrix<MV>::NoView _V = V;
226 # endif
227
228 //
229 // Call implementation
230 //
231 bak_generic(job, side, iLo, iHi, scale, V);
232
233 # ifdef CHECK_CXXLAPACK
234 //
235 // Compare results
236 //
237 bak_native(job, side, iLo, iHi, scale, _V);
238
239 if (! isIdentical(V, _V, " V", "_V")) {
240 std::cerr << "CXXLAPACK: V = " << V << std::endl;
241 std::cerr << "F77LAPACK: _V = " << _V << std::endl;
242 ASSERT(0);
243 }
244 # endif
245 }
246
247 //-- forwarding ----------------------------------------------------------------
248 template <typename IndexType, typename SCALE, typename MV>
249 void
250 bak(BALANCE::Balance job,
251 Side side,
252 IndexType iLo,
253 IndexType iHi,
254 const SCALE &scale,
255 MV &&V)
256 {
257 CHECKPOINT_ENTER;
258 bak(job, side, iLo, iHi, scale, V);
259 CHECKPOINT_LEAVE;
260 }
261
262 } } // namespace lapack, flens
263
264 #endif // FLENS_LAPACK_EIG_BAK_TCC