Subject: [Boost-commit] svn:boost r50177 - sandbox/libs/numeric/bindings/lapack/test
From: thomas.klimpel_at_[hidden]
Date: 2008-12-07 08:20:33
Author: klimpel
Date: 2008-12-07 08:20:32 EST (Sun, 07 Dec 2008)
New Revision: 50177
URL: http://svn.boost.org/trac/boost/changeset/50177
Log:
added Jesse Manning's gels bindings
Added:
sandbox/libs/numeric/bindings/lapack/test/convenience.h (contents, props changed)
sandbox/libs/numeric/bindings/lapack/test/ublas_gels.cpp (contents, props changed)
sandbox/libs/numeric/bindings/lapack/test/ublas_gelsd.cpp (contents, props changed)
sandbox/libs/numeric/bindings/lapack/test/ublas_gelss.cpp (contents, props changed)
Text files modified:
sandbox/libs/numeric/bindings/lapack/test/Jamfile.v2 | 4 ++++
1 files changed, 4 insertions(+), 0 deletions(-)
Modified: sandbox/libs/numeric/bindings/lapack/test/Jamfile.v2
==============================================================================
--- sandbox/libs/numeric/bindings/lapack/test/Jamfile.v2 (original)
+++ sandbox/libs/numeric/bindings/lapack/test/Jamfile.v2 2008-12-07 08:20:32 EST (Sun, 07 Dec 2008)
@@ -53,5 +53,9 @@
[ run ublas_sptrf_sptrs.cc : 20 ]
[ run ublas_sytrf_sytrs.cc : 20 ]
[ run ublas_gbsv.cpp ]
+
+ [ run ublas_gels.cpp ]
+ [ run ublas_gelss.cpp ]
+ [ run ublas_gelsd.cpp ]
;
Added: sandbox/libs/numeric/bindings/lapack/test/convenience.h
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/convenience.h 2008-12-07 08:20:32 EST (Sun, 07 Dec 2008)
@@ -0,0 +1,266 @@
+//
+// Copyright Jesse Manning 2007
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#ifndef CONVENIENCE_H
+#define CONVENIENCE_H
+
+#include <iostream>
+#include <fstream>
+#include <complex>
+#include <string>
+#include <iomanip>
+#include <sstream>
+
+#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
+#include <boost/numeric/bindings/traits/ublas_vector.hpp>
+#include <boost/numeric/bindings/traits/matrix_traits.hpp>
+#include <boost/numeric/bindings/traits/vector_traits.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+
+// included to implicitly convert a vector to an nx1 matrix
+// so that it is compatible with lapack binding
+#include <boost/numeric/bindings/traits/ublas_vector2.hpp>
+
+namespace ublas = boost::numeric::ublas;
+namespace traits = boost::numeric::bindings::traits;
+
+// single precision typedefs
+typedef float freal_t;
+typedef ublas::matrix<freal_t, ublas::column_major> fmat_t;
+typedef ublas::vector<freal_t> fvec_t;
+
+// double precision typedefs
+typedef double dreal_t;
+typedef ublas::matrix<dreal_t, ublas::column_major> dmat_t;
+typedef ublas::vector<dreal_t> dvec_t;
+
+// single precision complex typedefs
+typedef std::complex<float> fcmplx_t;
+typedef ublas::matrix<fcmplx_t, ublas::column_major> fcmat_t;
+typedef ublas::vector<fcmplx_t> fcvec_t;
+
+// double precision complex typedefs
+typedef std::complex<double> dcmplx_t;
+typedef ublas::matrix<dcmplx_t, ublas::column_major> dcmat_t;
+typedef ublas::vector<dcmplx_t> dcvec_t;
+
+// matrix/vector test sizes
+const int row_size = 3;
+const int col_size = 3;
+const int row_range = 2;
+const int col_range = 2;
+
+//////////////////////////////////////////////////////////////////
+//
+// Helper functions and structs to aid with testing
+//
+/////////////////////////////////////////////////////////////////
+template <typename StreamType, typename MatType>
+void matrix_print(StreamType& oss, const std::string& name, const MatType& mat)
+{
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+ oss << name << std::endl;
+ for (int i=0; i < m; ++i)
+ {
+ for (int j=0; j < n; ++j)
+ {
+ oss << mat(i,j) << std::setw(10);
+ }
+ oss << std::setw(0) << std::endl;
+ }
+}
+
+template <typename StreamType, typename VecType>
+void vector_print(StreamType& oss, const std::string& name, const VecType& vec)
+{
+ const int m = traits::vector_size(vec);
+
+ oss << name << std::endl;
+ for (int i=0; i < m; ++i)
+ {
+ oss << vec(i) << std::endl;
+ }
+}
+
+// structs to create matrices for testing
+template <typename T>
+struct MatrixGenerator {};
+
+template <>
+struct MatrixGenerator<fmat_t>
+{
+ typedef fmat_t Result;
+ inline Result operator() (size_t m, size_t n)
+ {
+ Result mat(row_size, col_size);
+ mat(0,0) = 8.;
+ mat(0,1) = 1.;
+ mat(0,2) = 6.;
+ mat(1,0) = 3.;
+ mat(1,1) = 5.;
+ mat(1,2) = 7.;
+ mat(2,0) = 4.;
+ mat(2,1) = 9.;
+ mat(2,2) = 2.;
+
+ return Result(ublas::project(mat, ublas::range(0, m), ublas::range(0, n)));
+ }
+};
+
+template <>
+struct MatrixGenerator<dmat_t>
+{
+ typedef dmat_t Result;
+ inline Result operator() (size_t m, size_t n)
+ {
+ Result mat(row_size, col_size);
+ mat(0,0) = 8.;
+ mat(0,1) = 1.;
+ mat(0,2) = 6.;
+ mat(1,0) = 3.;
+ mat(1,1) = 5.;
+ mat(1,2) = 7.;
+ mat(2,0) = 4.;
+ mat(2,1) = 9.;
+ mat(2,2) = 2.;
+
+ return Result(ublas::project(mat, ublas::range(0, m), ublas::range(0, n)));
+ }
+};
+
+template <>
+struct MatrixGenerator<fcmat_t>
+{
+ typedef fcmat_t Result;
+ inline Result operator() (size_t m, size_t n)
+ {
+ typedef Result::value_type val_t;
+
+ Result mat(row_size , col_size);
+ mat(0,0) = val_t(35.,1.);
+ mat(0,1) = val_t(6.,26.);
+ mat(0,2) = val_t(19.,24.);
+ mat(1,0) = val_t(3.,32.);
+ mat(1,1) = val_t(7.,21.);
+ mat(1,2) = val_t(23.,25.);
+ mat(2,0) = val_t(31.,9.);
+ mat(2,1) = val_t(2.,22.);
+ mat(2,2) = val_t(27.,20.);
+
+ return Result(ublas::project(mat, ublas::range(0, m), ublas::range(0, n)));
+ }
+};
+
+template <>
+struct MatrixGenerator<dcmat_t>
+{
+ typedef dcmat_t Result;
+ inline Result operator() (size_t m, size_t n)
+ {
+ typedef Result::value_type val_t;
+
+ Result mat(row_size , col_size);
+ mat(0,0) = val_t(35.,1.);
+ mat(0,1) = val_t(6.,26.);
+ mat(0,2) = val_t(19.,24.);
+ mat(1,0) = val_t(3.,32.);
+ mat(1,1) = val_t(7.,21.);
+ mat(1,2) = val_t(23.,25.);
+ mat(2,0) = val_t(31.,9.);
+ mat(2,1) = val_t(2.,22.);
+ mat(2,2) = val_t(27.,20.);
+
+ return Result(ublas::project(mat, ublas::range(0, m), ublas::range(0, n)));
+ }
+};
+
+
+// structs to create vectors for testing
+template <typename T>
+struct VectorGenerator {};
+
+template <>
+struct VectorGenerator<fvec_t>
+{
+ typedef fvec_t Result;
+ inline Result operator() (size_t m)
+ {
+ typedef Result::value_type val_t;
+ Result v(m);
+
+ for (size_t i=0; i < m; ++i)
+ {
+ val_t val = val_t(i+1);
+ v(i) = val;
+ }
+
+ return v;
+ }
+};
+
+template <>
+struct VectorGenerator<dvec_t>
+{
+ typedef dvec_t Result;
+ inline Result operator() (size_t m)
+ {
+ typedef Result::value_type val_t;
+ Result v(m);
+
+ for (size_t i=0; i < m; ++i)
+ {
+ val_t val = val_t(i+1);
+ v(i) = val;
+ }
+
+ return v;
+ }
+};
+
+template <>
+struct VectorGenerator<fcvec_t>
+{
+ typedef fcvec_t Result;
+ inline Result operator() (size_t m)
+ {
+ typedef Result::value_type val_t;
+ Result v(m);
+
+ for (size_t i=0; i < m; ++i)
+ {
+ val_t::value_type val = val_t::value_type(i);
+ v(i) = val_t(val+1, val+1);
+ }
+
+ return v;
+ }
+};
+
+template <>
+struct VectorGenerator<dcvec_t>
+{
+ typedef dcvec_t Result;
+ inline Result operator() (size_t m)
+ {
+ typedef Result::value_type val_t;
+ Result v(m);
+
+ for (size_t i=0; i < m; ++i)
+ {
+ val_t::value_type val = val_t::value_type(i);
+ v(i) = val_t(val+1, val+1);
+ }
+
+ return v;
+ }
+};
+
+#endif
Added: sandbox/libs/numeric/bindings/lapack/test/ublas_gels.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/ublas_gels.cpp 2008-12-07 08:20:32 EST (Sun, 07 Dec 2008)
@@ -0,0 +1,527 @@
+//
+// Copyright Jesse Manning 2007
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#include "convenience.h"
+#include <boost/numeric/bindings/lapack/gels.hpp>
+
+// set to 1 to write test output to file, otherwise outputs to console
+#define OUTPUT_TO_FILE 0
+
+// determines which tests to run
+#define TEST_SQUARE 1
+#define TEST_UNDERDETERMINED 1
+#define TEST_OVERDETERMINED 1
+#define TEST_MULTIPLE_SOLUTION_VECTORS 1
+#define TEST_TRANSPOSE 1
+
+// determines if optimal, minimal, or both workspaces are used for testing
+#define USE_OPTIMAL_WORKSPACE 1
+#define USE_MINIMAL_WORKSPACE 1
+
+namespace lapack = boost::numeric::bindings::lapack;
+
+// test function declarations
+template <typename StreamType, typename MatrType, typename VecType>
+int test_square_gels(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_under_gels(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_over_gels(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_multiple_gels(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_transpose_gels(StreamType& oss, const char& trans);
+
+int main()
+{
+ // stream for test output
+ typedef std::ostringstream stream_t;
+ stream_t oss;
+
+#if TEST_SQUARE
+ oss << "Start Square Matrix Least Squares Tests" << std::endl;
+ oss << "Testing sgels" << std::endl;
+ if (test_square_gels<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgels passed." << std::endl;
+ }
+ oss << "End sgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgels" << std::endl;
+ if (test_square_gels<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgels passed." << std::endl;
+ }
+ oss << "End dgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgels" << std::endl;
+ if (test_square_gels<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgels passed." << std::endl;
+ }
+ oss << "End cgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgels" << std::endl;
+ if (test_square_gels<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgels passed." << std::endl;
+ }
+ oss << "End zgels tests" << std::endl;
+ oss << std::endl;
+ oss << "End Square Matrix Least Squares Tests" << std::endl;
+#endif
+
+#if TEST_UNDERDETERMINED
+ oss << std::endl;
+ oss << "Start Under-determined Matrix Least Squares Test" << std::endl;
+ oss << "Testing sgels" << std::endl;
+ if (test_under_gels<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgels passed." << std::endl;
+ }
+ oss << "End sgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgels" << std::endl;
+ if (test_under_gels<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgels passed." << std::endl;
+ }
+ oss << "End dgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgels" << std::endl;
+ if (test_under_gels<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgels passed." << std::endl;
+ }
+ oss << "End cgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgels" << std::endl;
+ if (test_under_gels<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgels passed." << std::endl;
+ }
+ oss << "End zgels tests" << std::endl;
+ oss << std::endl;
+ oss << "End Underdetermined Matrix Least Squares Tests" << std::endl;
+#endif
+
+#if TEST_OVERDETERMINED
+ oss << std::endl;
+ oss << "Start Overdetermined Matrix Least Squares Test" << std::endl;
+ oss << "Testing sgels" << std::endl;
+ if (test_over_gels<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgels passed." << std::endl;
+ }
+ oss << "End sgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgels" << std::endl;
+ if (test_over_gels<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgels passed." << std::endl;
+ }
+ oss << "End dgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgels" << std::endl;
+ if (test_over_gels<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgels passed." << std::endl;
+ }
+ oss << "End cgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgels" << std::endl;
+ if (test_over_gels<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgels passed." << std::endl;
+ }
+ oss << "End zgels tests" << std::endl;
+ oss << std::endl;
+ oss << "End Overdetermined Matrix Least Squares Test" << std::endl;
+#endif
+
+#if TEST_MULTIPLE_SOLUTION_VECTORS
+ oss << std::endl;
+ oss << "Start Multiple Solution Vectors Least Squares Test" << std::endl;
+ oss << "Testing sgels" << std::endl;
+ if (test_multiple_gels<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgels passed." << std::endl;
+ }
+ oss << "End sgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgels" << std::endl;
+ if (test_multiple_gels<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgels passed." << std::endl;
+ }
+ oss << "End dgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgels" << std::endl;
+ if (test_multiple_gels<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgels passed." << std::endl;
+ }
+ oss << "End cgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgels" << std::endl;
+ if (test_multiple_gels<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgels passed." << std::endl;
+ }
+ oss << "End zgels tests" << std::endl;
+ oss << std::endl;
+ oss << "End Multiple Solution Vectors Least Squares Test" << std::endl;
+#endif
+
+#if TEST_TRANSPOSE
+ oss << std::endl;
+ oss << "Start Transpose Least Squares Test" << std::endl;
+ oss << "Testing sgels" << std::endl;
+ if (test_transpose_gels<stream_t, fmat_t, fvec_t>(oss, 'T') == 0)
+ {
+ oss << "sgels passed." << std::endl;
+ }
+ oss << "End sgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgels" << std::endl;
+ if (test_transpose_gels<stream_t, dmat_t, dvec_t>(oss, 'T') == 0)
+ {
+ oss << "dgels passed." << std::endl;
+ }
+ oss << "End dgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgels" << std::endl;
+ if (test_transpose_gels<stream_t, fcmat_t, fcvec_t>(oss, 'C') == 0)
+ {
+ oss << "cgels passed." << std::endl;
+ }
+ oss << "End cgels tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgels" << std::endl;
+ if (test_transpose_gels<stream_t, fcmat_t, fcvec_t>(oss, 'C') == 0)
+ {
+ oss << "zgels passed." << std::endl;
+ }
+ oss << "End zgels tests" << std::endl;
+ oss << std::endl;
+ oss << "End Transpose Least Squares Test" << std::endl;
+#endif
+
+ // Finished testing
+ std::cout << std::endl;
+ std::cout << "Tests Completed." << std::endl;
+ std::cout << std::endl;
+
+#if OUTPUT_TO_FILE
+ std::string filename;
+ std::cout << "Enter filename to write test results: ";
+ std::getline(std::cin, filename);
+
+ std::ofstream testFile(filename.c_str());
+
+ if (testFile)
+ {
+ testFile << oss.str();
+ testFile.close();
+ }
+#else
+ std::cout << oss.str() << std::endl;
+#endif
+
+ // wait for user to finish
+// std::string done;
+// std::cout << "Press Enter to exit";
+// std::getline(std::cin, done);
+
+}
+
+// tests square system (m-by-n where m == n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_square_gels(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // square matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_size));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gels('N', optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gels('N', minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+
+ // check A*x=B
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests overdetermined system (m-by-n where m < n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_under_gels(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // under-determined matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_range, col_size));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gels('N', optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gels('N', minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+
+ // check A*x=B
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests overdetermined system (m-by-n where m > n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_over_gels(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // overdetermined matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_range));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gels('N', optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gels('N', minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+
+ // check A*x=B
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests multiple solution vectors stored column-wise in B for equation A*x=B
+template <typename StreamType, typename MatType, typename VecType>
+int test_multiple_gels(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // multiple solutions vectors test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_size));
+ MatType vec(mat.size1(), 2);
+ ublas::column(vec, 0) = VectorGenerator<VecType>()(mat.size1());
+ ublas::column(vec, 1) = VectorGenerator<VecType>()(mat.size1());
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+ const int nrhs = traits::matrix_size2(vec);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ MatType optimalvec(vec);
+ err += lapack::gels('N', optimalmat, optimalvec, lapack::optimal_workspace());
+ MatType optimalanswer(ublas::project(optimalvec, ublas::range(0, n), ublas::range(0, nrhs)));
+ MatType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ MatType minimalvec(vec);
+ err += lapack::gels('N', minimalmat, minimalvec, lapack::minimal_workspace());
+ MatType minimalanswer(ublas::project(minimalvec, ublas::range(0, n), ublas::range(0, nrhs)));
+ MatType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ matrix_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ matrix_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ matrix_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+
+ // check A*x=B
+#if USE_OPTIMAL_WORKSPACE
+ matrix_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ matrix_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_transpose_gels(StreamType& oss, const char& trans)
+{
+ // return value
+ int err = 0;
+
+ // overdetermined matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_size));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gels(trans, optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gels(trans, minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+
+ // check A*x=B
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
Added: sandbox/libs/numeric/bindings/lapack/test/ublas_gelsd.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/ublas_gelsd.cpp 2008-12-07 08:20:32 EST (Sun, 07 Dec 2008)
@@ -0,0 +1,439 @@
+//
+// Copyright Jesse Manning 2007
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#include "convenience.h"
+#include <boost/numeric/bindings/lapack/gelsd.hpp>
+
+// set to 1 to write test output to file, otherwise outputs to console
+#define OUTPUT_TO_FILE 0
+
+// determines which tests to run
+#define TEST_SQUARE 1
+#define TEST_UNDERDETERMINED 1
+#define TEST_OVERDETERMINED 1
+#define TEST_MULTIPLE_SOLUTION_VECTORS 1
+
+// determines if optimal, minimal, or both workspaces are used for testing
+#define USE_OPTIMAL_WORKSPACE 1
+#define USE_MINIMAL_WORKSPACE 1
+
+namespace lapack = boost::numeric::bindings::lapack;
+
+// test function declarations
+template <typename StreamType, typename MatrType, typename VecType>
+int test_square_gelsd(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_under_gelsd(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_over_gelsd(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_multiple_gelsd(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_transpose_gel(StreamType& oss, const char& trans);
+
+
+int main()
+{
+ // stream for test output
+ typedef std::ostringstream stream_t;
+ stream_t oss;
+
+#if TEST_SQUARE
+ oss << "Start Square Matrix Least Squares Tests" << std::endl;
+ oss << "Testing sgelsd" << std::endl;
+ if (test_square_gelsd<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgelsd passed." << std::endl;
+ }
+ oss << "End sgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgelsd" << std::endl;
+ if (test_square_gelsd<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgelsd passed." << std::endl;
+ }
+ oss << "End dgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgelsd" << std::endl;
+ if (test_square_gelsd<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgelsd passed." << std::endl;
+ }
+ oss << "End cgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgelsd" << std::endl;
+ if (test_square_gelsd<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgelsd passed." << std::endl;
+ }
+ oss << "End zgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "End Square Matrix Least Squares Tests" << std::endl;
+#endif
+
+#if TEST_UNDERDETERMINED
+ oss << std::endl;
+ oss << "Start Under-determined Matrix Least Squares Test" << std::endl;
+ oss << "Testing sgelsd" << std::endl;
+ if (test_under_gelsd<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgelsd passed." << std::endl;
+ }
+ oss << "End sgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgelsd" << std::endl;
+ if (test_under_gelsd<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgelsd passed." << std::endl;
+ }
+ oss << "End dgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgelsd" << std::endl;
+ if (test_under_gelsd<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgelsd passed." << std::endl;
+ }
+ oss << "End cgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgelsd" << std::endl;
+ if (test_under_gelsd<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgelsd passed." << std::endl;
+ }
+ oss << "End zgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "End Underdetermined Matrix Least Squares Tests" << std::endl;
+#endif
+
+#if TEST_OVERDETERMINED
+ oss << std::endl;
+ oss << "Start Overdetermined Matrix Least Squares Test" << std::endl;
+ oss << "Testing sgelsd" << std::endl;
+ if (test_over_gelsd<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgelsd passed." << std::endl;
+ }
+ oss << "End sgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgelsd" << std::endl;
+ if (test_over_gelsd<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgelsd passed." << std::endl;
+ }
+ oss << "End dgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgelsd" << std::endl;
+ if (test_over_gelsd<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgelsd passed." << std::endl;
+ }
+ oss << "End cgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgelsd" << std::endl;
+ if (test_over_gelsd<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgelsd passed." << std::endl;
+ }
+ oss << "End zgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "End Overdetermined Matrix Least Squares Test" << std::endl;
+#endif
+
+#if TEST_MULTIPLE_SOLUTION_VECTORS
+ oss << std::endl;
+ oss << "Start Multiple Solution Vectors Least Squares Test" << std::endl;
+ oss << "Testing sgelsd" << std::endl;
+ if (test_multiple_gelsd<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgelsd passed." << std::endl;
+ }
+ oss << "End sgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgelsd" << std::endl;
+ if (test_multiple_gelsd<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgelsd passed." << std::endl;
+ }
+ oss << "End dgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgelsd" << std::endl;
+ if (test_multiple_gelsd<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgelsd passed." << std::endl;
+ }
+ oss << "End cgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgelsd" << std::endl;
+ if (test_multiple_gelsd<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgelsd passed." << std::endl;
+ }
+ oss << "End zgelsd tests" << std::endl;
+ oss << std::endl;
+ oss << "End Multiple Solution Vectors Least Squares Test" << std::endl;
+#endif
+
+#if OUTPUT_TO_FILE
+ // Finished testing
+ std::cout << std::endl;
+ std::cout << "Tests Completed." << std::endl;
+ std::cout << std::endl;
+
+ std::string filename;
+ std::cout << "Enter filename to write test results: ";
+ std::getline(std::cin, filename);
+
+ std::ofstream testFile(filename.c_str());
+
+ if (testFile)
+ {
+ testFile << oss.str();
+ testFile.close();
+ }
+#else
+ std::cout << oss.str() << std::endl;
+
+ // Finished testing
+ std::cout << std::endl;
+ std::cout << "Tests Completed." << std::endl;
+ std::cout << std::endl;
+#endif
+
+ // wait for user to finish
+// std::string done;
+// std::cout << "Press Enter to exit";
+// std::getline(std::cin, done);
+
+}
+
+// tests square system (m-by-n where m == n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_square_gelsd(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // square matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_size));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gelsd(optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gelsd(minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+#if USE_OPTIMAL_WORKSPACE
+ // check A*x=B
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests overdetermined system (m-by-n where m < n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_under_gelsd(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // under-determined matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_range, col_size));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gelsd(optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gelsd(minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+#if USE_OPTIMAL_WORKSPACE
+ // check A*x=B
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests overdetermined system (m-by-n where m > n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_over_gelsd(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // overdetermined matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_range));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gelsd(optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gelsd(minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+#if USE_OPTIMAL_WORKSPACE
+ // check A*x=B
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests multiple solution vectors stored column-wise in B for equation A*x=B
+template <typename StreamType, typename MatType, typename VecType>
+int test_multiple_gelsd(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // multiple solutions vectors test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_size));
+ MatType vec(mat.size1(), 2);
+ ublas::column(vec, 0) = VectorGenerator<VecType>()(mat.size1());
+ ublas::column(vec, 1) = VectorGenerator<VecType>()(mat.size1());
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+ const int nrhs = traits::matrix_size2(vec);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ MatType optimalvec(vec);
+ err += lapack::gelsd(optimalmat, optimalvec, lapack::optimal_workspace());
+ MatType optimalanswer(ublas::project(optimalvec, ublas::range(0, n), ublas::range(0, nrhs)));
+ MatType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ MatType minimalvec(vec);
+ err += lapack::gelsd(minimalmat, minimalvec, lapack::minimal_workspace());
+ MatType minimalanswer(ublas::project(minimalvec, ublas::range(0, n), ublas::range(0, nrhs)));
+ MatType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ matrix_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ matrix_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ matrix_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+#if USE_OPTIMAL_WORKSPACE
+ // check A*x=B
+ matrix_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ matrix_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
Added: sandbox/libs/numeric/bindings/lapack/test/ublas_gelss.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/ublas_gelss.cpp 2008-12-07 08:20:32 EST (Sun, 07 Dec 2008)
@@ -0,0 +1,439 @@
+//
+// Copyright Jesse Manning 2007
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#include "convenience.h"
+#include <boost/numeric/bindings/lapack/gelss.hpp>
+
+// set to 1 to write test output to file, otherwise outputs to console
+#define OUTPUT_TO_FILE 0
+
+// determines which tests to run
+#define TEST_SQUARE 1
+#define TEST_UNDERDETERMINED 1
+#define TEST_OVERDETERMINED 1
+#define TEST_MULTIPLE_SOLUTION_VECTORS 1
+
+// determines if optimal, minimal, or both workspaces are used for testing
+#define USE_OPTIMAL_WORKSPACE 1
+#define USE_MINIMAL_WORKSPACE 1
+
+namespace lapack = boost::numeric::bindings::lapack;
+
+// test function declarations
+template <typename StreamType, typename MatrType, typename VecType>
+int test_square_gelss(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_under_gelss(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_over_gelss(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_multiple_gelss(StreamType& oss);
+
+template <typename StreamType, typename MatType, typename VecType>
+int test_transpose_gel(StreamType& oss, const char& trans);
+
+
+int main()
+{
+ // stream for test output
+ typedef std::ostringstream stream_t;
+ stream_t oss;
+
+#if TEST_SQUARE
+ oss << "Start Square Matrix Least Squares Tests" << std::endl;
+ oss << "Testing sgelss" << std::endl;
+ if (test_square_gelss<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgelss passed." << std::endl;
+ }
+ oss << "End sgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgelss" << std::endl;
+ if (test_square_gelss<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgelss passed." << std::endl;
+ }
+ oss << "End dgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgelss" << std::endl;
+ if (test_square_gelss<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgelss passed." << std::endl;
+ }
+ oss << "End cgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgelss" << std::endl;
+ if (test_square_gelss<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgelss passed." << std::endl;
+ }
+ oss << "End zgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "End Square Matrix Least Squares Tests" << std::endl;
+#endif
+
+#if TEST_UNDERDETERMINED
+ oss << std::endl;
+ oss << "Start Under-determined Matrix Least Squares Test" << std::endl;
+ oss << "Testing sgelss" << std::endl;
+ if (test_under_gelss<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgelss passed." << std::endl;
+ }
+ oss << "End sgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgelss" << std::endl;
+ if (test_under_gelss<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgelss passed." << std::endl;
+ }
+ oss << "End dgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgelss" << std::endl;
+ if (test_under_gelss<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgelss passed." << std::endl;
+ }
+ oss << "End cgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgelss" << std::endl;
+ if (test_under_gelss<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgelss passed." << std::endl;
+ }
+ oss << "End zgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "End Underdetermined Matrix Least Squares Tests" << std::endl;
+#endif
+
+#if TEST_OVERDETERMINED
+ oss << std::endl;
+ oss << "Start Overdetermined Matrix Least Squares Test" << std::endl;
+ oss << "Testing sgelss" << std::endl;
+ if (test_over_gelss<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgelss passed." << std::endl;
+ }
+ oss << "End sgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgelss" << std::endl;
+ if (test_over_gelss<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgelss passed." << std::endl;
+ }
+ oss << "End dgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgelss" << std::endl;
+ if (test_over_gelss<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgelss passed." << std::endl;
+ }
+ oss << "End cgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgelss" << std::endl;
+ if (test_over_gelss<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgelss passed." << std::endl;
+ }
+ oss << "End zgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "End Overdetermined Matrix Least Squares Test" << std::endl;
+#endif
+
+#if TEST_MULTIPLE_SOLUTION_VECTORS
+ oss << std::endl;
+ oss << "Start Multiple Solution Vectors Least Squares Test" << std::endl;
+ oss << "Testing sgelss" << std::endl;
+ if (test_multiple_gelss<stream_t, fmat_t, fvec_t>(oss) == 0)
+ {
+ oss << "sgelss passed." << std::endl;
+ }
+ oss << "End sgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing dgelss" << std::endl;
+ if (test_multiple_gelss<stream_t, dmat_t, dvec_t>(oss) == 0)
+ {
+ oss << "dgelss passed." << std::endl;
+ }
+ oss << "End dgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing cgelss" << std::endl;
+ if (test_multiple_gelss<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "cgelss passed." << std::endl;
+ }
+ oss << "End cgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "Testing zgelss" << std::endl;
+ if (test_multiple_gelss<stream_t, fcmat_t, fcvec_t>(oss) == 0)
+ {
+ oss << "zgelss passed." << std::endl;
+ }
+ oss << "End zgelss tests" << std::endl;
+ oss << std::endl;
+ oss << "End Multiple Solution Vectors Least Squares Test" << std::endl;
+#endif
+
+#if OUTPUT_TO_FILE
+ // Finished testing
+ std::cout << std::endl;
+ std::cout << "Tests Completed." << std::endl;
+ std::cout << std::endl;
+
+ std::string filename;
+ std::cout << "Enter filename to write test results: ";
+ std::getline(std::cin, filename);
+
+ std::ofstream testFile(filename.c_str());
+
+ if (testFile)
+ {
+ testFile << oss.str();
+ testFile.close();
+ }
+#else
+ std::cout << oss.str() << std::endl;
+
+ // Finished testing
+ std::cout << std::endl;
+ std::cout << "Tests Completed." << std::endl;
+ std::cout << std::endl;
+#endif
+
+ // wait for user to finish
+// std::string done;
+// std::cout << "Press Enter to exit";
+// std::getline(std::cin, done);
+
+}
+
+// tests square system (m-by-n where m == n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_square_gelss(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // square matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_size));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gelss(optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gelss(minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+#if USE_OPTIMAL_WORKSPACE
+ // check A*x=B
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests overdetermined system (m-by-n where m < n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_under_gelss(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // under-determined matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_range, col_size));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gelss(optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gelss(minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+#if USE_OPTIMAL_WORKSPACE
+ // check A*x=B
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests overdetermined system (m-by-n where m > n)
+template <typename StreamType, typename MatType, typename VecType>
+int test_over_gelss(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // overdetermined matrix test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_range));
+ VecType vec(VectorGenerator<VecType>()(row_size));
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ VecType optimalvec(vec);
+ err += lapack::gelss(optimalmat, optimalvec, lapack::optimal_workspace());
+ VecType optimalanswer(ublas::project(optimalvec, ublas::range(0, n)));
+ VecType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ VecType minimalvec(vec);
+ err += lapack::gelss(minimalmat, minimalvec, lapack::minimal_workspace());
+ VecType minimalanswer(ublas::project(minimalvec, ublas::range(0, n)));
+ VecType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ vector_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ vector_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+#if USE_OPTIMAL_WORKSPACE
+ // check A*x=B
+ vector_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ vector_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}
+
+// tests multiple solution vectors stored column-wise in B for equation A*x=B
+template <typename StreamType, typename MatType, typename VecType>
+int test_multiple_gelss(StreamType& oss)
+{
+ // return value
+ int err = 0;
+
+ // multiple solutions vectors test
+ MatType mat(MatrixGenerator<MatType>()(row_size, col_size));
+ MatType vec(mat.size1(), 2);
+ ublas::column(vec, 0) = VectorGenerator<VecType>()(mat.size1());
+ ublas::column(vec, 1) = VectorGenerator<VecType>()(mat.size1());
+
+ const int m = traits::matrix_size1(mat);
+ const int n = traits::matrix_size2(mat);
+ const int nrhs = traits::matrix_size2(vec);
+
+#if USE_OPTIMAL_WORKSPACE
+ MatType optimalmat(mat);
+ MatType optimalvec(vec);
+ err += lapack::gelss(optimalmat, optimalvec, lapack::optimal_workspace());
+ MatType optimalanswer(ublas::project(optimalvec, ublas::range(0, n), ublas::range(0, nrhs)));
+ MatType optimal_check = ublas::prod(mat, optimalanswer);
+#endif
+#if USE_MINIMAL_WORKSPACE
+ MatType minimalmat(mat);
+ MatType minimalvec(vec);
+ err += lapack::gelss(minimalmat, minimalvec, lapack::minimal_workspace());
+ MatType minimalanswer(ublas::project(minimalvec, ublas::range(0, n), ublas::range(0, nrhs)));
+ MatType minimal_check = ublas::prod(mat, minimalanswer);
+#endif
+
+ matrix_print(oss, "A", mat);
+ oss << std::endl;
+ matrix_print(oss, "B", vec);
+ oss << std::endl;
+
+#if USE_OPTIMAL_WORKSPACE
+ matrix_print(oss, "optimal workspace x", optimalanswer);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ matrix_print(oss, "minimal workspace x", minimalanswer);
+ oss << std::endl;
+#endif
+#if USE_OPTIMAL_WORKSPACE
+ // check A*x=B
+ matrix_print(oss, "optimal A*x=B", optimal_check);
+ oss << std::endl;
+#endif
+#if USE_MINIMAL_WORKSPACE
+ matrix_print(oss, "minimal A*x=B", minimal_check);
+ oss << std::endl;
+#endif
+
+ return err;
+}