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 DLASSQ( N, X, INCX, SCALE, SUMSQ )
36 *
37 * -- LAPACK auxiliary 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_AUX_LASSQ_TCC
44 #define FLENS_LAPACK_AUX_LASSQ_TCC 1
45
46 #include <cmath>
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename VX, typename T>
55 void
56 lassq_generic(const DenseVector<VX> &x, T &scale, T &sumsq)
57 {
58 using std::abs;
59 using std::pow;
60
61 typedef typename DenseVector<VX>::IndexType IndexType;
62
63 const IndexType n = x.length();
64
65 if (n>0) {
66 for (IndexType i=x.firstIndex(); i<=x.lastIndex(); ++i) {
67 if (x(i)!=T(0)) {
68 const T absXi = abs(x(i));
69 if (scale<absXi) {
70 sumsq = 1 + sumsq*pow(scale/absXi, 2);
71 scale = absXi;
72 } else {
73 sumsq += pow(absXi/scale, 2);
74 }
75 }
76 }
77 }
78 }
79
80 //== interface for native lapack ===============================================
81
82 #ifdef CHECK_CXXLAPACK
83
84 template <typename VX, typename T>
85 void
86 lassq_native(const DenseVector<VX> &x, T &scale, T &sumsq)
87 {
88 const INTEGER N = x.length();
89 const INTEGER INCX = x.inc();
90 T SCALE = scale;
91 T SUMSQ = sumsq;
92
93 if (IsSame<T,DOUBLE>::value) {
94 LAPACK_IMPL(dlassq)(&N,
95 x.data(),
96 &INCX,
97 &SCALE,
98 &SUMSQ);
99 } else {
100 ASSERT(0);
101 }
102
103 scale = SCALE;
104 sumsq = SUMSQ;
105 }
106
107 #endif // CHECK_CXXLAPACK
108
109
110 //== public interface ==========================================================
111
112 template <typename VX, typename T>
113 void
114 lassq(const DenseVector<VX> &x, T &scale, T &sumsq)
115 {
116 LAPACK_DEBUG_OUT("lassq");
117
118 //
119 // Test the input parameters
120 //
121 ASSERT(x.inc()>0);
122
123 //
124 // Make copies of output arguments
125 //
126 # ifdef CHECK_CXXLAPACK
127 T _scale = scale;
128 T _sumsq = sumsq;
129 # endif
130
131 //
132 // Call implementation
133 //
134 lassq_generic(x, scale, sumsq);
135
136 # ifdef CHECK_CXXLAPACK
137 //
138 // Compare results
139 //
140 lassq_native(x, _scale, _sumsq);
141
142 bool failed = false;
143 if (! isIdentical(scale, _scale, " scale", "_scale")) {
144 std::cerr << "CXXLAPACK: scale = " << scale << std::endl;
145 std::cerr << "F77LAPACK: _scale = " << _scale << std::endl;
146 failed = true;
147 }
148
149 if (! isIdentical(sumsq, _sumsq, " sumsq", "_sumsq")) {
150 std::cerr << "CXXLAPACK: sumsq = " << sumsq << std::endl;
151 std::cerr << "F77LAPACK: _sumsq = " << _sumsq << std::endl;
152 failed = true;
153 }
154
155 if (failed) {
156 ASSERT(0);
157 } else {
158 // std::cerr << "passed: lassq.tcc" << std::endl;
159 }
160
161 # endif
162 }
163
164 } } // namespace lapack, flens
165
166 #endif // FLENS_LAPACK_AUX_LASSQ_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 DLASSQ( N, X, INCX, SCALE, SUMSQ )
36 *
37 * -- LAPACK auxiliary 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_AUX_LASSQ_TCC
44 #define FLENS_LAPACK_AUX_LASSQ_TCC 1
45
46 #include <cmath>
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename VX, typename T>
55 void
56 lassq_generic(const DenseVector<VX> &x, T &scale, T &sumsq)
57 {
58 using std::abs;
59 using std::pow;
60
61 typedef typename DenseVector<VX>::IndexType IndexType;
62
63 const IndexType n = x.length();
64
65 if (n>0) {
66 for (IndexType i=x.firstIndex(); i<=x.lastIndex(); ++i) {
67 if (x(i)!=T(0)) {
68 const T absXi = abs(x(i));
69 if (scale<absXi) {
70 sumsq = 1 + sumsq*pow(scale/absXi, 2);
71 scale = absXi;
72 } else {
73 sumsq += pow(absXi/scale, 2);
74 }
75 }
76 }
77 }
78 }
79
80 //== interface for native lapack ===============================================
81
82 #ifdef CHECK_CXXLAPACK
83
84 template <typename VX, typename T>
85 void
86 lassq_native(const DenseVector<VX> &x, T &scale, T &sumsq)
87 {
88 const INTEGER N = x.length();
89 const INTEGER INCX = x.inc();
90 T SCALE = scale;
91 T SUMSQ = sumsq;
92
93 if (IsSame<T,DOUBLE>::value) {
94 LAPACK_IMPL(dlassq)(&N,
95 x.data(),
96 &INCX,
97 &SCALE,
98 &SUMSQ);
99 } else {
100 ASSERT(0);
101 }
102
103 scale = SCALE;
104 sumsq = SUMSQ;
105 }
106
107 #endif // CHECK_CXXLAPACK
108
109
110 //== public interface ==========================================================
111
112 template <typename VX, typename T>
113 void
114 lassq(const DenseVector<VX> &x, T &scale, T &sumsq)
115 {
116 LAPACK_DEBUG_OUT("lassq");
117
118 //
119 // Test the input parameters
120 //
121 ASSERT(x.inc()>0);
122
123 //
124 // Make copies of output arguments
125 //
126 # ifdef CHECK_CXXLAPACK
127 T _scale = scale;
128 T _sumsq = sumsq;
129 # endif
130
131 //
132 // Call implementation
133 //
134 lassq_generic(x, scale, sumsq);
135
136 # ifdef CHECK_CXXLAPACK
137 //
138 // Compare results
139 //
140 lassq_native(x, _scale, _sumsq);
141
142 bool failed = false;
143 if (! isIdentical(scale, _scale, " scale", "_scale")) {
144 std::cerr << "CXXLAPACK: scale = " << scale << std::endl;
145 std::cerr << "F77LAPACK: _scale = " << _scale << std::endl;
146 failed = true;
147 }
148
149 if (! isIdentical(sumsq, _sumsq, " sumsq", "_sumsq")) {
150 std::cerr << "CXXLAPACK: sumsq = " << sumsq << std::endl;
151 std::cerr << "F77LAPACK: _sumsq = " << _sumsq << std::endl;
152 failed = true;
153 }
154
155 if (failed) {
156 ASSERT(0);
157 } else {
158 // std::cerr << "passed: lassq.tcc" << std::endl;
159 }
160
161 # endif
162 }
163
164 } } // namespace lapack, flens
165
166 #endif // FLENS_LAPACK_AUX_LASSQ_TCC