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 DRSCL( N, SA, SX, INCX )
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_EIG_RSCL_TCC
44 #define FLENS_LAPACK_EIG_RSCL_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 template <typename SA, typename VSX>
53 void
54 rscl_generic(const SA &sa, DenseVector<VSX> &sx)
55 {
56 using std::abs;
57
58 typedef typename DenseVector<VSX>::ElementType ElementType;
59 typedef typename DenseVector<VSX>::IndexType IndexType;
60
61 const IndexType n = sx.length();
62
63 const ElementType Zero(0), One(1);
64 //
65 // Quick return if possible
66 //
67 if (n==0) {
68 return;
69 }
70 //
71 // Get machine parameters
72 //
73 ElementType smallNum = lamch<ElementType>(SafeMin);
74 ElementType bigNum = One / smallNum;
75 labad(smallNum, bigNum);
76 //
77 // Initialize the denominator to SA and the numerator to 1.
78 //
79 ElementType cDen = sa;
80 ElementType cNum = One;
81
82 bool done = false;
83 do {
84 const ElementType cDen1 = cDen*smallNum;
85 const ElementType cNum1 = cNum / bigNum;
86
87 ElementType mul;
88
89 if (abs(cDen1)>abs(cNum) && cNum!=Zero) {
90 //
91 // Pre-multiply X by SMLNUM if CDEN is large compared to CNUM.
92 //
93 mul = smallNum;
94 done = false;
95 cDen = cDen1;
96 } else if (abs(cNum1)>abs(cDen)) {
97 //
98 // Pre-multiply X by BIGNUM if CDEN is small compared to CNUM.
99 //
100 mul = bigNum;
101 done = false;
102 cNum = cNum1;
103 } else {
104 //
105 // Multiply X by CNUM / CDEN and return.
106 //
107 mul = cNum / cDen;
108 done = true;
109 }
110 //
111 // Scale the vector X by MUL
112 //
113 sx *= mul;
114
115 } while (!done);
116 }
117
118 //== interface for native lapack ===============================================
119
120 #ifdef CHECK_CXXLAPACK
121 template <typename SA, typename VSX>
122 void
123 rscl_native(const SA &sa, DenseVector<VSX> &sx)
124 {
125 typedef typename DenseVector<VSX>::ElementType T;
126
127 const INTEGER N = sx.length();
128 const T _SA = sa;
129 const INTEGER INCX = sx.stride();
130
131 if (IsSame<T,double>::value) {
132 LAPACK_IMPL(drscl)(&N, &_SA, sx.data(), &INCX);
133 } else {
134 ASSERT(0);
135 }
136 }
137
138 #endif // CHECK_CXXLAPACK
139
140 //== public interface ==========================================================
141 template <typename SA, typename VSX>
142 void
143 rscl(const SA &sa, DenseVector<VSX> &sx)
144 {
145 //
146 // Test the input parameters
147 //
148 # ifndef NDEBUG
149 ASSERT(sx.firstIndex()==1);
150 # endif
151
152 # ifdef CHECK_CXXLAPACK
153 //
154 // Make copies of output arguments
155 //
156 typename DenseVector<VSX>::NoView sx_org = sx;
157 # endif
158
159 //
160 // Call implementation
161 //
162 rscl_generic(sa, sx);
163
164 # ifdef CHECK_CXXLAPACK
165 //
166 // Make copies of results computed by the generic implementation
167 //
168 typename DenseVector<VSX>::NoView sx_generic = sx;
169
170 //
171 // restore output arguments
172 //
173 sx = sx_org;
174
175 //
176 // Compare generic results with results from the native implementation
177 //
178 rscl_native(sa, sx);
179
180 bool failed = false;
181 if (! isIdentical(sx_generic, sx, "sx_generic", "sx")) {
182 std::cerr << "CXXLAPACK: sx_generic = " << sx_generic << std::endl;
183 std::cerr << "F77LAPACK: sx = " << sx << std::endl;
184 failed = true;
185 }
186
187 if (failed) {
188 ASSERT(0);
189 }
190 # endif
191 }
192
193 //-- forwarding ----------------------------------------------------------------
194 template <typename SA, typename VSX>
195 void
196 rscl(const SA &sa, VSX &&sx)
197 {
198 CHECKPOINT_ENTER;
199 rscl(sa, sx);
200 CHECKPOINT_LEAVE;
201 }
202
203 } } // namespace lapack, flens
204
205 #endif // FLENS_LAPACK_EIG_RSCL_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 DRSCL( N, SA, SX, INCX )
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_EIG_RSCL_TCC
44 #define FLENS_LAPACK_EIG_RSCL_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 template <typename SA, typename VSX>
53 void
54 rscl_generic(const SA &sa, DenseVector<VSX> &sx)
55 {
56 using std::abs;
57
58 typedef typename DenseVector<VSX>::ElementType ElementType;
59 typedef typename DenseVector<VSX>::IndexType IndexType;
60
61 const IndexType n = sx.length();
62
63 const ElementType Zero(0), One(1);
64 //
65 // Quick return if possible
66 //
67 if (n==0) {
68 return;
69 }
70 //
71 // Get machine parameters
72 //
73 ElementType smallNum = lamch<ElementType>(SafeMin);
74 ElementType bigNum = One / smallNum;
75 labad(smallNum, bigNum);
76 //
77 // Initialize the denominator to SA and the numerator to 1.
78 //
79 ElementType cDen = sa;
80 ElementType cNum = One;
81
82 bool done = false;
83 do {
84 const ElementType cDen1 = cDen*smallNum;
85 const ElementType cNum1 = cNum / bigNum;
86
87 ElementType mul;
88
89 if (abs(cDen1)>abs(cNum) && cNum!=Zero) {
90 //
91 // Pre-multiply X by SMLNUM if CDEN is large compared to CNUM.
92 //
93 mul = smallNum;
94 done = false;
95 cDen = cDen1;
96 } else if (abs(cNum1)>abs(cDen)) {
97 //
98 // Pre-multiply X by BIGNUM if CDEN is small compared to CNUM.
99 //
100 mul = bigNum;
101 done = false;
102 cNum = cNum1;
103 } else {
104 //
105 // Multiply X by CNUM / CDEN and return.
106 //
107 mul = cNum / cDen;
108 done = true;
109 }
110 //
111 // Scale the vector X by MUL
112 //
113 sx *= mul;
114
115 } while (!done);
116 }
117
118 //== interface for native lapack ===============================================
119
120 #ifdef CHECK_CXXLAPACK
121 template <typename SA, typename VSX>
122 void
123 rscl_native(const SA &sa, DenseVector<VSX> &sx)
124 {
125 typedef typename DenseVector<VSX>::ElementType T;
126
127 const INTEGER N = sx.length();
128 const T _SA = sa;
129 const INTEGER INCX = sx.stride();
130
131 if (IsSame<T,double>::value) {
132 LAPACK_IMPL(drscl)(&N, &_SA, sx.data(), &INCX);
133 } else {
134 ASSERT(0);
135 }
136 }
137
138 #endif // CHECK_CXXLAPACK
139
140 //== public interface ==========================================================
141 template <typename SA, typename VSX>
142 void
143 rscl(const SA &sa, DenseVector<VSX> &sx)
144 {
145 //
146 // Test the input parameters
147 //
148 # ifndef NDEBUG
149 ASSERT(sx.firstIndex()==1);
150 # endif
151
152 # ifdef CHECK_CXXLAPACK
153 //
154 // Make copies of output arguments
155 //
156 typename DenseVector<VSX>::NoView sx_org = sx;
157 # endif
158
159 //
160 // Call implementation
161 //
162 rscl_generic(sa, sx);
163
164 # ifdef CHECK_CXXLAPACK
165 //
166 // Make copies of results computed by the generic implementation
167 //
168 typename DenseVector<VSX>::NoView sx_generic = sx;
169
170 //
171 // restore output arguments
172 //
173 sx = sx_org;
174
175 //
176 // Compare generic results with results from the native implementation
177 //
178 rscl_native(sa, sx);
179
180 bool failed = false;
181 if (! isIdentical(sx_generic, sx, "sx_generic", "sx")) {
182 std::cerr << "CXXLAPACK: sx_generic = " << sx_generic << std::endl;
183 std::cerr << "F77LAPACK: sx = " << sx << std::endl;
184 failed = true;
185 }
186
187 if (failed) {
188 ASSERT(0);
189 }
190 # endif
191 }
192
193 //-- forwarding ----------------------------------------------------------------
194 template <typename SA, typename VSX>
195 void
196 rscl(const SA &sa, VSX &&sx)
197 {
198 CHECKPOINT_ENTER;
199 rscl(sa, sx);
200 CHECKPOINT_LEAVE;
201 }
202
203 } } // namespace lapack, flens
204
205 #endif // FLENS_LAPACK_EIG_RSCL_TCC