[Boost-commit] svn:boost r71114 - in sandbox/numeric_bindings: boost/numeric/bindings/lapack/driver libs/numeric/bindings/tools/templates/driver

Subject: [Boost-commit] svn:boost r71114 - in sandbox/numeric_bindings: boost/numeric/bindings/lapack/driver libs/numeric/bindings/tools/templates/driver
From: rutger_at_[hidden]
Date: 2011-04-08 04:55:56

Author: rutger
Date: 2011-04-08 04:55:55 EDT (Fri, 08 Apr 2011)
New Revision: 71114
URL: http://svn.boost.org/trac/boost/changeset/71114

Log:
Fixed the complex case for gesvd. 

Text files modified: 
   sandbox/numeric_bindings/boost/numeric/bindings/lapack/driver/gesvd.hpp         |    41 ++++++++++++++++++++++++++++++--------- 
   sandbox/numeric_bindings/libs/numeric/bindings/tools/templates/driver/gesvd.hpp |    30 ++++++++++++++++++++++++----            
   2 files changed, 56 insertions(+), 15 deletions(-)

Modified: sandbox/numeric_bindings/boost/numeric/bindings/lapack/driver/gesvd.hpp
==============================================================================
--- sandbox/numeric_bindings/boost/numeric/bindings/lapack/driver/gesvd.hpp	(original)
+++ sandbox/numeric_bindings/boost/numeric/bindings/lapack/driver/gesvd.hpp	2011-04-08 04:55:55 EDT (Fri, 08 Apr 2011)
@@ -245,17 +245,19 @@
         // Contributed by Marco Guazzone
         // Also see http://tinyurl.com/5rbpdc5
         //
-        if ( m >= n ) {
+        if ( m == 0 || n == 0 ) {
+            return 1;
+        } else if ( m >= n ) {
             if ( jobu == 'N' ) {
                 return 5*n;
             } else {
-        	return std::max< std::ptrdiff_t >(3*n+m,5*n);
+                return std::max< std::ptrdiff_t >(3*n+m,5*n);
             }
         } else {
             if ( jobvt == 'N' ) {
                 return 5*m;
             } else {
-        	return std::max< std::ptrdiff_t >(3*m+n,5*m);
+                return std::max< std::ptrdiff_t >(3*m+n,5*m);
             }
         }
     }
@@ -304,7 +306,7 @@
         BOOST_ASSERT( bindings::size(work.select(real_type())) >=
                 min_size_rwork( minmn ));
         BOOST_ASSERT( bindings::size(work.select(value_type())) >=
-                min_size_work( bindings::size_row(a),
+                min_size_work( jobu, jobvt, bindings::size_row(a),
                 bindings::size_column(a), minmn ));
         BOOST_ASSERT( bindings::size_column(a) >= 0 );
         BOOST_ASSERT( bindings::size_minor(a) == 1 ||
@@ -345,8 +347,9 @@
         namespace bindings = ::boost::numeric::bindings;
         std::ptrdiff_t minmn = std::min< std::ptrdiff_t >( size_row(a),
                 size_column(a) );
-        bindings::detail::array< value_type > tmp_work( min_size_work(
-                bindings::size_row(a), bindings::size_column(a), minmn ) );
+        bindings::detail::array< value_type > tmp_work( min_size_work( jobu,
+                jobvt, bindings::size_row(a), bindings::size_column(a),
+                minmn ) );
         bindings::detail::array< real_type > tmp_rwork( min_size_rwork(
                 minmn ) );
         return invoke( jobu, jobvt, a, s, u, vt, workspace( tmp_work,
@@ -387,10 +390,28 @@
     // Static member function that returns the minimum size of
     // workspace-array work.
     //
-    static std::ptrdiff_t min_size_work( const std::ptrdiff_t m,
-            const std::ptrdiff_t n, const std::ptrdiff_t minmn ) {
-        return std::max< std::ptrdiff_t >( 1, 2*minmn+std::max<
-                std::ptrdiff_t >(m,n) );
+    static std::ptrdiff_t min_size_work( const char jobu, const char jobvt,
+            const std::ptrdiff_t m, const std::ptrdiff_t n,
+            const std::ptrdiff_t minmn ) {
+        //
+        // Contributed by Marco Guazzone
+        // Also see http://tinyurl.com/5rbpdc5
+        //
+        if ( minmn == 0 ) {
+            return 1;
+        } else if ( m >= n ) {
+            if ( jobu == 'N' ) {
+                return 3*n;
+            } else {
+                return 2*n+m;
+            }
+        } else {
+            if ( jobvt == 'N' ) {
+                return 3*m;
+            } else {
+                return 2*m+n;
+            }
+        }
     }
 
     //

Modified: sandbox/numeric_bindings/libs/numeric/bindings/tools/templates/driver/gesvd.hpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/tools/templates/driver/gesvd.hpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/tools/templates/driver/gesvd.hpp	2011-04-08 04:55:55 EDT (Fri, 08 Apr 2011)
@@ -5,17 +5,19 @@
 // Contributed by Marco Guazzone
 // Also see http://tinyurl.com/5rbpdc5
 //
-if ( m >= n ) {
+if ( m == 0 || n == 0 ) {
+    return 1;
+} else if ( m >= n ) {
     if ( jobu == 'N' ) {
         return 5*n;
     } else {
-	return std::max< $INTEGER_TYPE >(3*n+m,5*n);
+        return std::max< $INTEGER_TYPE >(3*n+m,5*n);
     }
 } else {
     if ( jobvt == 'N' ) {
         return 5*m;
     } else {
-	return std::max< $INTEGER_TYPE >(3*m+n,5*m);
+        return std::max< $INTEGER_TYPE >(3*m+n,5*m);
     }
 }
 $TEMPLATE[gesvd.complex.extra_variables]
@@ -25,9 +27,27 @@
 $TEMPLATE[gesvd.complex.MINMN.init]
 $INTEGER_TYPE minmn = std::min< $INTEGER_TYPE >( size_row(a), size_column(a) );
 $TEMPLATE[gesvd.complex.min_size_work.args]
-M,N,MINMN
+JOBU,JOBVT,M,N,MINMN
 $TEMPLATE[gesvd.complex.min_size_work]
-return std::max< $INTEGER_TYPE >( 1, 2*minmn+std::max< $INTEGER_TYPE >(m,n) );
+//
+// Contributed by Marco Guazzone
+// Also see http://tinyurl.com/5rbpdc5
+//
+if ( minmn == 0 ) {
+    return 1;
+} else if ( m >= n ) {
+    if ( jobu == 'N' ) {
+        return 3*n;
+    } else {
+        return 2*n+m;
+    }
+} else {
+    if ( jobvt == 'N' ) {
+        return 3*m;
+    } else {
+        return 2*m+n;
+    }
+}
 $TEMPLATE[gesvd.complex.min_size_rwork.args]
 MINMN
 $TEMPLATE[gesvd.complex.min_size_rwork]