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 DLASCL( TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO )
36 *
37 * -- LAPACK auxiliary routine (version 3.3.0) --
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 2010
41 */
42
43 #ifndef FLENS_LAPACK_AUX_LASCL_TCC
44 #define FLENS_LAPACK_AUX_LASCL_TCC 1
45
46 #include <cmath>
47 #include <flens/matrixtypes/matrixtypes.h>
48 #include <flens/vectortypes/vectortypes.h>
49
50 ////////
51 ////////
52 ////////
53 ////////
54 ////////
55 ////////
56 ////////
57 ////////
58 //////// TODO: Completely redesign this crap!
59 ////////
60 ////////
61 ////////
62 //////// ... but it works.
63 ////////
64 ////////
65 ////////
66 ////////
67 ////////
68
69 namespace flens { namespace lapack {
70
71 //== generic lapack implementation =============================================
72
73 template <typename IndexType, typename T, typename MA>
74 typename RestrictTo<IsSame<typename MA::ElementType, T>::value, void>::Type
75 lascl_generic(LASCL::Type type,
76 IndexType kl,
77 IndexType ku,
78 const T &cFrom,
79 const T &cTo,
80 MA &A)
81 {
82 using namespace LASCL;
83
84 using std::abs;
85 using std::isnan;
86 using std::max;
87 using std::min;
88
89 const IndexType m = A.numRows();
90 const IndexType n = A.numCols();
91 const T Zero = 0;
92 const T One = 1;
93
94 if ((cFrom==Zero) || isnan(cFrom)) {
95 ASSERT(0);
96 } else if (isnan(cTo)) {
97 ASSERT(0);
98 } else if (type==SymmetricLowerBand) {
99 ASSERT(A.numRows()==A.numCols());
100 } else if (type==SymmetricUpperBand) {
101 ASSERT(A.numRows()==A.numCols());
102 }
103 //
104 // Quick return if possible
105 //
106 if ((n==0) || (m==0)) {
107 return;
108 }
109 //
110 // Get machine parameters
111 //
112 const T smallNum = lamch<T>(SafeMin);
113 const T bigNum = T(1) / smallNum;
114 //
115 // Make copies of cFrom, cTo
116 //
117 T cFromC = cFrom;
118 T cToC = cTo;
119
120 bool done = false;
121
122 do {
123 T cFrom1 = cFromC*smallNum;
124 T cTo1;
125 T mul;
126 if (cFrom1==cFromC) {
127 // cFromC is an inf. Multiply by a correctly signed zero for
128 // finite cToC, or a NaN if cToC is infinite.
129 mul = cToC / cFromC;
130 done = true;
131 cTo1 = cToC;
132 } else {
133 cTo1 = cToC / bigNum;
134 if (cTo1==cToC) {
135 // cToC is either 0 or an inf. In both cases, cToC itself
136 // serves as the correct multiplication factor.
137 mul = cToC;
138 done = true;
139 cFromC = One;
140 } else if (abs(cFrom1)>abs(cToC) && cToC!=Zero) {
141 mul = smallNum;
142 done = false;
143 cFromC = cFrom1;
144 } else if (abs(cTo1)>abs(cFromC)) {
145 mul = bigNum;
146 done = false;
147 cToC = cTo1;
148 } else {
149 mul = cToC / cFromC;
150 done = true;
151 }
152 }
153
154 if (type==FullMatrix) {
155 //
156 // Full matrix
157 //
158 for (IndexType j=1; j<=n; ++j) {
159 for (IndexType i=1; i<=m; ++i) {
160 A(i,j) *= mul;
161 }
162 }
163 } else if (type==LowerTriangular) {
164 //
165 // Lower triangular matrix
166 //
167 for (IndexType j=1; j<=n; ++j) {
168 for (IndexType i=j; i<=m; ++i) {
169 A(i,j) *= mul;
170 }
171 }
172 } else if (type==UpperTriangular) {
173 //
174 // Upper triangular matrix
175 //
176 for (IndexType j=1; j<=n; ++j) {
177 for (IndexType i=1; i<=min(j,m); ++i) {
178 A(i,j) *= mul;
179 }
180 }
181 } else if (type==UpperHessenberg) {
182 //
183 // Upper Hessenberg matrix
184 //
185 for (IndexType j=1; j<=n; ++j) {
186 for (IndexType i=1; i<=min(j+1,m); ++i) {
187 A(i,j) *= mul;
188 }
189 }
190 } else if (type==SymmetricLowerBand) {
191 //
192 // Lower half of a symmetric band matrix
193 //
194 // TODO: this only works for the internal fullstorage of a
195 // band matrix not for the external element access
196 // of SbMatrix
197 const IndexType k3 = kl + 1;
198 const IndexType k4 = n + 1;
199 for (IndexType j=1; j<=n; ++j) {
200 for (IndexType i=1; i<=min(k3, k4-j); ++i) {
201 A(i,j) *= mul;
202 }
203 }
204 } else if (type==SymmetricUpperBand) {
205 //
206 // Upper half of a symmetric band matrix
207 //
208 // TODO: this only works for the internal fullstorage of a
209 // band matrix not for the external element access
210 // of SbMatrix
211 const IndexType k1 = ku + 2;
212 const IndexType k3 = ku + 1;
213 for (IndexType j=1; j<=n; ++j) {
214 for (IndexType i=max(k1-j,IndexType(1)); i<=k3; ++i) {
215 A(i,j) *= mul;
216 }
217 }
218 } else if (type==GeneralBand) {
219 //
220 // Band matrix
221 //
222 // TODO: this only works for the internal fullstorage of a
223 // band matrix not for the external element access
224 // of GeMatrix
225 const IndexType k1 = kl + ku + 2;
226 const IndexType k2 = kl + 1;
227 const IndexType k3 = 2*kl + ku + 1;
228 const IndexType k4 = kl + ku + 1 + m;
229 for (IndexType j=1; j<=n; ++j) {
230 for (IndexType i=max(k1-j,k2); i<=min(k3,k4-j); ++i) {
231 A(i,j) *= mul;
232 }
233 }
234 } else {
235 ASSERT(0);
236 }
237 } while (!done);
238 }
239
240 //== interface for native lapack ===============================================
241
242 #ifdef CHECK_CXXLAPACK
243
244 template <typename IndexType, typename T, typename MA>
245 void
246 lascl_native(LASCL::Type type,
247 IndexType kl,
248 IndexType ku,
249 const T cFrom,
250 const T &cTo,
251 MA &A)
252 {
253 const char TYPE = getF77LapackChar<LASCL::Type>(type);
254 const INTEGER KL = kl;
255 const INTEGER KU = ku;
256 const INTEGER M = A.numRows();
257 const INTEGER N = A.numCols();
258 const INTEGER LDA = A.leadingDimension();
259 INTEGER INFO;
260
261 if (IsSame<T,DOUBLE>::value) {
262 LAPACK_IMPL(dlascl)(&TYPE,
263 &KL,
264 &KU,
265 &cFrom,
266 &cTo,
267 &M,
268 &N,
269 A.data(),
270 &LDA,
271 &INFO);
272 } else {
273 ASSERT(0);
274 }
275 ASSERT(INFO==0);
276 }
277
278 #endif // CHECK_CXXLAPACK
279
280 //== public interface ==========================================================
281
282 template <typename IndexType, typename T, typename MA>
283 typename RestrictTo<IsSame<typename MA::ElementType, T>::value, void>::Type
284 lascl(LASCL::Type type, IndexType kl, IndexType ku,
285 const T &cFrom, const T &cTo, MA &A)
286 {
287 LAPACK_DEBUG_OUT("lascl");
288
289 # ifdef CHECK_CXXLAPACK
290 typename MA::NoView _A = A;
291 # endif
292
293 lascl_generic(type, kl, ku, cFrom, cTo, A);
294
295 # ifdef CHECK_CXXLAPACK
296 lascl_native(type, kl, ku, cFrom, cTo, _A);
297 if (! isIdentical(A, _A, " A", "A_")) {
298 std::cerr << "CXXLAPACK: A = " << A << std::endl;
299 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
300 ASSERT(0);
301 }
302 # endif
303 }
304
305 //-- forwarding ----------------------------------------------------------------
306 template <typename IndexType, typename T, typename MA>
307 typename RestrictTo<IsSame<typename MA::ElementType, T>::value, void>::Type
308 lascl(LASCL::Type type, IndexType kl, IndexType ku,
309 const T &cFrom, const T &cTo, MA &&A)
310 {
311 lascl(type, kl, ku, cFrom, cTo, A);
312 }
313
314 //-- convert vector to matrix --------------------------------------------------
315 template <typename IndexType, typename T, typename VX>
316 void
317 lascl(LASCL::Type type, IndexType kl, IndexType ku,
318 const T &cFrom, const T &cTo, DenseVector<VX> &x)
319 {
320 using namespace LASCL;
321
322 typedef typename DenseVector<VX>::ElementType TX;
323 typedef FullStorage<TX, ColMajor> FS;
324
325 typename GeMatrix<FS>::View A(x.length(), 1, x);
326 ASSERT(type==FullMatrix);
327 lascl(type, kl, ku, cFrom, cTo, A);
328 }
329
330 //-- convert scalar to matrix --------------------------------------------------
331 template <typename IndexType, typename T, typename MA>
332 typename RestrictTo<IsSame<MA, T>::value, void>::Type
333 lascl(LASCL::Type type,
334 IndexType kl,
335 IndexType ku,
336 const T &cFrom,
337 const T &cTo,
338 MA &scalar)
339 {
340 using namespace LASCL;
341
342 GeMatrixView<MA> A = typename GeMatrixView<MA>::Engine(1, 1, &scalar, 1);
343
344 ASSERT(type==FullMatrix);
345 lascl(type, kl, ku, cFrom, cTo, A);
346 }
347
348 template <typename IndexType, typename T, typename VX>
349 void
350 lascl(LASCL::Type type, IndexType kl, IndexType ku,
351 const T &cFrom, const T &cTo, DenseVector<VX> &&x)
352 {
353 lascl(type, kl, ku, cFrom, cTo, x);
354 }
355
356 } } // namespace lapack, flens
357
358 #endif // FLENS_LAPACK_AUX_LASCL_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 DLASCL( TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO )
36 *
37 * -- LAPACK auxiliary routine (version 3.3.0) --
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 2010
41 */
42
43 #ifndef FLENS_LAPACK_AUX_LASCL_TCC
44 #define FLENS_LAPACK_AUX_LASCL_TCC 1
45
46 #include <cmath>
47 #include <flens/matrixtypes/matrixtypes.h>
48 #include <flens/vectortypes/vectortypes.h>
49
50 ////////
51 ////////
52 ////////
53 ////////
54 ////////
55 ////////
56 ////////
57 ////////
58 //////// TODO: Completely redesign this crap!
59 ////////
60 ////////
61 ////////
62 //////// ... but it works.
63 ////////
64 ////////
65 ////////
66 ////////
67 ////////
68
69 namespace flens { namespace lapack {
70
71 //== generic lapack implementation =============================================
72
73 template <typename IndexType, typename T, typename MA>
74 typename RestrictTo<IsSame<typename MA::ElementType, T>::value, void>::Type
75 lascl_generic(LASCL::Type type,
76 IndexType kl,
77 IndexType ku,
78 const T &cFrom,
79 const T &cTo,
80 MA &A)
81 {
82 using namespace LASCL;
83
84 using std::abs;
85 using std::isnan;
86 using std::max;
87 using std::min;
88
89 const IndexType m = A.numRows();
90 const IndexType n = A.numCols();
91 const T Zero = 0;
92 const T One = 1;
93
94 if ((cFrom==Zero) || isnan(cFrom)) {
95 ASSERT(0);
96 } else if (isnan(cTo)) {
97 ASSERT(0);
98 } else if (type==SymmetricLowerBand) {
99 ASSERT(A.numRows()==A.numCols());
100 } else if (type==SymmetricUpperBand) {
101 ASSERT(A.numRows()==A.numCols());
102 }
103 //
104 // Quick return if possible
105 //
106 if ((n==0) || (m==0)) {
107 return;
108 }
109 //
110 // Get machine parameters
111 //
112 const T smallNum = lamch<T>(SafeMin);
113 const T bigNum = T(1) / smallNum;
114 //
115 // Make copies of cFrom, cTo
116 //
117 T cFromC = cFrom;
118 T cToC = cTo;
119
120 bool done = false;
121
122 do {
123 T cFrom1 = cFromC*smallNum;
124 T cTo1;
125 T mul;
126 if (cFrom1==cFromC) {
127 // cFromC is an inf. Multiply by a correctly signed zero for
128 // finite cToC, or a NaN if cToC is infinite.
129 mul = cToC / cFromC;
130 done = true;
131 cTo1 = cToC;
132 } else {
133 cTo1 = cToC / bigNum;
134 if (cTo1==cToC) {
135 // cToC is either 0 or an inf. In both cases, cToC itself
136 // serves as the correct multiplication factor.
137 mul = cToC;
138 done = true;
139 cFromC = One;
140 } else if (abs(cFrom1)>abs(cToC) && cToC!=Zero) {
141 mul = smallNum;
142 done = false;
143 cFromC = cFrom1;
144 } else if (abs(cTo1)>abs(cFromC)) {
145 mul = bigNum;
146 done = false;
147 cToC = cTo1;
148 } else {
149 mul = cToC / cFromC;
150 done = true;
151 }
152 }
153
154 if (type==FullMatrix) {
155 //
156 // Full matrix
157 //
158 for (IndexType j=1; j<=n; ++j) {
159 for (IndexType i=1; i<=m; ++i) {
160 A(i,j) *= mul;
161 }
162 }
163 } else if (type==LowerTriangular) {
164 //
165 // Lower triangular matrix
166 //
167 for (IndexType j=1; j<=n; ++j) {
168 for (IndexType i=j; i<=m; ++i) {
169 A(i,j) *= mul;
170 }
171 }
172 } else if (type==UpperTriangular) {
173 //
174 // Upper triangular matrix
175 //
176 for (IndexType j=1; j<=n; ++j) {
177 for (IndexType i=1; i<=min(j,m); ++i) {
178 A(i,j) *= mul;
179 }
180 }
181 } else if (type==UpperHessenberg) {
182 //
183 // Upper Hessenberg matrix
184 //
185 for (IndexType j=1; j<=n; ++j) {
186 for (IndexType i=1; i<=min(j+1,m); ++i) {
187 A(i,j) *= mul;
188 }
189 }
190 } else if (type==SymmetricLowerBand) {
191 //
192 // Lower half of a symmetric band matrix
193 //
194 // TODO: this only works for the internal fullstorage of a
195 // band matrix not for the external element access
196 // of SbMatrix
197 const IndexType k3 = kl + 1;
198 const IndexType k4 = n + 1;
199 for (IndexType j=1; j<=n; ++j) {
200 for (IndexType i=1; i<=min(k3, k4-j); ++i) {
201 A(i,j) *= mul;
202 }
203 }
204 } else if (type==SymmetricUpperBand) {
205 //
206 // Upper half of a symmetric band matrix
207 //
208 // TODO: this only works for the internal fullstorage of a
209 // band matrix not for the external element access
210 // of SbMatrix
211 const IndexType k1 = ku + 2;
212 const IndexType k3 = ku + 1;
213 for (IndexType j=1; j<=n; ++j) {
214 for (IndexType i=max(k1-j,IndexType(1)); i<=k3; ++i) {
215 A(i,j) *= mul;
216 }
217 }
218 } else if (type==GeneralBand) {
219 //
220 // Band matrix
221 //
222 // TODO: this only works for the internal fullstorage of a
223 // band matrix not for the external element access
224 // of GeMatrix
225 const IndexType k1 = kl + ku + 2;
226 const IndexType k2 = kl + 1;
227 const IndexType k3 = 2*kl + ku + 1;
228 const IndexType k4 = kl + ku + 1 + m;
229 for (IndexType j=1; j<=n; ++j) {
230 for (IndexType i=max(k1-j,k2); i<=min(k3,k4-j); ++i) {
231 A(i,j) *= mul;
232 }
233 }
234 } else {
235 ASSERT(0);
236 }
237 } while (!done);
238 }
239
240 //== interface for native lapack ===============================================
241
242 #ifdef CHECK_CXXLAPACK
243
244 template <typename IndexType, typename T, typename MA>
245 void
246 lascl_native(LASCL::Type type,
247 IndexType kl,
248 IndexType ku,
249 const T cFrom,
250 const T &cTo,
251 MA &A)
252 {
253 const char TYPE = getF77LapackChar<LASCL::Type>(type);
254 const INTEGER KL = kl;
255 const INTEGER KU = ku;
256 const INTEGER M = A.numRows();
257 const INTEGER N = A.numCols();
258 const INTEGER LDA = A.leadingDimension();
259 INTEGER INFO;
260
261 if (IsSame<T,DOUBLE>::value) {
262 LAPACK_IMPL(dlascl)(&TYPE,
263 &KL,
264 &KU,
265 &cFrom,
266 &cTo,
267 &M,
268 &N,
269 A.data(),
270 &LDA,
271 &INFO);
272 } else {
273 ASSERT(0);
274 }
275 ASSERT(INFO==0);
276 }
277
278 #endif // CHECK_CXXLAPACK
279
280 //== public interface ==========================================================
281
282 template <typename IndexType, typename T, typename MA>
283 typename RestrictTo<IsSame<typename MA::ElementType, T>::value, void>::Type
284 lascl(LASCL::Type type, IndexType kl, IndexType ku,
285 const T &cFrom, const T &cTo, MA &A)
286 {
287 LAPACK_DEBUG_OUT("lascl");
288
289 # ifdef CHECK_CXXLAPACK
290 typename MA::NoView _A = A;
291 # endif
292
293 lascl_generic(type, kl, ku, cFrom, cTo, A);
294
295 # ifdef CHECK_CXXLAPACK
296 lascl_native(type, kl, ku, cFrom, cTo, _A);
297 if (! isIdentical(A, _A, " A", "A_")) {
298 std::cerr << "CXXLAPACK: A = " << A << std::endl;
299 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
300 ASSERT(0);
301 }
302 # endif
303 }
304
305 //-- forwarding ----------------------------------------------------------------
306 template <typename IndexType, typename T, typename MA>
307 typename RestrictTo<IsSame<typename MA::ElementType, T>::value, void>::Type
308 lascl(LASCL::Type type, IndexType kl, IndexType ku,
309 const T &cFrom, const T &cTo, MA &&A)
310 {
311 lascl(type, kl, ku, cFrom, cTo, A);
312 }
313
314 //-- convert vector to matrix --------------------------------------------------
315 template <typename IndexType, typename T, typename VX>
316 void
317 lascl(LASCL::Type type, IndexType kl, IndexType ku,
318 const T &cFrom, const T &cTo, DenseVector<VX> &x)
319 {
320 using namespace LASCL;
321
322 typedef typename DenseVector<VX>::ElementType TX;
323 typedef FullStorage<TX, ColMajor> FS;
324
325 typename GeMatrix<FS>::View A(x.length(), 1, x);
326 ASSERT(type==FullMatrix);
327 lascl(type, kl, ku, cFrom, cTo, A);
328 }
329
330 //-- convert scalar to matrix --------------------------------------------------
331 template <typename IndexType, typename T, typename MA>
332 typename RestrictTo<IsSame<MA, T>::value, void>::Type
333 lascl(LASCL::Type type,
334 IndexType kl,
335 IndexType ku,
336 const T &cFrom,
337 const T &cTo,
338 MA &scalar)
339 {
340 using namespace LASCL;
341
342 GeMatrixView<MA> A = typename GeMatrixView<MA>::Engine(1, 1, &scalar, 1);
343
344 ASSERT(type==FullMatrix);
345 lascl(type, kl, ku, cFrom, cTo, A);
346 }
347
348 template <typename IndexType, typename T, typename VX>
349 void
350 lascl(LASCL::Type type, IndexType kl, IndexType ku,
351 const T &cFrom, const T &cTo, DenseVector<VX> &&x)
352 {
353 lascl(type, kl, ku, cFrom, cTo, x);
354 }
355
356 } } // namespace lapack, flens
357
358 #endif // FLENS_LAPACK_AUX_LASCL_TCC