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 /* Baesed on
34 *
35 SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
36 *
37 * -- LAPACK auxiliary routine (version 3.3.1) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * -- April 2011 --
41 */
42
43 #ifndef FLENS_LAPACK_AUX_LARFG_TCC
44 #define FLENS_LAPACK_AUX_LARFG_TCC 1
45
46 #include <flens/blas/blas.h>
47 #include <flens/lapack/lapack.h>
48
49 namespace flens { namespace lapack {
50
51 //== generic lapack implementation =============================================
52
53 template <typename IndexType, typename ALPHA, typename VX, typename TAU>
54 void
55 larfg_generic(IndexType n, ALPHA &alpha, DenseVector<VX> &x, TAU &tau)
56 {
57 using std::abs;
58
59 typedef typename DenseVector<VX>::ElementType T;
60
61 if (n<=1) {
62 tau = TAU(0);
63 return;
64 }
65
66 T xNorm = blas::nrm2(x);
67 if (xNorm==T(0)) {
68 //
69 // H = I
70 //
71 tau = TAU(0);
72 } else {
73 //
74 // general case
75 //
76 T beta = -sign(lapy2(alpha, xNorm), alpha);
77 T safeMin = lamch<T>(SafeMin) / lamch<T>(Eps);
78
79 IndexType count=0;
80 if (abs(beta)<safeMin) {
81 //
82 // XNORM, BETA may be inaccurate; scale X and recompute them
83 //
84 T rSafeMin = T(1)/safeMin;
85 do {
86 ++count;
87 blas::scal(rSafeMin, x);
88 beta *= rSafeMin;
89 alpha *= rSafeMin;
90 } while (abs(beta)<safeMin);
91 //
92 // New BETA is at most 1, at least SAFMIN
93 //
94 xNorm = blas::nrm2(x);
95 beta = -sign(lapy2(alpha, xNorm), alpha);
96 }
97 tau = (beta-alpha) / beta;
98 blas::scal(T(1)/(alpha-beta), x);
99 //
100 // If ALPHA is subnormal, it may lose relative accuracy
101 //
102 for (IndexType j=1; j<=count; ++j) {
103 beta *= safeMin;
104 }
105 alpha = beta;
106 }
107 }
108
109 //== interface for native lapack ===============================================
110
111 #ifdef CHECK_CXXLAPACK
112
113 template <typename IndexType, typename ALPHA, typename VX, typename TAU>
114 void
115 larfg_native(IndexType n, ALPHA &alpha, DenseVector<VX> &x, TAU &tau)
116 {
117 typedef typename DenseVector<VX>::ElementType T;
118
119 INTEGER N = n;
120 INTEGER INCX = x.inc();
121
122 if (IsSame<T, DOUBLE>::value) {
123 LAPACK_IMPL(dlarfg)(&N,
124 &alpha,
125 x.data(),
126 &INCX,
127 &tau);
128 } else {
129 ASSERT(0);
130 }
131 }
132
133 #endif // CHECK_CXXLAPACK
134
135 //== public interface ==========================================================
136
137 template <typename IndexType, typename ALPHA, typename VX, typename TAU>
138 void
139 larfg(IndexType n, ALPHA &alpha, DenseVector<VX> &x, TAU &tau)
140 {
141 LAPACK_DEBUG_OUT("larfg");
142
143 //
144 // Test the input parameters
145 //
146 ASSERT(x.firstIndex()==1);
147 ASSERT(x.inc()>0);
148 ASSERT(x.length()<=n);
149
150 # ifdef CHECK_CXXLAPACK
151 //
152 // Make copies of output arguments
153 //
154 ALPHA _alpha = alpha;
155 typename DenseVector<VX>::NoView _x = x;
156 TAU _tau = tau;
157 # endif
158
159 //
160 // Call implementation
161 //
162 larfg_generic(n, alpha, x, tau);
163
164 # ifdef CHECK_CXXLAPACK
165 //
166 // Compare results
167 //
168 larfg_native(n, _alpha, _x, _tau);
169
170 bool failed = false;
171 if (alpha!=_alpha) {
172 std::cerr << "CXXLAPACK: alpha = " << alpha << std::endl;
173 std::cerr << "F77LAPACK: _alpha = " << _alpha << std::endl;
174 std::cerr << "CXXLAPACK: alpha - _alpha= "
175 << alpha - _alpha << std::endl;
176 failed = true;
177 }
178
179 if (! isIdentical(x, _x, " x", "x_")) {
180 std::cerr << "CXXLAPACK: x = " << x << std::endl;
181 std::cerr << "F77LAPACK: _x = " << _x << std::endl;
182 failed = true;
183 }
184
185 if (tau!=_tau) {
186 std::cerr << "CXXLAPACK: tau = " << tau << std::endl;
187 std::cerr << "F77LAPACK: _tau = " << _tau << std::endl;
188 std::cerr << "CXXLAPACK: tau - _tau= "
189 << tau - _tau << std::endl;
190 failed = true;
191 }
192
193 if (failed) {
194 ASSERT(0);
195 }
196 # endif
197 }
198
199 //-- forwarding ----------------------------------------------------------------
200 template <typename IndexType, typename ALPHA, typename VX, typename TAU>
201 void
202 larfg(IndexType n, ALPHA &alpha, VX &&x, TAU &tau)
203 {
204 larfg(n, alpha, x, tau);
205 }
206
207 } } // namespace lapack, flens
208
209 #endif // FLENS_LAPACK_AUX_LARFG_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 /* Baesed on
34 *
35 SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
36 *
37 * -- LAPACK auxiliary routine (version 3.3.1) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * -- April 2011 --
41 */
42
43 #ifndef FLENS_LAPACK_AUX_LARFG_TCC
44 #define FLENS_LAPACK_AUX_LARFG_TCC 1
45
46 #include <flens/blas/blas.h>
47 #include <flens/lapack/lapack.h>
48
49 namespace flens { namespace lapack {
50
51 //== generic lapack implementation =============================================
52
53 template <typename IndexType, typename ALPHA, typename VX, typename TAU>
54 void
55 larfg_generic(IndexType n, ALPHA &alpha, DenseVector<VX> &x, TAU &tau)
56 {
57 using std::abs;
58
59 typedef typename DenseVector<VX>::ElementType T;
60
61 if (n<=1) {
62 tau = TAU(0);
63 return;
64 }
65
66 T xNorm = blas::nrm2(x);
67 if (xNorm==T(0)) {
68 //
69 // H = I
70 //
71 tau = TAU(0);
72 } else {
73 //
74 // general case
75 //
76 T beta = -sign(lapy2(alpha, xNorm), alpha);
77 T safeMin = lamch<T>(SafeMin) / lamch<T>(Eps);
78
79 IndexType count=0;
80 if (abs(beta)<safeMin) {
81 //
82 // XNORM, BETA may be inaccurate; scale X and recompute them
83 //
84 T rSafeMin = T(1)/safeMin;
85 do {
86 ++count;
87 blas::scal(rSafeMin, x);
88 beta *= rSafeMin;
89 alpha *= rSafeMin;
90 } while (abs(beta)<safeMin);
91 //
92 // New BETA is at most 1, at least SAFMIN
93 //
94 xNorm = blas::nrm2(x);
95 beta = -sign(lapy2(alpha, xNorm), alpha);
96 }
97 tau = (beta-alpha) / beta;
98 blas::scal(T(1)/(alpha-beta), x);
99 //
100 // If ALPHA is subnormal, it may lose relative accuracy
101 //
102 for (IndexType j=1; j<=count; ++j) {
103 beta *= safeMin;
104 }
105 alpha = beta;
106 }
107 }
108
109 //== interface for native lapack ===============================================
110
111 #ifdef CHECK_CXXLAPACK
112
113 template <typename IndexType, typename ALPHA, typename VX, typename TAU>
114 void
115 larfg_native(IndexType n, ALPHA &alpha, DenseVector<VX> &x, TAU &tau)
116 {
117 typedef typename DenseVector<VX>::ElementType T;
118
119 INTEGER N = n;
120 INTEGER INCX = x.inc();
121
122 if (IsSame<T, DOUBLE>::value) {
123 LAPACK_IMPL(dlarfg)(&N,
124 &alpha,
125 x.data(),
126 &INCX,
127 &tau);
128 } else {
129 ASSERT(0);
130 }
131 }
132
133 #endif // CHECK_CXXLAPACK
134
135 //== public interface ==========================================================
136
137 template <typename IndexType, typename ALPHA, typename VX, typename TAU>
138 void
139 larfg(IndexType n, ALPHA &alpha, DenseVector<VX> &x, TAU &tau)
140 {
141 LAPACK_DEBUG_OUT("larfg");
142
143 //
144 // Test the input parameters
145 //
146 ASSERT(x.firstIndex()==1);
147 ASSERT(x.inc()>0);
148 ASSERT(x.length()<=n);
149
150 # ifdef CHECK_CXXLAPACK
151 //
152 // Make copies of output arguments
153 //
154 ALPHA _alpha = alpha;
155 typename DenseVector<VX>::NoView _x = x;
156 TAU _tau = tau;
157 # endif
158
159 //
160 // Call implementation
161 //
162 larfg_generic(n, alpha, x, tau);
163
164 # ifdef CHECK_CXXLAPACK
165 //
166 // Compare results
167 //
168 larfg_native(n, _alpha, _x, _tau);
169
170 bool failed = false;
171 if (alpha!=_alpha) {
172 std::cerr << "CXXLAPACK: alpha = " << alpha << std::endl;
173 std::cerr << "F77LAPACK: _alpha = " << _alpha << std::endl;
174 std::cerr << "CXXLAPACK: alpha - _alpha= "
175 << alpha - _alpha << std::endl;
176 failed = true;
177 }
178
179 if (! isIdentical(x, _x, " x", "x_")) {
180 std::cerr << "CXXLAPACK: x = " << x << std::endl;
181 std::cerr << "F77LAPACK: _x = " << _x << std::endl;
182 failed = true;
183 }
184
185 if (tau!=_tau) {
186 std::cerr << "CXXLAPACK: tau = " << tau << std::endl;
187 std::cerr << "F77LAPACK: _tau = " << _tau << std::endl;
188 std::cerr << "CXXLAPACK: tau - _tau= "
189 << tau - _tau << std::endl;
190 failed = true;
191 }
192
193 if (failed) {
194 ASSERT(0);
195 }
196 # endif
197 }
198
199 //-- forwarding ----------------------------------------------------------------
200 template <typename IndexType, typename ALPHA, typename VX, typename TAU>
201 void
202 larfg(IndexType n, ALPHA &alpha, VX &&x, TAU &tau)
203 {
204 larfg(n, alpha, x, tau);
205 }
206
207 } } // namespace lapack, flens
208
209 #endif // FLENS_LAPACK_AUX_LARFG_TCC