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 DLAQGE( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
36 $ EQUED )
37 *
38 * -- LAPACK auxiliary 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_LAQ_TCC
45 #define FLENS_LAPACK_EIG_LAQ_TCC 1
46
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53 template <typename MA, typename VR, typename VC,
54 typename ROWCOND, typename COLCOND,
55 typename AMAX>
56 LAQ::Equilibration
57 laq_generic(GeMatrix<MA> &A,
58 const DenseVector<VR> &r,
59 const DenseVector<VC> &c,
60 const ROWCOND &rowCond,
61 const COLCOND &colCond,
62 const AMAX &amax)
63 {
64 using namespace LAQ;
65
66 typedef typename GeMatrix<MA>::ElementType T;
67 typedef typename GeMatrix<MA>::IndexType IndexType;
68
69 const T One(1), Thresh(0.1);
70
71 const Underscore<IndexType> _;
72
73 const IndexType m = A.numRows();
74 const IndexType n = A.numCols();
75 //
76 // Quick return if possible
77 //
78 if (m==0 || n==0) {
79 return None;
80 }
81 //
82 // Initialize LARGE and SMALL.
83 //
84 const T small = lamch<T>(SafeMin) / lamch<T>(Precision);
85 const T large = One / small;
86
87 Equilibration equed;
88
89 if (rowCond>=Thresh && amax>=small && amax<=large) {
90 //
91 // No row scaling
92 //
93 if (colCond>=Thresh) {
94 //
95 // No column scaling
96 //
97 equed = None;
98 } else {
99 //
100 // Column scaling
101 //
102 for (IndexType j=1; j<=n; ++j) {
103 A(_,j) *= c(j);
104 }
105 equed = Column;
106 }
107 } else if (colCond>=Thresh) {
108 //
109 // Row scaling, no column scaling
110 //
111 for (IndexType j=1; j<=n; ++j) {
112 for (IndexType i=1; i<=m; ++i) {
113 A(i,j) *= r(i);
114 }
115 }
116 equed = Row;
117 } else {
118 //
119 // Row and column scaling
120 //
121 for (IndexType j=1; j<=n; ++j) {
122 const T cj = c(j);
123 for (IndexType i=1; i<=m; ++i) {
124 A(i,j) *= cj*r(i);
125 }
126 }
127 equed = Both;
128 }
129 return equed;
130 }
131
132
133 //== interface for native lapack ===============================================
134
135 #ifdef CHECK_CXXLAPACK
136
137 template <typename MA, typename VR, typename VC,
138 typename ROWCOND, typename COLCOND,
139 typename AMAX>
140 LAQ::Equilibration
141 laq_native(GeMatrix<MA> &A,
142 const DenseVector<VR> &r,
143 const DenseVector<VC> &c,
144 const ROWCOND &rowCond,
145 const COLCOND &colCond,
146 const AMAX &amax)
147 {
148 typedef typename GeMatrix<MA>::ElementType T;
149
150 const INTEGER M = A.numRows();
151 const INTEGER N = A.numCols();
152 const INTEGER LDA = A.leadingDimension();
153 char EQUED;
154
155 if (IsSame<T,double>::value) {
156 LAPACK_IMPL(dlaqge)(&M,
157 &N,
158 A.data(),
159 &LDA,
160 r.data(),
161 c.data(),
162 &rowCond,
163 &colCond,
164 &amax,
165 &EQUED);
166 } else {
167 ASSERT(0);
168 }
169
170 if (EQUED=='B') {
171 return LAQ::Both;
172 } else if (EQUED=='R') {
173 return LAQ::Row;
174 } else if (EQUED=='C') {
175 return LAQ::Column;
176 }
177
178 ASSERT(EQUED=='N');
179 return LAQ::None;
180 }
181
182 #endif // CHECK_CXXLAPACK
183
184 //== public interface ==========================================================
185 template <typename MA, typename VR, typename VC,
186 typename ROWCOND, typename COLCOND,
187 typename AMAX>
188 LAQ::Equilibration
189 laq(GeMatrix<MA> &A,
190 const DenseVector<VR> &r,
191 const DenseVector<VC> &c,
192 const ROWCOND &rowCond,
193 const COLCOND &colCond,
194 const AMAX &amax)
195 {
196 typedef typename GeMatrix<MA>::IndexType IndexType;
197 //
198 // Test the input parameters
199 //
200 # ifndef NDEBUG
201 ASSERT(A.firstRow()==1);
202 ASSERT(A.firstCol()==1);
203
204 const IndexType m = A.numRows();
205 const IndexType n = A.numCols();
206
207 ASSERT(r.firstIndex()==1);
208 ASSERT(r.length()==m);
209
210 ASSERT(c.firstIndex()==1);
211 ASSERT(c.length()==n);
212 # endif
213
214 # ifdef CHECK_CXXLAPACK
215 //
216 // Make copies of output arguments
217 //
218 typename GeMatrix<MA>::NoView A_org = A;
219 # endif
220
221 //
222 // Call implementation
223 //
224 const auto equed = laq_generic(A, r, c, rowCond, colCond, amax);
225
226 # ifdef CHECK_CXXLAPACK
227 //
228 // Make copies of results computed by the generic implementation
229 //
230 typename GeMatrix<MA>::NoView A_generic = A;
231 //
232 // restore output arguments
233 //
234 A = A_org;
235
236 //
237 // Compare generic results with results from the native implementation
238 //
239 const auto equed_ = laq_native(A, r, c, rowCond, colCond, amax);
240
241 bool failed = false;
242 if (! isIdentical(A_generic, A, "A_generic", "A")) {
243 std::cerr << "CXXLAPACK: A_generic = "
244 << A_generic << std::endl;
245 std::cerr << "F77LAPACK: A = " << A << std::endl;
246 failed = true;
247 }
248
249 if (equed!=equed_) {
250 std::cerr << "CXXLAPACK: equed = " << equed << std::endl;
251 std::cerr << "F77LAPACK: equed_ = " << equed_ << std::endl;
252 failed = true;
253 }
254
255 if (failed) {
256 ASSERT(0);
257 } else {
258 // std::cerr << "passed: laq.tcc" << std::endl;
259 }
260 # endif
261
262 return equed;
263 }
264
265
266 //-- forwarding ----------------------------------------------------------------
267 template <typename MA, typename VR, typename VC,
268 typename ROWCOND, typename COLCOND,
269 typename AMAX>
270 LAQ::Equilibration
271 laq(MA &&A,
272 const VR &r,
273 const VC &c,
274 const ROWCOND &rowCond,
275 const COLCOND &colCond,
276 const AMAX &amax)
277 {
278 CHECKPOINT_ENTER;
279 const auto equed = laq(A, r, c, rowCond, colCond, amax);
280 CHECKPOINT_LEAVE;
281
282 return equed;
283 }
284
285 } } // namespace lapack, flens
286
287 #endif // FLENS_LAPACK_EIG_LAQ_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 DLAQGE( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
36 $ EQUED )
37 *
38 * -- LAPACK auxiliary 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_LAQ_TCC
45 #define FLENS_LAPACK_EIG_LAQ_TCC 1
46
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53 template <typename MA, typename VR, typename VC,
54 typename ROWCOND, typename COLCOND,
55 typename AMAX>
56 LAQ::Equilibration
57 laq_generic(GeMatrix<MA> &A,
58 const DenseVector<VR> &r,
59 const DenseVector<VC> &c,
60 const ROWCOND &rowCond,
61 const COLCOND &colCond,
62 const AMAX &amax)
63 {
64 using namespace LAQ;
65
66 typedef typename GeMatrix<MA>::ElementType T;
67 typedef typename GeMatrix<MA>::IndexType IndexType;
68
69 const T One(1), Thresh(0.1);
70
71 const Underscore<IndexType> _;
72
73 const IndexType m = A.numRows();
74 const IndexType n = A.numCols();
75 //
76 // Quick return if possible
77 //
78 if (m==0 || n==0) {
79 return None;
80 }
81 //
82 // Initialize LARGE and SMALL.
83 //
84 const T small = lamch<T>(SafeMin) / lamch<T>(Precision);
85 const T large = One / small;
86
87 Equilibration equed;
88
89 if (rowCond>=Thresh && amax>=small && amax<=large) {
90 //
91 // No row scaling
92 //
93 if (colCond>=Thresh) {
94 //
95 // No column scaling
96 //
97 equed = None;
98 } else {
99 //
100 // Column scaling
101 //
102 for (IndexType j=1; j<=n; ++j) {
103 A(_,j) *= c(j);
104 }
105 equed = Column;
106 }
107 } else if (colCond>=Thresh) {
108 //
109 // Row scaling, no column scaling
110 //
111 for (IndexType j=1; j<=n; ++j) {
112 for (IndexType i=1; i<=m; ++i) {
113 A(i,j) *= r(i);
114 }
115 }
116 equed = Row;
117 } else {
118 //
119 // Row and column scaling
120 //
121 for (IndexType j=1; j<=n; ++j) {
122 const T cj = c(j);
123 for (IndexType i=1; i<=m; ++i) {
124 A(i,j) *= cj*r(i);
125 }
126 }
127 equed = Both;
128 }
129 return equed;
130 }
131
132
133 //== interface for native lapack ===============================================
134
135 #ifdef CHECK_CXXLAPACK
136
137 template <typename MA, typename VR, typename VC,
138 typename ROWCOND, typename COLCOND,
139 typename AMAX>
140 LAQ::Equilibration
141 laq_native(GeMatrix<MA> &A,
142 const DenseVector<VR> &r,
143 const DenseVector<VC> &c,
144 const ROWCOND &rowCond,
145 const COLCOND &colCond,
146 const AMAX &amax)
147 {
148 typedef typename GeMatrix<MA>::ElementType T;
149
150 const INTEGER M = A.numRows();
151 const INTEGER N = A.numCols();
152 const INTEGER LDA = A.leadingDimension();
153 char EQUED;
154
155 if (IsSame<T,double>::value) {
156 LAPACK_IMPL(dlaqge)(&M,
157 &N,
158 A.data(),
159 &LDA,
160 r.data(),
161 c.data(),
162 &rowCond,
163 &colCond,
164 &amax,
165 &EQUED);
166 } else {
167 ASSERT(0);
168 }
169
170 if (EQUED=='B') {
171 return LAQ::Both;
172 } else if (EQUED=='R') {
173 return LAQ::Row;
174 } else if (EQUED=='C') {
175 return LAQ::Column;
176 }
177
178 ASSERT(EQUED=='N');
179 return LAQ::None;
180 }
181
182 #endif // CHECK_CXXLAPACK
183
184 //== public interface ==========================================================
185 template <typename MA, typename VR, typename VC,
186 typename ROWCOND, typename COLCOND,
187 typename AMAX>
188 LAQ::Equilibration
189 laq(GeMatrix<MA> &A,
190 const DenseVector<VR> &r,
191 const DenseVector<VC> &c,
192 const ROWCOND &rowCond,
193 const COLCOND &colCond,
194 const AMAX &amax)
195 {
196 typedef typename GeMatrix<MA>::IndexType IndexType;
197 //
198 // Test the input parameters
199 //
200 # ifndef NDEBUG
201 ASSERT(A.firstRow()==1);
202 ASSERT(A.firstCol()==1);
203
204 const IndexType m = A.numRows();
205 const IndexType n = A.numCols();
206
207 ASSERT(r.firstIndex()==1);
208 ASSERT(r.length()==m);
209
210 ASSERT(c.firstIndex()==1);
211 ASSERT(c.length()==n);
212 # endif
213
214 # ifdef CHECK_CXXLAPACK
215 //
216 // Make copies of output arguments
217 //
218 typename GeMatrix<MA>::NoView A_org = A;
219 # endif
220
221 //
222 // Call implementation
223 //
224 const auto equed = laq_generic(A, r, c, rowCond, colCond, amax);
225
226 # ifdef CHECK_CXXLAPACK
227 //
228 // Make copies of results computed by the generic implementation
229 //
230 typename GeMatrix<MA>::NoView A_generic = A;
231 //
232 // restore output arguments
233 //
234 A = A_org;
235
236 //
237 // Compare generic results with results from the native implementation
238 //
239 const auto equed_ = laq_native(A, r, c, rowCond, colCond, amax);
240
241 bool failed = false;
242 if (! isIdentical(A_generic, A, "A_generic", "A")) {
243 std::cerr << "CXXLAPACK: A_generic = "
244 << A_generic << std::endl;
245 std::cerr << "F77LAPACK: A = " << A << std::endl;
246 failed = true;
247 }
248
249 if (equed!=equed_) {
250 std::cerr << "CXXLAPACK: equed = " << equed << std::endl;
251 std::cerr << "F77LAPACK: equed_ = " << equed_ << std::endl;
252 failed = true;
253 }
254
255 if (failed) {
256 ASSERT(0);
257 } else {
258 // std::cerr << "passed: laq.tcc" << std::endl;
259 }
260 # endif
261
262 return equed;
263 }
264
265
266 //-- forwarding ----------------------------------------------------------------
267 template <typename MA, typename VR, typename VC,
268 typename ROWCOND, typename COLCOND,
269 typename AMAX>
270 LAQ::Equilibration
271 laq(MA &&A,
272 const VR &r,
273 const VC &c,
274 const ROWCOND &rowCond,
275 const COLCOND &colCond,
276 const AMAX &amax)
277 {
278 CHECKPOINT_ENTER;
279 const auto equed = laq(A, r, c, rowCond, colCond, amax);
280 CHECKPOINT_LEAVE;
281
282 return equed;
283 }
284
285 } } // namespace lapack, flens
286
287 #endif // FLENS_LAPACK_EIG_LAQ_TCC