| 29 | |
| 30 | template<typename Tr> |
| 31 | Tr calc_scale(Tr From, Tr To) { |
| 32 | // FIXME: I am not sure this is correct, removing this for now |
| 33 | #if 0 |
| 34 | //http://www.netlib.org/lapack/explore-3.1.1-html/dlascl.f.html |
| 35 | cpu_lapack_lamch_func<Tr> cpu_lapack_lamch; |
| 36 | |
| 37 | Tr S = cpu_lapack_lamch('S'); |
| 38 | Tr B = 1.0 / S; |
| 39 | |
| 40 | Tr FromCopy = From, ToCopy = To; |
| 41 | |
| 42 | Tr Mul = 1; |
| 43 | |
| 44 | while (true) { |
| 45 | Tr From1 = FromCopy * S, To1 = ToCopy / B; |
| 46 | if (std::abs(From1) > std::abs(ToCopy) && ToCopy != 0) { |
| 47 | Mul *= S; |
| 48 | FromCopy = From1; |
| 49 | } else if (std::abs(To1) > std::abs(FromCopy)) { |
| 50 | Mul *= B; |
| 51 | ToCopy = To1; |
| 52 | } else { |
| 53 | Mul *= (ToCopy) / (FromCopy); |
| 54 | break; |
| 55 | } |
| 56 | } |
| 57 | |
| 58 | return Mul; |
| 59 | #else |
| 60 | return To / From; |
| 61 | #endif |
| 62 | } |
| 63 | |
| 64 | template<typename T, typename Tr> |
| 65 | void svd(Array<T> &arrU, Array<Tr> &arrS, Array<T> &arrVT, Array<T> &arrA, |