dlib C++ Library - matrix3.cpp

// Copyright (C) 2006 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#include <dlib/matrix.h>
#include <sstream>
#include <string>
#include <cstdlib>
#include <ctime>
#include <vector>
#include "../stl_checked.h"
#include "../array.h"
#include "../rand.h"
#include "tester.h"
#include <dlib/memory_manager_stateless.h>
#include <dlib/array2d.h>
namespace 
{
 using namespace test;
 using namespace dlib;
 using namespace std;
 logger dlog("test.matrix3");
 const double eps_mul = 200;
 template <typename T, typename U>
 void check_equal (
 const T& a,
 const U& b
 )
 {
 DLIB_TEST(a.nr() == b.nr());
 DLIB_TEST(a.nc() == b.nc());
 typedef typename T::type type;
 for (long r = 0; r < a.nr(); ++r)
 {
 for (long c = 0; c < a.nc(); ++c)
 {
 type error = std::abs(a(r,c) - b(r,c));
 DLIB_TEST_MSG(error < std::sqrt(std::numeric_limits<type>::epsilon())*eps_mul, "error: " << error <<
 " eps: " << std::sqrt(std::numeric_limits<type>::epsilon())*eps_mul);
 }
 }
 }
 template <typename T, typename U>
 void c_check_equal (
 const T& a,
 const U& b
 )
 {
 DLIB_TEST(a.nr() == b.nr());
 DLIB_TEST(a.nc() == b.nc());
 typedef typename T::type type;
 for (long r = 0; r < a.nr(); ++r)
 {
 for (long c = 0; c < a.nc(); ++c)
 {
 typename type::value_type error = std::abs(a(r,c) - b(r,c));
 DLIB_TEST_MSG(error < std::sqrt(std::numeric_limits<typename type::value_type>::epsilon())*eps_mul, "error: " << error <<
 " eps: " << std::sqrt(std::numeric_limits<typename type::value_type>::epsilon())*eps_mul);
 }
 }
 }
 template <typename T, typename U>
 void assign_no_blas (
 const T& a_,
 const U& b
 )
 {
 T& a = const_cast<T&>(a_);
 DLIB_TEST(a.nr() == b.nr());
 DLIB_TEST(a.nc() == b.nc());
 for (long r = 0; r < a.nr(); ++r)
 {
 for (long c = 0; c < a.nc(); ++c)
 {
 a(r,c) = b(r,c);
 }
 }
 }
 template <typename type>
 type rnd_num (dlib::rand& rnd)
 {
 return static_cast<type>(10*rnd.get_random_double());
 }
 template <typename type>
 void test_blas( long rows, long cols)
 {
 // The tests in this function exercise the BLAS bindings located in the matrix/matrix_blas_bindings.h file.
 // It does this by performing an assignment that is subject to BLAS bindings and comparing the
 // results directly to an unevaluated matrix_exp that should be equal.
 dlib::rand rnd;
 matrix<type> a(rows,cols), temp, temp2, temp3;
 for (int k = 0; k < 6; ++k)
 {
 for (long r= 0; r < a.nr(); ++r)
 {
 for (long c = 0; c < a.nc(); ++c)
 {
 a(r,c) = rnd_num<type>(rnd);
 }
 }
 matrix<type> at;
 at = trans(a);
 matrix<complex<type> > c_a(rows,cols), c_at, c_sqr;
 for (long r= 0; r < a.nr(); ++r)
 {
 for (long c = 0; c < a.nc(); ++c)
 {
 c_a(r,c) = complex<type>(rnd_num<type>(rnd),rnd_num<type>(rnd));
 }
 }
 c_at = trans(c_a);
 const int size = max(rows,cols);
 c_sqr = 10*matrix_cast<complex<type> >(complex_matrix(randm(size,size,rnd), randm(size,size,rnd)));
 matrix<complex<type> > c_temp(cols,cols), c_temp2(cols,cols);
 const complex<type> i(0,1);
 const type one = 1;
 const type two = 1;
 const type num1 = static_cast<type>(3.6);
 const type num2 = static_cast<type>(6.6);
 const type num3 = static_cast<type>(8.6);
 matrix<complex<type>,0,1> c_cv4(cols), c_cv3(rows);
 matrix<complex<type>,1,0> c_rv4(cols), c_rv3(rows);
 matrix<type,0,1> cv4(cols);
 for (long idx = 0; idx < cv4.size(); ++idx)
 cv4(idx) = rnd_num<type>(rnd);
 for (long idx = 0; idx < c_cv4.size(); ++idx)
 c_cv4(idx) = complex<type>(rnd_num<type>(rnd),rnd_num<type>(rnd));
 matrix<type,1,0> rv3(rows);
 for (long idx = 0; idx < rv3.size(); ++idx)
 rv3(idx) = rnd_num<type>(rnd);
 for (long idx = 0; idx < c_rv3.size(); ++idx)
 c_rv3(idx) = complex<type>(rnd_num<type>(rnd),rnd_num<type>(rnd));
 matrix<type,0,1> cv3(rows);
 for (long idx = 0; idx < cv3.size(); ++idx)
 cv3(idx) = rnd_num<type>(rnd);
 for (long idx = 0; idx < c_cv3.size(); ++idx)
 c_cv3(idx) = complex<type>(rnd_num<type>(rnd),rnd_num<type>(rnd));
 matrix<type,1,0> rv4(cols);
 for (long idx = 0; idx < rv4.size(); ++idx)
 rv4(idx) = rnd_num<type>(rnd);
 for (long idx = 0; idx < c_rv4.size(); ++idx)
 c_rv4(idx) = complex<type>(rnd_num<type>(rnd),rnd_num<type>(rnd));
 // GEMM tests
 dlog << LTRACE << "1.1";
 check_equal(tmp(at*a), at*a);
 check_equal(tmp(trans(at*a)), trans(at*a));
 check_equal(tmp(2.4*trans(4*trans(at*a) + at*3*a)), 2.4*trans(4*trans(at*a) + at*3*a));
 dlog << LTRACE << "1.2";
 check_equal(tmp(trans(a)*a), trans(a)*a);
 check_equal(tmp(trans(trans(a)*a)), trans(trans(a)*a));
 dlog << LTRACE << "1.3";
 check_equal(tmp(at*trans(at)), at*trans(at));
 check_equal(tmp(trans(at*trans(at))), trans(at*trans(at)));
 dlog << LTRACE << "1.4";
 check_equal(tmp(trans(at)*trans(a)), a*at);
 check_equal(tmp(trans(trans(at)*trans(a))), trans(a*at));
 dlog << LTRACE << "1.5";
 print_spinner();
 c_check_equal(tmp(conj(trans(c_a))*c_a), trans(conj(c_a))*c_a);
 dlog << LTRACE << "1.5.1";
 c_check_equal(tmp(trans(conj(trans(c_a))*c_a)), trans(trans(conj(c_a))*c_a));
 dlog << LTRACE << "1.5.2";
 c_check_equal(tmp((conj(trans(c_sqr))*trans(c_sqr))), (trans(conj(c_sqr))*trans(c_sqr)));
 dlog << LTRACE << "1.5.3";
 c_check_equal(tmp(trans(conj(trans(c_sqr))*trans(c_sqr))), trans(trans(conj(c_sqr))*trans(c_sqr)));
 dlog << LTRACE << "1.6";
 c_check_equal(tmp(c_at*trans(conj(c_at))), c_at*conj(trans(c_at)));
 dlog << LTRACE << "1.6.1";
 c_check_equal(tmp(trans(c_at*trans(conj(c_at)))), trans(c_at*conj(trans(c_at))));
 dlog << LTRACE << "1.6.2";
 c_check_equal(tmp((c_sqr)*trans(conj(c_sqr))), (c_sqr)*conj(trans(c_sqr)));
 dlog << LTRACE << "1.6.2.1";
 c_check_equal(tmp(trans(c_sqr)*trans(conj(c_sqr))), trans(c_sqr)*conj(trans(c_sqr)));
 dlog << LTRACE << "1.6.3";
 c_check_equal(tmp(trans(trans(c_sqr)*trans(conj(c_sqr)))), trans(trans(c_sqr)*conj(trans(c_sqr))));
 dlog << LTRACE << "1.7";
 c_check_equal(tmp(conj(trans(c_at))*trans(conj(c_a))), conj(trans(c_at))*trans(conj(c_a)));
 c_check_equal(tmp(trans(conj(trans(c_at))*trans(conj(c_a)))), trans(conj(trans(c_at))*trans(conj(c_a))));
 dlog << LTRACE << "1.8";
 check_equal(tmp(a*trans(rowm(a,1))) , a*trans(rowm(a,1)));
 check_equal(tmp(a*colm(at,1)) , a*colm(at,1));
 check_equal(tmp(subm(a,1,1,2,2)*subm(a,1,2,2,2)), subm(a,1,1,2,2)*subm(a,1,2,2,2));
 dlog << LTRACE << "1.9";
 check_equal(tmp(trans(a*trans(rowm(a,1)))) , trans(a*trans(rowm(a,1))));
 dlog << LTRACE << "1.10";
 check_equal(tmp(trans(a*colm(at,1))) , trans(a*colm(at,1)));
 dlog << LTRACE << "1.11";
 check_equal(tmp(trans(subm(a,1,1,2,2)*subm(a,1,2,2,2))), trans(subm(a,1,1,2,2)*subm(a,1,2,2,2)));
 dlog << LTRACE << "1.12";
 {
 temp = at*a;
 temp2 = temp;
 temp += 3.5*at*a;
 assign_no_blas(temp2, temp2 + 3.5*at*a);
 check_equal(temp, temp2);
 temp -= at*3.5*a;
 assign_no_blas(temp2, temp2 - at*3.5*a);
 check_equal(temp, temp2);
 temp = temp + 4*at*a;
 assign_no_blas(temp2, temp2 + 4*at*a);
 check_equal(temp, temp2);
 temp = temp - 2.4*at*a;
 assign_no_blas(temp2, temp2 - 2.4*at*a);
 check_equal(temp, temp2);
 }
 dlog << LTRACE << "1.13";
 {
 temp = trans(at*a);
 temp2 = temp;
 temp3 = temp;
 dlog << LTRACE << "1.14";
 temp += trans(3.5*at*a);
 assign_no_blas(temp2, temp2 + trans(3.5*at*a));
 check_equal(temp, temp2);
 dlog << LTRACE << "1.15";
 temp -= trans(at*3.5*a);
 assign_no_blas(temp2, temp2 - trans(at*3.5*a));
 check_equal(temp, temp2);
 dlog << LTRACE << "1.16";
 temp = trans(temp + 4*at*a);
 assign_no_blas(temp3, trans(temp2 + 4*at*a));
 check_equal(temp, temp3);
 temp2 = temp;
 dlog << LTRACE << "1.17";
 temp = trans(temp - 2.4*at*a);
 assign_no_blas(temp3, trans(temp2 - 2.4*at*a));
 check_equal(temp, temp3);
 }
 dlog << LTRACE << "1.17.1";
 {
 matrix<type> m1, m2;
 m1 = matrix_cast<type>(randm(rows, cols, rnd));
 m2 = matrix_cast<type>(randm(cols, rows + 8, rnd));
 check_equal(tmp(m1*m2), m1*m2);
 check_equal(tmp(trans(m1*m2)), trans(m1*m2));
 m1 = trans(m1);
 check_equal(tmp(trans(m1)*m2), trans(m1)*m2);
 check_equal(tmp(trans(trans(m1)*m2)), trans(trans(m1)*m2));
 m2 = trans(m2);
 check_equal(tmp(trans(m1)*trans(m2)), trans(m1)*trans(m2));
 check_equal(tmp(trans(trans(m1)*trans(m2))), trans(trans(m1)*trans(m2)));
 m1 = trans(m1);
 check_equal(tmp(m1*trans(m2)), m1*trans(m2));
 check_equal(tmp(trans(m1*trans(m2))), trans(m1*trans(m2)));
 }
 dlog << LTRACE << "1.17.5";
 {
 matrix<type,1,0> r;
 matrix<type,0,1> c;
 r = matrix_cast<type>(randm(1, rows+9, rnd));
 c = matrix_cast<type>(randm(rows, 1, rnd));
 check_equal(tmp(c*r), c*r);
 check_equal(tmp(trans(c*r)), trans(c*r));
 check_equal(tmp(trans(r)*trans(c)), trans(r)*trans(c));
 check_equal(tmp(trans(trans(r)*trans(c))), trans(trans(r)*trans(c)));
 }
 dlog << LTRACE << "1.18";
 // GEMV tests
 check_equal(tmp(a*cv4), a*cv4);
 check_equal(tmp(trans(a*cv4)), trans(a*cv4));
 check_equal(tmp(rv3*a), rv3*a);
 check_equal(tmp(trans(cv4)*at), trans(cv4)*at);
 check_equal(tmp(a*trans(rv4)), a*trans(rv4));
 check_equal(tmp(trans(a*trans(rv4))), trans(a*trans(rv4)));
 check_equal(tmp(trans(a)*cv3), trans(a)*cv3);
 check_equal(tmp(rv4*trans(a)), rv4*trans(a));
 check_equal(tmp(trans(cv3)*trans(at)), trans(cv3)*trans(at));
 check_equal(tmp(trans(cv3)*a), trans(cv3)*a);
 check_equal(tmp(trans(a)*trans(rv3)), trans(a)*trans(rv3));
 c_check_equal(tmp(trans(conj(c_a))*c_cv3), trans(conj(c_a))*c_cv3);
 c_check_equal(tmp(c_rv4*trans(conj(c_a))), c_rv4*trans(conj(c_a)));
 c_check_equal(tmp(trans(c_cv3)*trans(conj(c_at))), trans(c_cv3)*trans(conj(c_at)));
 c_check_equal(tmp(conj(trans(c_a))*trans(c_rv3)), trans(conj(c_a))*trans(c_rv3));
 c_check_equal(tmp(c_rv4*conj(c_at)), c_rv4*conj(c_at));
 c_check_equal(tmp(trans(c_cv4)*conj(c_at)), trans(c_cv4)*conj(c_at));
 dlog << LTRACE << "2.00";
 c_check_equal(tmp(trans(trans(conj(c_a))*c_cv3)), trans(trans(conj(c_a))*c_cv3));
 c_check_equal(tmp(trans(c_rv4*trans(conj(c_a)))), trans(c_rv4*trans(conj(c_a))));
 c_check_equal(tmp(trans(trans(c_cv3)*trans(conj(c_at)))), trans(trans(c_cv3)*trans(conj(c_at))));
 dlog << LTRACE << "2.20";
 c_check_equal(tmp(trans(conj(trans(c_a))*trans(c_rv3))), trans(trans(conj(c_a))*trans(c_rv3)));
 c_check_equal(tmp(trans(c_rv4*conj(c_at))), trans(c_rv4*conj(c_at)));
 c_check_equal(tmp(trans(trans(c_cv4)*conj(c_at))), trans(trans(c_cv4)*conj(c_at)));
 dlog << LTRACE << "6";
 temp = a*at;
 check_equal(temp, a*at);
 temp = temp + a*at + trans(at)*at + trans(at)*sin(at);
 check_equal(temp, a*at + a*at+ trans(at)*at + trans(at)*sin(at));
 dlog << LTRACE << "6.1";
 temp = a*at;
 check_equal(temp, a*at);
 temp = a*at + temp;
 check_equal(temp, a*at + a*at);
 print_spinner();
 dlog << LTRACE << "6.2";
 temp = a*at;
 check_equal(temp, a*at);
 dlog << LTRACE << "6.2.3";
 temp = temp - a*at;
 dlog << LTRACE << "6.2.4";
 check_equal(temp, a*at-a*at);
 dlog << LTRACE << "6.3";
 temp = a*at;
 dlog << LTRACE << "6.3.5";
 check_equal(temp, a*at);
 dlog << LTRACE << "6.3.6";
 temp = a*at - temp;
 dlog << LTRACE << "6.4";
 check_equal(temp, a*at-a*at);
 const long d = min(rows,cols);
 rectangle rect(1,1,d,d);
 temp.set_size(max(rows,cols)+4,max(rows,cols)+4);
 set_all_elements(temp,4);
 temp2 = temp;
 dlog << LTRACE << "7";
 set_subm(temp,rect) = a*at;
 assign_no_blas( set_subm(temp2,rect) , a*at);
 check_equal(temp, temp2);
 temp = a;
 temp2 = a;
 set_colm(temp,1) = a*cv4;
 assign_no_blas( set_colm(temp2,1) , a*cv4);
 check_equal(temp, temp2);
 set_rowm(temp,1) = rv3*a;
 assign_no_blas( set_rowm(temp2,1) , rv3*a);
 check_equal(temp, temp2);
 // Test BLAS GER
 {
 temp.set_size(cols,cols);
 set_all_elements(temp,3);
 temp2 = temp;
 dlog << LTRACE << "8";
 temp += cv4*rv4;
 assign_no_blas(temp2, temp2 + cv4*rv4);
 check_equal(temp, temp2);
 dlog << LTRACE << "8.3";
 temp = temp + cv4*rv4;
 assign_no_blas(temp2, temp2 + cv4*rv4);
 check_equal(temp, temp2);
 dlog << LTRACE << "8.9";
 }
 {
 temp.set_size(cols,cols);
 set_all_elements(temp,3);
 temp2 = temp;
 temp3 = 0;
 dlog << LTRACE << "8.10";
 temp += trans(cv4*rv4);
 assign_no_blas(temp3, temp2 + trans(cv4*rv4));
 check_equal(temp, temp3);
 temp3 = 0;
 dlog << LTRACE << "8.11";
 temp2 = temp;
 temp = trans(temp + cv4*rv4);
 assign_no_blas(temp3, trans(temp2 + cv4*rv4));
 check_equal(temp, temp3);
 dlog << LTRACE << "8.12";
 }
 {
 matrix<complex<type> > temp, temp2, temp3;
 matrix<complex<type>,0,1 > cv4;
 matrix<complex<type>,1,0 > rv4;
 cv4.set_size(cols);
 rv4.set_size(cols);
 temp.set_size(cols,cols);
 set_all_elements(temp,complex<type>(3,5));
 temp(cols-1, cols-4) = 9;
 temp2 = temp;
 temp3.set_size(cols,cols);
 temp3 = 0;
 for (long i = 0; i < rv4.size(); ++i)
 {
 rv4(i) = complex<type>(rnd_num<type>(rnd),rnd_num<type>(rnd));
 cv4(i) = complex<type>(rnd_num<type>(rnd),rnd_num<type>(rnd));
 }
 dlog << LTRACE << "8.13";
 temp += trans(cv4*rv4);
 assign_no_blas(temp3, temp2 + trans(cv4*rv4));
 c_check_equal(temp, temp3);
 temp3 = 0;
 dlog << LTRACE << "8.14";
 temp2 = temp;
 temp = trans(temp + cv4*rv4);
 assign_no_blas(temp3, trans(temp2 + cv4*rv4));
 c_check_equal(temp, temp3);
 dlog << LTRACE << "8.15";
 }
 set_all_elements(c_temp, one + num1*i);
 c_temp2 = c_temp;
 set_all_elements(c_rv4, one + num2*i);
 set_all_elements(c_cv4, two + num3*i);
 dlog << LTRACE << "9";
 c_temp += c_cv4*c_rv4;
 assign_no_blas(c_temp2, c_temp2 + c_cv4*c_rv4);
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "9.1";
 c_temp += c_cv4*conj(c_rv4);
 assign_no_blas(c_temp2, c_temp2 + c_cv4*conj(c_rv4));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "9.2";
 c_temp = c_cv4*conj(c_rv4) + c_temp;
 assign_no_blas(c_temp2, c_temp2 + c_cv4*conj(c_rv4));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "9.3";
 c_temp = trans(c_rv4)*trans(conj(c_cv4)) + c_temp;
 assign_no_blas(c_temp2, c_temp2 + trans(c_rv4)*trans(conj(c_cv4)));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "9.4";
 c_temp += conj(c_cv4)*c_rv4;
 assign_no_blas(c_temp2, c_temp2 + conj(c_cv4)*c_rv4);
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "9.5";
 c_temp += conj(c_cv4)*conj(c_rv4);
 assign_no_blas(c_temp2, c_temp2 + conj(c_cv4)*conj(c_rv4));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "9.6";
 c_temp = conj(c_cv4)*conj(c_rv4) + c_temp;
 assign_no_blas(c_temp2, c_temp2 + conj(c_cv4)*conj(c_rv4));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "9.7";
 c_temp = conj(trans(c_rv4))*trans(conj(c_cv4)) + c_temp;
 assign_no_blas(c_temp2, c_temp2 + conj(trans(c_rv4))*trans(conj(c_cv4)));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10";
 c_temp += trans(c_cv4*c_rv4);
 assign_no_blas(c_temp2, c_temp2 + trans(c_cv4*c_rv4));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10.1";
 c_temp += trans(c_cv4*conj(c_rv4));
 assign_no_blas(c_temp2, c_temp2 + trans(c_cv4*conj(c_rv4)));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10.2";
 c_temp = trans(c_cv4*conj(c_rv4)) + c_temp;
 assign_no_blas(c_temp2, c_temp2 + trans(c_cv4*conj(c_rv4)));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10.3";
 c_temp = trans(trans(c_rv4)*trans(conj(c_cv4))) + c_temp;
 assign_no_blas(c_temp2, c_temp2 + trans(trans(c_rv4)*trans(conj(c_cv4))));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10.4";
 c_temp += trans(conj(c_cv4)*c_rv4);
 assign_no_blas(c_temp2, c_temp2 + trans(conj(c_cv4)*c_rv4));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10.5";
 c_temp += trans(conj(c_cv4)*conj(c_rv4));
 assign_no_blas(c_temp2, c_temp2 + trans(conj(c_cv4)*conj(c_rv4)));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10.6";
 c_temp = trans(conj(c_cv4)*conj(c_rv4)) + c_temp;
 assign_no_blas(c_temp2, c_temp2 + trans(conj(c_cv4)*conj(c_rv4)));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10.7";
 c_temp = trans(conj(trans(c_rv4))*trans(conj(c_cv4))) + c_temp;
 assign_no_blas(c_temp2, c_temp2 + trans(conj(trans(c_rv4))*trans(conj(c_cv4))));
 c_check_equal(c_temp, c_temp2);
 dlog << LTRACE << "10.8";
 print_spinner();
 // Test DOT
 check_equal( tmp(rv4*cv4), rv4*cv4);
 check_equal( tmp(trans(rv4*cv4)), trans(rv4*cv4));
 check_equal( tmp(trans(cv4)*trans(rv4)), trans(cv4)*trans(rv4));
 check_equal( tmp(rv4*3.9*cv4), rv4*3.9*cv4);
 check_equal( tmp(trans(cv4)*3.9*trans(rv4)), trans(cv4)*3.9*trans(rv4));
 check_equal( tmp(rv4*cv4*3.9), rv4*3.9*cv4);
 check_equal( tmp(trans(cv4)*trans(rv4)*3.9), trans(cv4)*3.9*trans(rv4));
 check_equal( tmp(trans(rv4*cv4)), trans(rv4*cv4));
 check_equal( tmp(trans(trans(rv4*cv4))), trans(trans(rv4*cv4)));
 check_equal( tmp(trans(trans(cv4)*trans(rv4))), trans(trans(cv4)*trans(rv4)));
 check_equal( tmp(trans(rv4*3.9*cv4)), trans(rv4*3.9*cv4));
 check_equal( tmp(trans(trans(cv4)*3.9*trans(rv4))), trans(trans(cv4)*3.9*trans(rv4)));
 check_equal( tmp(trans(rv4*cv4*3.9)), trans(rv4*3.9*cv4));
 check_equal( tmp(trans(trans(cv4)*trans(rv4)*3.9)), trans(trans(cv4)*3.9*trans(rv4)));
 temp.set_size(1,1);
 temp = 4;
 check_equal( tmp(temp + rv4*cv4), temp + rv4*cv4);
 check_equal( tmp(temp + trans(cv4)*trans(rv4)), temp + trans(cv4)*trans(rv4));
 dlog << LTRACE << "11";
 c_check_equal( tmp(conj(c_rv4)*c_cv4), conj(c_rv4)*c_cv4);
 c_check_equal( tmp(conj(trans(c_cv4))*trans(c_rv4)), trans(conj(c_cv4))*trans(c_rv4));
 c_check_equal( tmp(conj(c_rv4)*i*c_cv4), conj(c_rv4)*i*c_cv4);
 c_check_equal( tmp(conj(trans(c_cv4))*i*trans(c_rv4)), trans(conj(c_cv4))*i*trans(c_rv4));
 c_temp.set_size(1,1);
 c_temp = 4;
 c_check_equal( tmp(c_temp + conj(c_rv4)*c_cv4), c_temp + conj(c_rv4)*c_cv4);
 c_check_equal( tmp(c_temp + trans(conj(c_cv4))*trans(c_rv4)), c_temp + trans(conj(c_cv4))*trans(c_rv4));
 complex<type> tmp = c_rv4*c_cv4;
 DLIB_TEST(abs((tmp + i) - ((c_rv4*c_cv4)(0) + i)) < std::sqrt(std::numeric_limits<type>::epsilon())*eps_mul );
 DLIB_TEST(max(abs((rv4*cv4 + 1.0) - ((rv4*cv4)(0) + 1.0))) < std::sqrt(std::numeric_limits<type>::epsilon())*eps_mul);
 }
 {
 matrix<int> m(2,3), m2(6,1);
 m = 1,2,3,
 4,5,6;
 m2 = 1,2,3,4,5,6;
 DLIB_TEST(reshape_to_column_vector(m) == m2);
 DLIB_TEST(reshape_to_column_vector(m+m) == m2+m2);
 }
 {
 matrix<int,2,3> m(2,3);
 matrix<int> m2(6,1);
 m = 1,2,3,
 4,5,6;
 m2 = 1,2,3,4,5,6;
 DLIB_TEST(reshape_to_column_vector(m) == m2);
 DLIB_TEST(reshape_to_column_vector(m+m) == m2+m2);
 }
 }
 void matrix_test (
 )
 /*!
 ensures
 - runs tests on the matrix stuff compliance with the specs
 !*/
 { 
 print_spinner();
 {
 matrix<long> m1(2,2), m2(2,2);
 m1 = 1, 2,
 3, 4;
 m2 = 4, 5,
 6, 7;
 DLIB_TEST(subm(tensor_product(m1,m2),range(0,1), range(0,1)) == 1*m2);
 DLIB_TEST(subm(tensor_product(m1,m2),range(0,1), range(2,3)) == 2*m2);
 DLIB_TEST(subm(tensor_product(m1,m2),range(2,3), range(0,1)) == 3*m2);
 DLIB_TEST(subm(tensor_product(m1,m2),range(2,3), range(2,3)) == 4*m2);
 }
 {
 print_spinner();
 dlog << LTRACE << "testing blas stuff";
 dlog << LTRACE << " \nsmall double";
 test_blas<double>(3,4);
 print_spinner();
 dlog << LTRACE << " \nsmall float";
 test_blas<float>(3,4);
 print_spinner();
 dlog << LTRACE << " \nbig double";
 test_blas<double>(120,131);
 print_spinner();
 dlog << LTRACE << " \nbig float";
 test_blas<float>(120,131);
 print_spinner();
 dlog << LTRACE << "testing done";
 }
 {
 matrix<long> m(3,4), ml(3,4), mu(3,4);
 m = 1,2,3,4,
 4,5,6,7,
 7,8,9,0;
 ml = 1,0,0,0,
 4,5,0,0,
 7,8,9,0;
 mu = 1,2,3,4,
 0,5,6,7,
 0,0,9,0;
 DLIB_TEST(lowerm(m) == ml);
 DLIB_TEST(upperm(m) == mu);
 ml = 3,0,0,0,
 4,3,0,0,
 7,8,3,0;
 mu = 4,2,3,4,
 0,4,6,7,
 0,0,4,0;
 DLIB_TEST(lowerm(m,3) == ml);
 DLIB_TEST(upperm(m,4) == mu);
 }
 {
 matrix<long> m(3,4), row(1,3), col(2,1);
 m = 1,2,3,4,
 4,5,6,7,
 7,8,9,0;
 row = 4,5,6;
 col = 3,6;
 DLIB_TEST(rowm(m, 1, 3) == row);
 DLIB_TEST(colm(m, 2, 2) == col);
 }
 {
 std::vector<double> v(34, 8);
 std::vector<double> v2(34, 9);
 DLIB_TEST(mat(&v[0], v.size()) == mat(v));
 DLIB_TEST(mat(&v2[0], v.size()) != mat(v));
 }
 {
 std::vector<long> v(1, 3);
 std::vector<long> v2(1, 2);
 DLIB_TEST(mat(&v[0], v.size()) == mat(v));
 DLIB_TEST(mat(&v2[0], v.size()) != mat(v));
 }
 {
 matrix<double> a(3,3), b(3,3);
 a = 1, 2.5, 1,
 3, 4, 5,
 0.5, 2.2, 3;
 b = 0, 1, 0,
 1, 1, 1,
 0, 1, 1;
 DLIB_TEST((a>1) == b);
 DLIB_TEST((1<a) == b);
 b = 1, 1, 1,
 1, 1, 1,
 0, 1, 1;
 DLIB_TEST((a>=1) == b);
 DLIB_TEST((1<=a) == b);
 b = 0, 0, 0,
 0, 0, 0,
 0, 1, 0;
 DLIB_TEST((a==2.2) == b);
 DLIB_TEST((a!=2.2) == (b==0));
 DLIB_TEST((2.2==a) == b);
 DLIB_TEST((2.2!=a) == (0==b));
 b = 0, 0, 0,
 0, 0, 0,
 1, 0, 0;
 DLIB_TEST((a<1) == b);
 DLIB_TEST((1>a) == b);
 b = 1, 0, 1,
 0, 0, 0,
 1, 0, 0;
 DLIB_TEST((a<=1) == b);
 DLIB_TEST((1>=a) == b);
 }
 {
 matrix<double> a, b, c;
 a = randm(4,2);
 b += a;
 c -= a;
 DLIB_TEST(equal(a, b));
 DLIB_TEST(equal(-a, c));
 b += a;
 c -= a;
 DLIB_TEST(equal(2*a, b));
 DLIB_TEST(equal(-2*a, c));
 b += a + a;
 c -= a + a;
 DLIB_TEST(equal(4*a, b));
 DLIB_TEST(equal(-4*a, c));
 b.set_size(0,0);
 c.set_size(0,0);
 b += a + a;
 c -= a + a;
 DLIB_TEST(equal(2*a, b));
 DLIB_TEST(equal(-2*a, c));
 }
 {
 matrix<int> a, b, c;
 a.set_size(2, 3);
 b.set_size(2, 6);
 c.set_size(4, 3);
 a = 1, 2, 3,
 4, 5, 6;
 b = 1, 2, 3, 1, 2, 3,
 4, 5, 6, 4, 5, 6;
 c = 1, 2, 3,
 4, 5, 6,
 1, 2, 3,
 4, 5, 6;
 DLIB_TEST(join_rows(a,a) == b);
 DLIB_TEST(join_rows(a,abs(a)) == b);
 DLIB_TEST(join_cols(trans(a), trans(a)) == trans(b));
 DLIB_TEST(join_cols(a,a) == c);
 DLIB_TEST(join_cols(a,abs(a)) == c);
 DLIB_TEST(join_rows(trans(a),trans(a)) == trans(c));
 }
 {
 matrix<int, 2, 3> a;
 matrix<int, 2, 6> b;
 matrix<int, 4, 3> c;
 a = 1, 2, 3,
 4, 5, 6;
 b = 1, 2, 3, 1, 2, 3,
 4, 5, 6, 4, 5, 6;
 c = 1, 2, 3,
 4, 5, 6,
 1, 2, 3,
 4, 5, 6;
 DLIB_TEST(join_rows(a,a) == b);
 DLIB_TEST(join_rows(a,abs(a)) == b);
 DLIB_TEST(join_cols(trans(a), trans(a)) == trans(b));
 DLIB_TEST(join_cols(a,a) == c);
 DLIB_TEST(join_cols(a,abs(a)) == c);
 DLIB_TEST(join_rows(trans(a),trans(a)) == trans(c));
 }
 {
 matrix<int, 2, 3> a;
 matrix<int> a2;
 matrix<int, 2, 6> b;
 matrix<int, 4, 3> c;
 a = 1, 2, 3,
 4, 5, 6;
 a2 = a;
 b = 1, 2, 3, 1, 2, 3,
 4, 5, 6, 4, 5, 6;
 c = 1, 2, 3,
 4, 5, 6,
 1, 2, 3,
 4, 5, 6;
 DLIB_TEST(join_rows(a,a2) == b);
 DLIB_TEST(join_rows(a2,a) == b);
 DLIB_TEST(join_cols(trans(a2), trans(a)) == trans(b));
 DLIB_TEST(join_cols(a2,a) == c);
 DLIB_TEST(join_cols(a,a2) == c);
 DLIB_TEST(join_rows(trans(a2),trans(a)) == trans(c));
 }
 {
 matrix<int> a, b;
 a.set_size(2,3);
 a = 1, 2, 3,
 4, 5, 6;
 b.set_size(3,2);
 b = 1, 2,
 3, 4,
 5, 6;
 DLIB_TEST(reshape(a, 3, 2) == b);
 b.set_size(2,3);
 b = 1, 4, 2,
 5, 3, 6;
 DLIB_TEST(reshape(trans(a), 2, 3) == b);
 }
 {
 matrix<int,2,3> a;
 matrix<int> b;
 a = 1, 2, 3,
 4, 5, 6;
 b.set_size(3,2);
 b = 1, 2,
 3, 4,
 5, 6;
 DLIB_TEST(reshape(a, 3, 2) == b);
 b.set_size(2,3);
 b = 1, 4, 2,
 5, 3, 6;
 DLIB_TEST(reshape(trans(a), 2, 3) == b);
 }
 {
 std::vector<int> v(6);
 for (unsigned long i = 0; i < v.size(); ++i)
 v[i] = i;
 matrix<int,2,3> a;
 a = 0, 1, 2, 
 3, 4, 5;
 DLIB_TEST(mat(&v[0], 2, 3) == a);
 }
 {
 matrix<int> a(3,4);
 matrix<int> b(3,1), c(1,4);
 a = 1, 2, 3, 6,
 4, 5, 6, 9,
 1, 1, 1, 3;
 b(0) = sum(rowm(a,0));
 b(1) = sum(rowm(a,1));
 b(2) = sum(rowm(a,2));
 c(0) = sum(colm(a,0));
 c(1) = sum(colm(a,1));
 c(2) = sum(colm(a,2));
 c(3) = sum(colm(a,3));
 DLIB_TEST(sum_cols(a) == b);
 DLIB_TEST(sum_rows(a) == c);
 }
 {
 matrix<int> m(3,3);
 m = 1, 2, 3,
 4, 5, 6,
 7, 8, 9;
 DLIB_TEST(make_symmetric(m) == trans(make_symmetric(m)));
 DLIB_TEST(lowerm(make_symmetric(m)) == lowerm(m));
 DLIB_TEST(upperm(make_symmetric(m)) == trans(lowerm(m)));
 }
 {
 matrix<int,3,4> a;
 matrix<int> b(3,1), c(1,4);
 a = 1, 2, 3, 6,
 4, 5, 6, 9,
 1, 1, 1, 3;
 b(0) = sum(rowm(a,0));
 b(1) = sum(rowm(a,1));
 b(2) = sum(rowm(a,2));
 c(0) = sum(colm(a,0));
 c(1) = sum(colm(a,1));
 c(2) = sum(colm(a,2));
 c(3) = sum(colm(a,3));
 DLIB_TEST(sum_cols(a) == b);
 DLIB_TEST(sum_rows(a) == c);
 }
 {
 matrix<int> m(3,4), s(3,4);
 m = -2, 1, 5, -5,
 5, 5, 5, 5,
 9, 0, -4, -2;
 s = -1, 1, 1, -1,
 1, 1, 1, 1, 
 1, 1, -1, -1;
 DLIB_TEST(sign(m) == s);
 DLIB_TEST(sign(matrix_cast<double>(m)) == matrix_cast<double>(s));
 }
 }
 void test_matrix_IO()
 {
 dlib::rand rnd;
 print_spinner();
 for (int i = 0; i < 400; ++i)
 {
 ostringstream sout;
 sout.precision(20);
 matrix<double> m1, m2, m3;
 const long r = rnd.get_random_32bit_number()%7+1;
 const long c = rnd.get_random_32bit_number()%7+1;
 const long num = rnd.get_random_32bit_number()%2+1;
 m1 = randm(r,c,rnd);
 sout << m1;
 if (num != 1)
 sout << "\n" << m1;
 if (rnd.get_random_double() < 0.3)
 sout << " \n";
 else if (rnd.get_random_double() < 0.3)
 sout << " \n\n 3 3 3 3";
 else if (rnd.get_random_double() < 0.3)
 sout << " \n \n v 3 3 3 3 3";
 istringstream sin(sout.str());
 sin >> m2;
 DLIB_TEST_MSG(equal(m1,m2), m1 << "\n***********\n" << m2);
 if (num != 1)
 {
 sin >> m3;
 DLIB_TEST_MSG(equal(m1,m3), m1 << "\n***********\n" << m3);
 }
 }
 {
 istringstream sin(" 1 2\n3");
 matrix<double> m;
 DLIB_TEST(sin.good());
 sin >> m;
 DLIB_TEST(!sin.good());
 }
 {
 istringstream sin("");
 matrix<double> m;
 DLIB_TEST(sin.good());
 sin >> m;
 DLIB_TEST(!sin.good());
 }
 }
 void test_axpy()
 {
 const int n = 4;
 matrix<double> B = dlib::randm(n,n);
 matrix<double> g = dlib::uniform_matrix<double>(n,1,0.0);
 const double tau = 1;
 matrix<double> p = g + tau*dlib::colm(B,0);
 matrix<double> q = dlib::colm(B,0);
 DLIB_TEST(max(abs(p-q)) < 1e-14);
 p = tau*dlib::colm(B,0);
 q = dlib::colm(B,0);
 DLIB_TEST(max(abs(p-q)) < 1e-14);
 g = dlib::uniform_matrix<double>(n,n,0.0);
 p = g + tau*B;
 DLIB_TEST(max(abs(p-B)) < 1e-14);
 p = g + tau*subm(B,get_rect(B));
 DLIB_TEST(max(abs(p-B)) < 1e-14);
 g = dlib::uniform_matrix<double>(2,2,0.0);
 p = g + tau*subm(B,1,1,2,2);
 DLIB_TEST(max(abs(p-subm(B,1,1,2,2))) < 1e-14);
 set_subm(p,0,0,2,2) = g + tau*subm(B,1,1,2,2);
 DLIB_TEST(max(abs(p-subm(B,1,1,2,2))) < 1e-14);
 }
 class matrix_tester : public tester
 {
 public:
 matrix_tester (
 ) :
 tester ("test_matrix3",
 "Runs tests on the matrix component.")
 {}
 void perform_test (
 )
 {
 test_axpy();
 test_matrix_IO();
 matrix_test();
 }
 } a;
}