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       SUBROUTINE DLABAD( SMALL, LARGE )
 35  *
 36  *  -- LAPACK auxiliary routine (version 3.2) --
 37  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 38  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 39  *     November 2006
 40  */
 41 
 42 #ifndef FLENS_LAPACK_AUX_LABAD_TCC
 43 #define FLENS_LAPACK_AUX_LABAD_TCC 1
 44 
 45 #include <cmath>
 46 #include <flens/lapack/lapack.h>
 47 
 48 namespace flens { namespace lapack {
 49 
 50 //== generic lapack implementation =============================================
 51 
 52 template <typename T>
 53 void
 54 labad_generic(T &small, T &large)
 55 {
 56 //
 57 //  If it looks like we're on a Cray, take the square root of
 58 //  SMALL and LARGE to avoid overflow and underflow problems.
 59 //
 60     if (log10(large)>2000) {
 61         small = sqrt(small);
 62         large = sqrt(large);
 63     }
 64 }
 65 
 66 //== interface for native lapack ===============================================
 67 
 68 #ifdef CHECK_CXXLAPACK
 69 
 70 template <typename T>
 71 void
 72 labad_native(T &small, T &large)
 73 {
 74     if (IsSame<T, DOUBLE>::value) {
 75         LAPACK_IMPL(dlabad)(&small,
 76                             &large);
 77     } else {
 78         ASSERT(0);
 79     }
 80 }
 81 
 82 #endif // CHECK_CXXLAPACK
 83 
 84 //== public interface ==========================================================
 85 
 86 template <typename T>
 87 void
 88 labad(T &small, T &large)
 89 {
 90     LAPACK_DEBUG_OUT("BEGIN: labad");
 91 
 92 #   ifdef CHECK_CXXLAPACK
 93 //
 94 //  Make copies of output arguments
 95 //
 96     T _small    = small;
 97     T _large    = large;
 98 #   endif
 99 
100 //
101 //  Call implementation
102 //
103     labad_generic(small, large);
104 
105 #   ifdef CHECK_CXXLAPACK
106 //
107 //  Compare results
108 //
109     labad_native(_small, _large);
110 
111     bool failed = false;
112     if (! isIdentical(small, _small, " small""_small")) {
113         std::cerr << "CXXLAPACK:  small = " << small << std::endl;
114         std::cerr << "F77LAPACK: _small = " << _small << std::endl;
115         failed = true;
116     }
117 
118     if (! isIdentical(large, _large, " large""_large")) {
119         std::cerr << "CXXLAPACK:  large = " << large << std::endl;
120         std::cerr << "F77LAPACK: _large = " << _large << std::endl;
121         failed = true;
122     }
123 
124     if (failed) {
125         ASSERT(0);
126     }
127 #   endif
128 
129     LAPACK_DEBUG_OUT("END: labad");
130 }
131 
132 } } // namespace lapack, flens
133 
134 #endif // FLENS_LAPACK_AUX_LABAD_TCC