LCOV - code coverage report
Current view: top level - tests - catch_matrix.cpp (source / functions) Hit Total Coverage
Test: coverage.info Lines: 3235 3244 99.7 %
Date: 2022-02-13 10:56:14 Functions: 8 8 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : // Copyright (c) 2018-2022  Made to Order Software Corp.  All Rights Reserved
       2             : //
       3             : // https://snapwebsites.org/
       4             : // contact@m2osw.com
       5             : //
       6             : // This program is free software; you can redistribute it and/or modify
       7             : // it under the terms of the GNU General Public License as published by
       8             : // the Free Software Foundation; either version 2 of the License, or
       9             : // (at your option) any later version.
      10             : //
      11             : // This program is distributed in the hope that it will be useful,
      12             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      13             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      14             : // GNU General Public License for more details.
      15             : //
      16             : // You should have received a copy of the GNU General Public License along
      17             : // with this program; if not, write to the Free Software Foundation, Inc.,
      18             : // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
      19             : 
      20             : 
      21             : /** \file
      22             :  * \brief Verify the various matrix implementation.
      23             :  *
      24             :  * This file implements tests to verify that the matrix templates
      25             :  * work as expected.
      26             :  */
      27             : 
      28             : // self
      29             : //
      30             : #include    "catch_main.h"
      31             : 
      32             : 
      33             : 
      34             : // ignore the == and != against float warnings
      35             : //
      36             : #pragma GCC diagnostic ignored "-Wfloat-equal"
      37             : 
      38             : 
      39             : // snapdev lib
      40             : //
      41             : #include    <snapdev/matrix.h>
      42             : 
      43             : 
      44             : // last include
      45             : //
      46             : #include    <snapdev/poison.h>
      47             : 
      48             : 
      49             : 
      50             : 
      51             : 
      52             : namespace
      53             : {
      54         700 :     auto frand = []()
      55             :     {
      56         700 :         return static_cast<double>(rand()) / static_cast<double>(rand());
      57         700 :     };
      58             : }
      59             : // no name namespace
      60             : 
      61             : 
      62             : 
      63           6 : CATCH_TEST_CASE("matrix_init", "[matrix]")
      64             : {
      65             :     // constructor/copy
      66             :     // and
      67             :     // zero/identity
      68             :     //
      69           8 :     CATCH_GIVEN("constructor")
      70             :     {
      71           8 :         CATCH_START_SECTION("matrix: empty")
      72             :         {
      73           2 :             snapdev::matrix<double> empty;
      74             : 
      75           1 :             CATCH_REQUIRE(empty.rows() == 0);
      76           1 :             CATCH_REQUIRE(empty.columns() == 0);
      77             : 
      78           2 :             snapdev::matrix<double> copy(empty);
      79             : 
      80           1 :             CATCH_REQUIRE(empty.rows() == 0);
      81           1 :             CATCH_REQUIRE(empty.columns() == 0);
      82             :         }
      83             :         CATCH_END_SECTION()
      84             : 
      85           8 :         CATCH_START_SECTION("matrix: 2x2")
      86             :         {
      87           2 :             snapdev::matrix<double> m(2, 2);
      88             : 
      89           1 :             CATCH_REQUIRE(m.rows() == 2);
      90           1 :             CATCH_REQUIRE(m.columns() == 2);
      91             : 
      92             :             // by default we get an identity
      93             :             //
      94           1 :             CATCH_REQUIRE(m[0][0] == 1.0);
      95           1 :             CATCH_REQUIRE(m[0][1] == 0.0);
      96           1 :             CATCH_REQUIRE(m[1][0] == 0.0);
      97           1 :             CATCH_REQUIRE(m[1][1] == 1.0);
      98             : 
      99           1 :             double const r00(frand());
     100           1 :             double const r01(frand());
     101           1 :             double const r10(frand());
     102           1 :             double const r11(frand());
     103             : 
     104           1 :             m[0][0] = r00;
     105           1 :             m[0][1] = r01;
     106           1 :             m[1][0] = r10;
     107           1 :             m[1][1] = r11;
     108             : 
     109           1 :             CATCH_REQUIRE(m[0][0] == r00);
     110           1 :             CATCH_REQUIRE(m[0][1] == r01);
     111           1 :             CATCH_REQUIRE(m[1][0] == r10);
     112           1 :             CATCH_REQUIRE(m[1][1] == r11);
     113             : 
     114           2 :             snapdev::matrix<double> copy(m);
     115           2 :             snapdev::matrix<double> c2;
     116             : 
     117           1 :             CATCH_REQUIRE(copy[0][0] == r00);
     118           1 :             CATCH_REQUIRE(copy[0][1] == r01);
     119           1 :             CATCH_REQUIRE(copy[1][0] == r10);
     120           1 :             CATCH_REQUIRE(copy[1][1] == r11);
     121             : 
     122           1 :             m.clear();
     123             : 
     124           1 :             CATCH_REQUIRE(m[0][0] == 0.0);
     125           1 :             CATCH_REQUIRE(m[0][1] == 0.0);
     126           1 :             CATCH_REQUIRE(m[1][0] == 0.0);
     127           1 :             CATCH_REQUIRE(m[1][1] == 0.0);
     128             : 
     129             :             // copy is still intact
     130             :             //
     131           1 :             CATCH_REQUIRE(copy[0][0] == r00);
     132           1 :             CATCH_REQUIRE(copy[0][1] == r01);
     133           1 :             CATCH_REQUIRE(copy[1][0] == r10);
     134           1 :             CATCH_REQUIRE(copy[1][1] == r11);
     135             : 
     136           1 :             CATCH_REQUIRE(c2.rows() == 0);
     137           1 :             CATCH_REQUIRE(c2.columns() == 0);
     138             : 
     139           1 :             c2 = copy;
     140             : 
     141           1 :             CATCH_REQUIRE(c2[0][0] == r00);
     142           1 :             CATCH_REQUIRE(c2[0][1] == r01);
     143           1 :             CATCH_REQUIRE(c2[1][0] == r10);
     144           1 :             CATCH_REQUIRE(c2[1][1] == r11);
     145             : 
     146           1 :             c2.swap(m);
     147             : 
     148           1 :             CATCH_REQUIRE(m[0][0] == r00);
     149           1 :             CATCH_REQUIRE(m[0][1] == r01);
     150           1 :             CATCH_REQUIRE(m[1][0] == r10);
     151           1 :             CATCH_REQUIRE(m[1][1] == r11);
     152             : 
     153           1 :             CATCH_REQUIRE(c2[0][0] == 0.0);
     154           1 :             CATCH_REQUIRE(c2[0][1] == 0.0);
     155           1 :             CATCH_REQUIRE(c2[1][0] == 0.0);
     156           1 :             CATCH_REQUIRE(c2[1][1] == 0.0);
     157             :         }
     158             :         CATCH_END_SECTION()
     159             : 
     160           8 :         CATCH_START_SECTION("matrix: 3x3")
     161             :         {
     162           2 :             snapdev::matrix<double> m(3, 3);
     163             : 
     164           1 :             CATCH_REQUIRE(m.rows() == 3);
     165           1 :             CATCH_REQUIRE(m.columns() == 3);
     166             : 
     167             :             // by default we get an identity
     168             :             //
     169           1 :             CATCH_REQUIRE(m[0][0] == 1.0);
     170           1 :             CATCH_REQUIRE(m[0][1] == 0.0);
     171           1 :             CATCH_REQUIRE(m[0][2] == 0.0);
     172           1 :             CATCH_REQUIRE(m[1][0] == 0.0);
     173           1 :             CATCH_REQUIRE(m[1][1] == 1.0);
     174           1 :             CATCH_REQUIRE(m[1][2] == 0.0);
     175           1 :             CATCH_REQUIRE(m[2][0] == 0.0);
     176           1 :             CATCH_REQUIRE(m[2][1] == 0.0);
     177           1 :             CATCH_REQUIRE(m[2][2] == 1.0);
     178             : 
     179           1 :             double const r00(frand());
     180           1 :             double const r01(frand());
     181           1 :             double const r02(frand());
     182           1 :             double const r10(frand());
     183           1 :             double const r11(frand());
     184           1 :             double const r12(frand());
     185           1 :             double const r20(frand());
     186           1 :             double const r21(frand());
     187           1 :             double const r22(frand());
     188             : 
     189           1 :             m[0][0] = r00;
     190           1 :             m[0][1] = r01;
     191           1 :             m[0][2] = r02;
     192           1 :             m[1][0] = r10;
     193           1 :             m[1][1] = r11;
     194           1 :             m[1][2] = r12;
     195           1 :             m[2][0] = r20;
     196           1 :             m[2][1] = r21;
     197           1 :             m[2][2] = r22;
     198             : 
     199           1 :             CATCH_REQUIRE(m[0][0] == r00);
     200           1 :             CATCH_REQUIRE(m[0][1] == r01);
     201           1 :             CATCH_REQUIRE(m[0][2] == r02);
     202           1 :             CATCH_REQUIRE(m[1][0] == r10);
     203           1 :             CATCH_REQUIRE(m[1][1] == r11);
     204           1 :             CATCH_REQUIRE(m[1][2] == r12);
     205           1 :             CATCH_REQUIRE(m[2][0] == r20);
     206           1 :             CATCH_REQUIRE(m[2][1] == r21);
     207           1 :             CATCH_REQUIRE(m[2][2] == r22);
     208             : 
     209           2 :             snapdev::matrix<double> copy(m);
     210           2 :             snapdev::matrix<double> c2;
     211             : 
     212           1 :             CATCH_REQUIRE(copy[0][0] == r00);
     213           1 :             CATCH_REQUIRE(copy[0][1] == r01);
     214           1 :             CATCH_REQUIRE(copy[0][2] == r02);
     215           1 :             CATCH_REQUIRE(copy[1][0] == r10);
     216           1 :             CATCH_REQUIRE(copy[1][1] == r11);
     217           1 :             CATCH_REQUIRE(copy[1][2] == r12);
     218           1 :             CATCH_REQUIRE(copy[2][0] == r20);
     219           1 :             CATCH_REQUIRE(copy[2][1] == r21);
     220           1 :             CATCH_REQUIRE(copy[2][2] == r22);
     221             : 
     222           1 :             m.clear();
     223             : 
     224           1 :             CATCH_REQUIRE(m[0][0] == 0.0);
     225           1 :             CATCH_REQUIRE(m[0][1] == 0.0);
     226           1 :             CATCH_REQUIRE(m[0][2] == 0.0);
     227           1 :             CATCH_REQUIRE(m[1][0] == 0.0);
     228           1 :             CATCH_REQUIRE(m[1][1] == 0.0);
     229           1 :             CATCH_REQUIRE(m[1][2] == 0.0);
     230           1 :             CATCH_REQUIRE(m[2][0] == 0.0);
     231           1 :             CATCH_REQUIRE(m[2][1] == 0.0);
     232           1 :             CATCH_REQUIRE(m[2][2] == 0.0);
     233             : 
     234             :             // copy is still intact
     235             :             //
     236           1 :             CATCH_REQUIRE(copy[0][0] == r00);
     237           1 :             CATCH_REQUIRE(copy[0][1] == r01);
     238           1 :             CATCH_REQUIRE(copy[0][2] == r02);
     239           1 :             CATCH_REQUIRE(copy[1][0] == r10);
     240           1 :             CATCH_REQUIRE(copy[1][1] == r11);
     241           1 :             CATCH_REQUIRE(copy[1][2] == r12);
     242           1 :             CATCH_REQUIRE(copy[2][0] == r20);
     243           1 :             CATCH_REQUIRE(copy[2][1] == r21);
     244           1 :             CATCH_REQUIRE(copy[2][2] == r22);
     245             : 
     246           1 :             CATCH_REQUIRE(c2.rows() == 0);
     247           1 :             CATCH_REQUIRE(c2.columns() == 0);
     248             : 
     249           1 :             c2 = copy;
     250             : 
     251           1 :             CATCH_REQUIRE(c2[0][0] == r00);
     252           1 :             CATCH_REQUIRE(c2[0][1] == r01);
     253           1 :             CATCH_REQUIRE(c2[0][2] == r02);
     254           1 :             CATCH_REQUIRE(c2[1][0] == r10);
     255           1 :             CATCH_REQUIRE(c2[1][1] == r11);
     256           1 :             CATCH_REQUIRE(c2[1][2] == r12);
     257           1 :             CATCH_REQUIRE(c2[2][0] == r20);
     258           1 :             CATCH_REQUIRE(c2[2][1] == r21);
     259           1 :             CATCH_REQUIRE(c2[2][2] == r22);
     260             : 
     261           1 :             c2.swap(m);
     262             : 
     263           1 :             CATCH_REQUIRE(m[0][0] == r00);
     264           1 :             CATCH_REQUIRE(m[0][1] == r01);
     265           1 :             CATCH_REQUIRE(m[0][2] == r02);
     266           1 :             CATCH_REQUIRE(m[1][0] == r10);
     267           1 :             CATCH_REQUIRE(m[1][1] == r11);
     268           1 :             CATCH_REQUIRE(m[1][2] == r12);
     269           1 :             CATCH_REQUIRE(m[2][0] == r20);
     270           1 :             CATCH_REQUIRE(m[2][1] == r21);
     271           1 :             CATCH_REQUIRE(m[2][2] == r22);
     272             : 
     273           1 :             CATCH_REQUIRE(c2[0][0] == 0.0);
     274           1 :             CATCH_REQUIRE(c2[0][1] == 0.0);
     275           1 :             CATCH_REQUIRE(c2[0][2] == 0.0);
     276           1 :             CATCH_REQUIRE(c2[1][0] == 0.0);
     277           1 :             CATCH_REQUIRE(c2[1][1] == 0.0);
     278           1 :             CATCH_REQUIRE(c2[1][2] == 0.0);
     279           1 :             CATCH_REQUIRE(c2[2][0] == 0.0);
     280           1 :             CATCH_REQUIRE(c2[2][1] == 0.0);
     281           1 :             CATCH_REQUIRE(c2[2][2] == 0.0);
     282             :         }
     283             :         CATCH_END_SECTION()
     284             : 
     285           8 :         CATCH_START_SECTION("matrix: 4x4")
     286             :         {
     287           2 :             snapdev::matrix<double> m(4, 4);
     288             : 
     289           1 :             CATCH_REQUIRE(m.rows() == 4);
     290           1 :             CATCH_REQUIRE(m.columns() == 4);
     291             : 
     292             :             // by default we get an identity
     293             :             //
     294           1 :             CATCH_REQUIRE(m[0][0] == 1.0);
     295           1 :             CATCH_REQUIRE(m[0][1] == 0.0);
     296           1 :             CATCH_REQUIRE(m[0][2] == 0.0);
     297           1 :             CATCH_REQUIRE(m[0][3] == 0.0);
     298           1 :             CATCH_REQUIRE(m[1][0] == 0.0);
     299           1 :             CATCH_REQUIRE(m[1][1] == 1.0);
     300           1 :             CATCH_REQUIRE(m[1][2] == 0.0);
     301           1 :             CATCH_REQUIRE(m[1][3] == 0.0);
     302           1 :             CATCH_REQUIRE(m[2][0] == 0.0);
     303           1 :             CATCH_REQUIRE(m[2][1] == 0.0);
     304           1 :             CATCH_REQUIRE(m[2][2] == 1.0);
     305           1 :             CATCH_REQUIRE(m[2][3] == 0.0);
     306           1 :             CATCH_REQUIRE(m[3][0] == 0.0);
     307           1 :             CATCH_REQUIRE(m[3][1] == 0.0);
     308           1 :             CATCH_REQUIRE(m[3][2] == 0.0);
     309           1 :             CATCH_REQUIRE(m[3][3] == 1.0);
     310             : 
     311           1 :             double const r00(frand());
     312           1 :             double const r01(frand());
     313           1 :             double const r02(frand());
     314           1 :             double const r03(frand());
     315           1 :             double const r10(frand());
     316           1 :             double const r11(frand());
     317           1 :             double const r12(frand());
     318           1 :             double const r13(frand());
     319           1 :             double const r20(frand());
     320           1 :             double const r21(frand());
     321           1 :             double const r22(frand());
     322           1 :             double const r23(frand());
     323           1 :             double const r30(frand());
     324           1 :             double const r31(frand());
     325           1 :             double const r32(frand());
     326           1 :             double const r33(frand());
     327             : 
     328           1 :             m[0][0] = r00;
     329           1 :             m[0][1] = r01;
     330           1 :             m[0][2] = r02;
     331           1 :             m[0][3] = r03;
     332           1 :             m[1][0] = r10;
     333           1 :             m[1][1] = r11;
     334           1 :             m[1][2] = r12;
     335           1 :             m[1][3] = r13;
     336           1 :             m[2][0] = r20;
     337           1 :             m[2][1] = r21;
     338           1 :             m[2][2] = r22;
     339           1 :             m[2][3] = r23;
     340           1 :             m[3][0] = r30;
     341           1 :             m[3][1] = r31;
     342           1 :             m[3][2] = r32;
     343           1 :             m[3][3] = r33;
     344             : 
     345           1 :             CATCH_REQUIRE(m[0][0] == r00);
     346           1 :             CATCH_REQUIRE(m[0][1] == r01);
     347           1 :             CATCH_REQUIRE(m[0][2] == r02);
     348           1 :             CATCH_REQUIRE(m[0][3] == r03);
     349           1 :             CATCH_REQUIRE(m[1][0] == r10);
     350           1 :             CATCH_REQUIRE(m[1][1] == r11);
     351           1 :             CATCH_REQUIRE(m[1][2] == r12);
     352           1 :             CATCH_REQUIRE(m[1][3] == r13);
     353           1 :             CATCH_REQUIRE(m[2][0] == r20);
     354           1 :             CATCH_REQUIRE(m[2][1] == r21);
     355           1 :             CATCH_REQUIRE(m[2][2] == r22);
     356           1 :             CATCH_REQUIRE(m[2][3] == r23);
     357           1 :             CATCH_REQUIRE(m[3][0] == r30);
     358           1 :             CATCH_REQUIRE(m[3][1] == r31);
     359           1 :             CATCH_REQUIRE(m[3][2] == r32);
     360           1 :             CATCH_REQUIRE(m[3][3] == r33);
     361             : 
     362           2 :             snapdev::matrix<double> copy(m);
     363           2 :             snapdev::matrix<double> c2;
     364             : 
     365           1 :             CATCH_REQUIRE(copy[0][0] == r00);
     366           1 :             CATCH_REQUIRE(copy[0][1] == r01);
     367           1 :             CATCH_REQUIRE(copy[0][2] == r02);
     368           1 :             CATCH_REQUIRE(copy[0][3] == r03);
     369           1 :             CATCH_REQUIRE(copy[1][0] == r10);
     370           1 :             CATCH_REQUIRE(copy[1][1] == r11);
     371           1 :             CATCH_REQUIRE(copy[1][2] == r12);
     372           1 :             CATCH_REQUIRE(copy[1][3] == r13);
     373           1 :             CATCH_REQUIRE(copy[2][0] == r20);
     374           1 :             CATCH_REQUIRE(copy[2][1] == r21);
     375           1 :             CATCH_REQUIRE(copy[2][2] == r22);
     376           1 :             CATCH_REQUIRE(copy[2][3] == r23);
     377           1 :             CATCH_REQUIRE(copy[3][0] == r30);
     378           1 :             CATCH_REQUIRE(copy[3][1] == r31);
     379           1 :             CATCH_REQUIRE(copy[3][2] == r32);
     380           1 :             CATCH_REQUIRE(copy[3][3] == r33);
     381             : 
     382           1 :             m.clear();
     383             : 
     384           1 :             CATCH_REQUIRE(m[0][0] == 0.0);
     385           1 :             CATCH_REQUIRE(m[0][1] == 0.0);
     386           1 :             CATCH_REQUIRE(m[0][2] == 0.0);
     387           1 :             CATCH_REQUIRE(m[0][3] == 0.0);
     388           1 :             CATCH_REQUIRE(m[1][0] == 0.0);
     389           1 :             CATCH_REQUIRE(m[1][1] == 0.0);
     390           1 :             CATCH_REQUIRE(m[1][2] == 0.0);
     391           1 :             CATCH_REQUIRE(m[1][3] == 0.0);
     392           1 :             CATCH_REQUIRE(m[2][0] == 0.0);
     393           1 :             CATCH_REQUIRE(m[2][1] == 0.0);
     394           1 :             CATCH_REQUIRE(m[2][2] == 0.0);
     395           1 :             CATCH_REQUIRE(m[2][3] == 0.0);
     396             : 
     397             :             // copy is still intact
     398             :             //
     399           1 :             CATCH_REQUIRE(copy[0][0] == r00);
     400           1 :             CATCH_REQUIRE(copy[0][1] == r01);
     401           1 :             CATCH_REQUIRE(copy[0][2] == r02);
     402           1 :             CATCH_REQUIRE(copy[0][3] == r03);
     403           1 :             CATCH_REQUIRE(copy[1][0] == r10);
     404           1 :             CATCH_REQUIRE(copy[1][1] == r11);
     405           1 :             CATCH_REQUIRE(copy[1][2] == r12);
     406           1 :             CATCH_REQUIRE(copy[1][3] == r13);
     407           1 :             CATCH_REQUIRE(copy[2][0] == r20);
     408           1 :             CATCH_REQUIRE(copy[2][1] == r21);
     409           1 :             CATCH_REQUIRE(copy[2][2] == r22);
     410           1 :             CATCH_REQUIRE(copy[2][3] == r23);
     411           1 :             CATCH_REQUIRE(copy[3][0] == r30);
     412           1 :             CATCH_REQUIRE(copy[3][1] == r31);
     413           1 :             CATCH_REQUIRE(copy[3][2] == r32);
     414           1 :             CATCH_REQUIRE(copy[3][3] == r33);
     415             : 
     416           1 :             CATCH_REQUIRE(c2.rows() == 0);
     417           1 :             CATCH_REQUIRE(c2.columns() == 0);
     418             : 
     419           1 :             c2 = copy;
     420             : 
     421           1 :             CATCH_REQUIRE(c2[0][0] == r00);
     422           1 :             CATCH_REQUIRE(c2[0][1] == r01);
     423           1 :             CATCH_REQUIRE(c2[0][2] == r02);
     424           1 :             CATCH_REQUIRE(c2[0][3] == r03);
     425           1 :             CATCH_REQUIRE(c2[1][0] == r10);
     426           1 :             CATCH_REQUIRE(c2[1][1] == r11);
     427           1 :             CATCH_REQUIRE(c2[1][2] == r12);
     428           1 :             CATCH_REQUIRE(c2[1][3] == r13);
     429           1 :             CATCH_REQUIRE(c2[2][0] == r20);
     430           1 :             CATCH_REQUIRE(c2[2][1] == r21);
     431           1 :             CATCH_REQUIRE(c2[2][2] == r22);
     432           1 :             CATCH_REQUIRE(c2[2][3] == r23);
     433           1 :             CATCH_REQUIRE(c2[3][0] == r30);
     434           1 :             CATCH_REQUIRE(c2[3][1] == r31);
     435           1 :             CATCH_REQUIRE(c2[3][2] == r32);
     436           1 :             CATCH_REQUIRE(c2[3][3] == r33);
     437             : 
     438             :             //std::cout << c2 << std::endl;
     439             : 
     440           1 :             c2.swap(m);
     441             : 
     442           1 :             CATCH_REQUIRE(m[0][0] == r00);
     443           1 :             CATCH_REQUIRE(m[0][1] == r01);
     444           1 :             CATCH_REQUIRE(m[0][2] == r02);
     445           1 :             CATCH_REQUIRE(m[0][3] == r03);
     446           1 :             CATCH_REQUIRE(m[1][0] == r10);
     447           1 :             CATCH_REQUIRE(m[1][1] == r11);
     448           1 :             CATCH_REQUIRE(m[1][2] == r12);
     449           1 :             CATCH_REQUIRE(m[1][3] == r13);
     450           1 :             CATCH_REQUIRE(m[2][0] == r20);
     451           1 :             CATCH_REQUIRE(m[2][1] == r21);
     452           1 :             CATCH_REQUIRE(m[2][2] == r22);
     453           1 :             CATCH_REQUIRE(m[2][3] == r23);
     454           1 :             CATCH_REQUIRE(m[3][0] == r30);
     455           1 :             CATCH_REQUIRE(m[3][1] == r31);
     456           1 :             CATCH_REQUIRE(m[3][2] == r32);
     457           1 :             CATCH_REQUIRE(m[3][3] == r33);
     458             : 
     459           1 :             CATCH_REQUIRE(c2[0][0] == 0.0);
     460           1 :             CATCH_REQUIRE(c2[0][1] == 0.0);
     461           1 :             CATCH_REQUIRE(c2[0][2] == 0.0);
     462           1 :             CATCH_REQUIRE(c2[0][3] == 0.0);
     463           1 :             CATCH_REQUIRE(c2[1][0] == 0.0);
     464           1 :             CATCH_REQUIRE(c2[1][1] == 0.0);
     465           1 :             CATCH_REQUIRE(c2[1][2] == 0.0);
     466           1 :             CATCH_REQUIRE(c2[1][3] == 0.0);
     467           1 :             CATCH_REQUIRE(c2[2][0] == 0.0);
     468           1 :             CATCH_REQUIRE(c2[2][1] == 0.0);
     469           1 :             CATCH_REQUIRE(c2[2][2] == 0.0);
     470           1 :             CATCH_REQUIRE(c2[2][3] == 0.0);
     471             :         }
     472             :         CATCH_END_SECTION()
     473             :     }
     474           4 : }
     475             : 
     476             : 
     477          10 : CATCH_TEST_CASE("matrix_additive", "[matrix]")
     478             : {
     479             :     // create two random 4x4 matrices and make sure the add works
     480             :     //
     481          16 :     CATCH_GIVEN("add")
     482             :     {
     483           8 :         CATCH_START_SECTION("matrix: a+=<scalar>")
     484             :         {
     485             :             // setup A
     486             :             //
     487           2 :             snapdev::matrix<double> a(4, 4);
     488             : 
     489           1 :             CATCH_REQUIRE(a.rows() == 4);
     490           1 :             CATCH_REQUIRE(a.columns() == 4);
     491             : 
     492           1 :             double const a00(frand());
     493           1 :             double const a01(frand());
     494           1 :             double const a02(frand());
     495           1 :             double const a03(frand());
     496           1 :             double const a10(frand());
     497           1 :             double const a11(frand());
     498           1 :             double const a12(frand());
     499           1 :             double const a13(frand());
     500           1 :             double const a20(frand());
     501           1 :             double const a21(frand());
     502           1 :             double const a22(frand());
     503           1 :             double const a23(frand());
     504           1 :             double const a30(frand());
     505           1 :             double const a31(frand());
     506           1 :             double const a32(frand());
     507           1 :             double const a33(frand());
     508             : 
     509           1 :             a[0][0] = a00;
     510           1 :             a[0][1] = a01;
     511           1 :             a[0][2] = a02;
     512           1 :             a[0][3] = a03;
     513           1 :             a[1][0] = a10;
     514           1 :             a[1][1] = a11;
     515           1 :             a[1][2] = a12;
     516           1 :             a[1][3] = a13;
     517           1 :             a[2][0] = a20;
     518           1 :             a[2][1] = a21;
     519           1 :             a[2][2] = a22;
     520           1 :             a[2][3] = a23;
     521           1 :             a[3][0] = a30;
     522           1 :             a[3][1] = a31;
     523           1 :             a[3][2] = a32;
     524           1 :             a[3][3] = a33;
     525             : 
     526           1 :             CATCH_REQUIRE(a[0][0] == a00);
     527           1 :             CATCH_REQUIRE(a[0][1] == a01);
     528           1 :             CATCH_REQUIRE(a[0][2] == a02);
     529           1 :             CATCH_REQUIRE(a[0][3] == a03);
     530           1 :             CATCH_REQUIRE(a[1][0] == a10);
     531           1 :             CATCH_REQUIRE(a[1][1] == a11);
     532           1 :             CATCH_REQUIRE(a[1][2] == a12);
     533           1 :             CATCH_REQUIRE(a[1][3] == a13);
     534           1 :             CATCH_REQUIRE(a[2][0] == a20);
     535           1 :             CATCH_REQUIRE(a[2][1] == a21);
     536           1 :             CATCH_REQUIRE(a[2][2] == a22);
     537           1 :             CATCH_REQUIRE(a[2][3] == a23);
     538           1 :             CATCH_REQUIRE(a[3][0] == a30);
     539           1 :             CATCH_REQUIRE(a[3][1] == a31);
     540           1 :             CATCH_REQUIRE(a[3][2] == a32);
     541           1 :             CATCH_REQUIRE(a[3][3] == a33);
     542             : 
     543             :             // create a scalar for our test
     544             :             //
     545           1 :             double scalar(frand());
     546           1 :             while(scalar == 0.0)
     547             :             {
     548           0 :                 scalar = frand();
     549             :             }
     550             : 
     551             :             // run operation A += <scalar>
     552             :             //
     553           1 :             a += scalar;
     554             : 
     555             :             // this can't fail because we ensure scalar != 0
     556             :             //
     557           1 :             CATCH_REQUIRE_FALSE(a[0][0] == a00);
     558           1 :             CATCH_REQUIRE_FALSE(a[0][1] == a01);
     559           1 :             CATCH_REQUIRE_FALSE(a[0][2] == a02);
     560           1 :             CATCH_REQUIRE_FALSE(a[0][3] == a03);
     561           1 :             CATCH_REQUIRE_FALSE(a[1][0] == a10);
     562           1 :             CATCH_REQUIRE_FALSE(a[1][1] == a11);
     563           1 :             CATCH_REQUIRE_FALSE(a[1][2] == a12);
     564           1 :             CATCH_REQUIRE_FALSE(a[1][3] == a13);
     565           1 :             CATCH_REQUIRE_FALSE(a[2][0] == a20);
     566           1 :             CATCH_REQUIRE_FALSE(a[2][1] == a21);
     567           1 :             CATCH_REQUIRE_FALSE(a[2][2] == a22);
     568           1 :             CATCH_REQUIRE_FALSE(a[2][3] == a23);
     569           1 :             CATCH_REQUIRE_FALSE(a[3][0] == a30);
     570           1 :             CATCH_REQUIRE_FALSE(a[3][1] == a31);
     571           1 :             CATCH_REQUIRE_FALSE(a[3][2] == a32);
     572           1 :             CATCH_REQUIRE_FALSE(a[3][3] == a33);
     573             : 
     574           1 :             CATCH_REQUIRE(a[0][0] == a00 + scalar);
     575           1 :             CATCH_REQUIRE(a[0][1] == a01 + scalar);
     576           1 :             CATCH_REQUIRE(a[0][2] == a02 + scalar);
     577           1 :             CATCH_REQUIRE(a[0][3] == a03 + scalar);
     578           1 :             CATCH_REQUIRE(a[1][0] == a10 + scalar);
     579           1 :             CATCH_REQUIRE(a[1][1] == a11 + scalar);
     580           1 :             CATCH_REQUIRE(a[1][2] == a12 + scalar);
     581           1 :             CATCH_REQUIRE(a[1][3] == a13 + scalar);
     582           1 :             CATCH_REQUIRE(a[2][0] == a20 + scalar);
     583           1 :             CATCH_REQUIRE(a[2][1] == a21 + scalar);
     584           1 :             CATCH_REQUIRE(a[2][2] == a22 + scalar);
     585           1 :             CATCH_REQUIRE(a[2][3] == a23 + scalar);
     586           1 :             CATCH_REQUIRE(a[3][0] == a30 + scalar);
     587           1 :             CATCH_REQUIRE(a[3][1] == a31 + scalar);
     588           1 :             CATCH_REQUIRE(a[3][2] == a32 + scalar);
     589           1 :             CATCH_REQUIRE(a[3][3] == a33 + scalar);
     590             :         }
     591             :         CATCH_END_SECTION()
     592             : 
     593           8 :         CATCH_START_SECTION("matrix: b=a+<scalar>")
     594             :         {
     595             :             // setup A
     596             :             //
     597           2 :             snapdev::matrix<double> a(4, 4);
     598             : 
     599           1 :             CATCH_REQUIRE(a.rows() == 4);
     600           1 :             CATCH_REQUIRE(a.columns() == 4);
     601             : 
     602           1 :             double const a00(frand());
     603           1 :             double const a01(frand());
     604           1 :             double const a02(frand());
     605           1 :             double const a03(frand());
     606           1 :             double const a10(frand());
     607           1 :             double const a11(frand());
     608           1 :             double const a12(frand());
     609           1 :             double const a13(frand());
     610           1 :             double const a20(frand());
     611           1 :             double const a21(frand());
     612           1 :             double const a22(frand());
     613           1 :             double const a23(frand());
     614           1 :             double const a30(frand());
     615           1 :             double const a31(frand());
     616           1 :             double const a32(frand());
     617           1 :             double const a33(frand());
     618             : 
     619           1 :             a[0][0] = a00;
     620           1 :             a[0][1] = a01;
     621           1 :             a[0][2] = a02;
     622           1 :             a[0][3] = a03;
     623           1 :             a[1][0] = a10;
     624           1 :             a[1][1] = a11;
     625           1 :             a[1][2] = a12;
     626           1 :             a[1][3] = a13;
     627           1 :             a[2][0] = a20;
     628           1 :             a[2][1] = a21;
     629           1 :             a[2][2] = a22;
     630           1 :             a[2][3] = a23;
     631           1 :             a[3][0] = a30;
     632           1 :             a[3][1] = a31;
     633           1 :             a[3][2] = a32;
     634           1 :             a[3][3] = a33;
     635             : 
     636           1 :             CATCH_REQUIRE(a[0][0] == a00);
     637           1 :             CATCH_REQUIRE(a[0][1] == a01);
     638           1 :             CATCH_REQUIRE(a[0][2] == a02);
     639           1 :             CATCH_REQUIRE(a[0][3] == a03);
     640           1 :             CATCH_REQUIRE(a[1][0] == a10);
     641           1 :             CATCH_REQUIRE(a[1][1] == a11);
     642           1 :             CATCH_REQUIRE(a[1][2] == a12);
     643           1 :             CATCH_REQUIRE(a[1][3] == a13);
     644           1 :             CATCH_REQUIRE(a[2][0] == a20);
     645           1 :             CATCH_REQUIRE(a[2][1] == a21);
     646           1 :             CATCH_REQUIRE(a[2][2] == a22);
     647           1 :             CATCH_REQUIRE(a[2][3] == a23);
     648           1 :             CATCH_REQUIRE(a[3][0] == a30);
     649           1 :             CATCH_REQUIRE(a[3][1] == a31);
     650           1 :             CATCH_REQUIRE(a[3][2] == a32);
     651           1 :             CATCH_REQUIRE(a[3][3] == a33);
     652             : 
     653             :             // setup B
     654             :             //
     655           2 :             snapdev::matrix<double> b(4, 4);
     656             : 
     657           1 :             CATCH_REQUIRE(b.rows() == 4);
     658           1 :             CATCH_REQUIRE(b.columns() == 4);
     659             : 
     660           1 :             double const b00(frand());
     661           1 :             double const b01(frand());
     662           1 :             double const b02(frand());
     663           1 :             double const b03(frand());
     664           1 :             double const b10(frand());
     665           1 :             double const b11(frand());
     666           1 :             double const b12(frand());
     667           1 :             double const b13(frand());
     668           1 :             double const b20(frand());
     669           1 :             double const b21(frand());
     670           1 :             double const b22(frand());
     671           1 :             double const b23(frand());
     672           1 :             double const b30(frand());
     673           1 :             double const b31(frand());
     674           1 :             double const b32(frand());
     675           1 :             double const b33(frand());
     676             : 
     677           1 :             b[0][0] = b00;
     678           1 :             b[0][1] = b01;
     679           1 :             b[0][2] = b02;
     680           1 :             b[0][3] = b03;
     681           1 :             b[1][0] = b10;
     682           1 :             b[1][1] = b11;
     683           1 :             b[1][2] = b12;
     684           1 :             b[1][3] = b13;
     685           1 :             b[2][0] = b20;
     686           1 :             b[2][1] = b21;
     687           1 :             b[2][2] = b22;
     688           1 :             b[2][3] = b23;
     689           1 :             b[3][0] = b30;
     690           1 :             b[3][1] = b31;
     691           1 :             b[3][2] = b32;
     692           1 :             b[3][3] = b33;
     693             : 
     694           1 :             CATCH_REQUIRE(b[0][0] == b00);
     695           1 :             CATCH_REQUIRE(b[0][1] == b01);
     696           1 :             CATCH_REQUIRE(b[0][2] == b02);
     697           1 :             CATCH_REQUIRE(b[0][3] == b03);
     698           1 :             CATCH_REQUIRE(b[1][0] == b10);
     699           1 :             CATCH_REQUIRE(b[1][1] == b11);
     700           1 :             CATCH_REQUIRE(b[1][2] == b12);
     701           1 :             CATCH_REQUIRE(b[1][3] == b13);
     702           1 :             CATCH_REQUIRE(b[2][0] == b20);
     703           1 :             CATCH_REQUIRE(b[2][1] == b21);
     704           1 :             CATCH_REQUIRE(b[2][2] == b22);
     705           1 :             CATCH_REQUIRE(b[2][3] == b23);
     706           1 :             CATCH_REQUIRE(b[3][0] == b30);
     707           1 :             CATCH_REQUIRE(b[3][1] == b31);
     708           1 :             CATCH_REQUIRE(b[3][2] == b32);
     709           1 :             CATCH_REQUIRE(b[3][3] == b33);
     710             : 
     711             :             // create a scalar for our test
     712             :             //
     713           1 :             double scalar(frand());
     714           1 :             while(scalar == 0.0)
     715             :             {
     716           0 :                 scalar = frand();
     717             :             }
     718             : 
     719             :             // run operation B = A + <scalar>
     720             :             //
     721           1 :             b = a + scalar;
     722             : 
     723           1 :             CATCH_REQUIRE(a[0][0] == a00);
     724           1 :             CATCH_REQUIRE(a[0][1] == a01);
     725           1 :             CATCH_REQUIRE(a[0][2] == a02);
     726           1 :             CATCH_REQUIRE(a[0][3] == a03);
     727           1 :             CATCH_REQUIRE(a[1][0] == a10);
     728           1 :             CATCH_REQUIRE(a[1][1] == a11);
     729           1 :             CATCH_REQUIRE(a[1][2] == a12);
     730           1 :             CATCH_REQUIRE(a[1][3] == a13);
     731           1 :             CATCH_REQUIRE(a[2][0] == a20);
     732           1 :             CATCH_REQUIRE(a[2][1] == a21);
     733           1 :             CATCH_REQUIRE(a[2][2] == a22);
     734           1 :             CATCH_REQUIRE(a[2][3] == a23);
     735           1 :             CATCH_REQUIRE(a[3][0] == a30);
     736           1 :             CATCH_REQUIRE(a[3][1] == a31);
     737           1 :             CATCH_REQUIRE(a[3][2] == a32);
     738           1 :             CATCH_REQUIRE(a[3][3] == a33);
     739             : 
     740             :             // this can't fail because we ensure scalar != 0
     741             :             //
     742           1 :             CATCH_REQUIRE_FALSE(b[0][0] == b00);
     743           1 :             CATCH_REQUIRE_FALSE(b[0][1] == b01);
     744           1 :             CATCH_REQUIRE_FALSE(b[0][2] == b02);
     745           1 :             CATCH_REQUIRE_FALSE(b[0][3] == b03);
     746           1 :             CATCH_REQUIRE_FALSE(b[1][0] == b10);
     747           1 :             CATCH_REQUIRE_FALSE(b[1][1] == b11);
     748           1 :             CATCH_REQUIRE_FALSE(b[1][2] == b12);
     749           1 :             CATCH_REQUIRE_FALSE(b[1][3] == b13);
     750           1 :             CATCH_REQUIRE_FALSE(b[2][0] == b20);
     751           1 :             CATCH_REQUIRE_FALSE(b[2][1] == b21);
     752           1 :             CATCH_REQUIRE_FALSE(b[2][2] == b22);
     753           1 :             CATCH_REQUIRE_FALSE(b[2][3] == b23);
     754           1 :             CATCH_REQUIRE_FALSE(b[3][0] == b30);
     755           1 :             CATCH_REQUIRE_FALSE(b[3][1] == b31);
     756           1 :             CATCH_REQUIRE_FALSE(b[3][2] == b32);
     757           1 :             CATCH_REQUIRE_FALSE(b[3][3] == b33);
     758             : 
     759           1 :             CATCH_REQUIRE(b[0][0] == a00 + scalar);
     760           1 :             CATCH_REQUIRE(b[0][1] == a01 + scalar);
     761           1 :             CATCH_REQUIRE(b[0][2] == a02 + scalar);
     762           1 :             CATCH_REQUIRE(b[0][3] == a03 + scalar);
     763           1 :             CATCH_REQUIRE(b[1][0] == a10 + scalar);
     764           1 :             CATCH_REQUIRE(b[1][1] == a11 + scalar);
     765           1 :             CATCH_REQUIRE(b[1][2] == a12 + scalar);
     766           1 :             CATCH_REQUIRE(b[1][3] == a13 + scalar);
     767           1 :             CATCH_REQUIRE(b[2][0] == a20 + scalar);
     768           1 :             CATCH_REQUIRE(b[2][1] == a21 + scalar);
     769           1 :             CATCH_REQUIRE(b[2][2] == a22 + scalar);
     770           1 :             CATCH_REQUIRE(b[2][3] == a23 + scalar);
     771           1 :             CATCH_REQUIRE(b[3][0] == a30 + scalar);
     772           1 :             CATCH_REQUIRE(b[3][1] == a31 + scalar);
     773           1 :             CATCH_REQUIRE(b[3][2] == a32 + scalar);
     774           1 :             CATCH_REQUIRE(b[3][3] == a33 + scalar);
     775             :         }
     776             :         CATCH_END_SECTION()
     777             : 
     778           8 :         CATCH_START_SECTION("matrix: c=a+b")
     779             :         {
     780             :             // setup A
     781             :             //
     782           2 :             snapdev::matrix<double> a(4, 4);
     783             : 
     784           1 :             CATCH_REQUIRE(a.rows() == 4);
     785           1 :             CATCH_REQUIRE(a.columns() == 4);
     786             : 
     787           1 :             double const a00(frand());
     788           1 :             double const a01(frand());
     789           1 :             double const a02(frand());
     790           1 :             double const a03(frand());
     791           1 :             double const a10(frand());
     792           1 :             double const a11(frand());
     793           1 :             double const a12(frand());
     794           1 :             double const a13(frand());
     795           1 :             double const a20(frand());
     796           1 :             double const a21(frand());
     797           1 :             double const a22(frand());
     798           1 :             double const a23(frand());
     799           1 :             double const a30(frand());
     800           1 :             double const a31(frand());
     801           1 :             double const a32(frand());
     802           1 :             double const a33(frand());
     803             : 
     804           1 :             a[0][0] = a00;
     805           1 :             a[0][1] = a01;
     806           1 :             a[0][2] = a02;
     807           1 :             a[0][3] = a03;
     808           1 :             a[1][0] = a10;
     809           1 :             a[1][1] = a11;
     810           1 :             a[1][2] = a12;
     811           1 :             a[1][3] = a13;
     812           1 :             a[2][0] = a20;
     813           1 :             a[2][1] = a21;
     814           1 :             a[2][2] = a22;
     815           1 :             a[2][3] = a23;
     816           1 :             a[3][0] = a30;
     817           1 :             a[3][1] = a31;
     818           1 :             a[3][2] = a32;
     819           1 :             a[3][3] = a33;
     820             : 
     821           1 :             CATCH_REQUIRE(a[0][0] == a00);
     822           1 :             CATCH_REQUIRE(a[0][1] == a01);
     823           1 :             CATCH_REQUIRE(a[0][2] == a02);
     824           1 :             CATCH_REQUIRE(a[0][3] == a03);
     825           1 :             CATCH_REQUIRE(a[1][0] == a10);
     826           1 :             CATCH_REQUIRE(a[1][1] == a11);
     827           1 :             CATCH_REQUIRE(a[1][2] == a12);
     828           1 :             CATCH_REQUIRE(a[1][3] == a13);
     829           1 :             CATCH_REQUIRE(a[2][0] == a20);
     830           1 :             CATCH_REQUIRE(a[2][1] == a21);
     831           1 :             CATCH_REQUIRE(a[2][2] == a22);
     832           1 :             CATCH_REQUIRE(a[2][3] == a23);
     833           1 :             CATCH_REQUIRE(a[3][0] == a30);
     834           1 :             CATCH_REQUIRE(a[3][1] == a31);
     835           1 :             CATCH_REQUIRE(a[3][2] == a32);
     836           1 :             CATCH_REQUIRE(a[3][3] == a33);
     837             : 
     838             :             // setup B
     839             :             //
     840           2 :             snapdev::matrix<double> b(4, 4);
     841             : 
     842           1 :             CATCH_REQUIRE(b.rows() == 4);
     843           1 :             CATCH_REQUIRE(b.columns() == 4);
     844             : 
     845           1 :             double const b00(frand());
     846           1 :             double const b01(frand());
     847           1 :             double const b02(frand());
     848           1 :             double const b03(frand());
     849           1 :             double const b10(frand());
     850           1 :             double const b11(frand());
     851           1 :             double const b12(frand());
     852           1 :             double const b13(frand());
     853           1 :             double const b20(frand());
     854           1 :             double const b21(frand());
     855           1 :             double const b22(frand());
     856           1 :             double const b23(frand());
     857           1 :             double const b30(frand());
     858           1 :             double const b31(frand());
     859           1 :             double const b32(frand());
     860           1 :             double const b33(frand());
     861             : 
     862           1 :             b[0][0] = b00;
     863           1 :             b[0][1] = b01;
     864           1 :             b[0][2] = b02;
     865           1 :             b[0][3] = b03;
     866           1 :             b[1][0] = b10;
     867           1 :             b[1][1] = b11;
     868           1 :             b[1][2] = b12;
     869           1 :             b[1][3] = b13;
     870           1 :             b[2][0] = b20;
     871           1 :             b[2][1] = b21;
     872           1 :             b[2][2] = b22;
     873           1 :             b[2][3] = b23;
     874           1 :             b[3][0] = b30;
     875           1 :             b[3][1] = b31;
     876           1 :             b[3][2] = b32;
     877           1 :             b[3][3] = b33;
     878             : 
     879           1 :             CATCH_REQUIRE(b[0][0] == b00);
     880           1 :             CATCH_REQUIRE(b[0][1] == b01);
     881           1 :             CATCH_REQUIRE(b[0][2] == b02);
     882           1 :             CATCH_REQUIRE(b[0][3] == b03);
     883           1 :             CATCH_REQUIRE(b[1][0] == b10);
     884           1 :             CATCH_REQUIRE(b[1][1] == b11);
     885           1 :             CATCH_REQUIRE(b[1][2] == b12);
     886           1 :             CATCH_REQUIRE(b[1][3] == b13);
     887           1 :             CATCH_REQUIRE(b[2][0] == b20);
     888           1 :             CATCH_REQUIRE(b[2][1] == b21);
     889           1 :             CATCH_REQUIRE(b[2][2] == b22);
     890           1 :             CATCH_REQUIRE(b[2][3] == b23);
     891           1 :             CATCH_REQUIRE(b[3][0] == b30);
     892           1 :             CATCH_REQUIRE(b[3][1] == b31);
     893           1 :             CATCH_REQUIRE(b[3][2] == b32);
     894           1 :             CATCH_REQUIRE(b[3][3] == b33);
     895             : 
     896             :             // setup C
     897             :             //
     898           2 :             snapdev::matrix<double> c(4, 4);
     899             : 
     900           1 :             CATCH_REQUIRE(c.rows() == 4);
     901           1 :             CATCH_REQUIRE(c.columns() == 4);
     902             : 
     903           1 :             double const c00(frand());
     904           1 :             double const c01(frand());
     905           1 :             double const c02(frand());
     906           1 :             double const c03(frand());
     907           1 :             double const c10(frand());
     908           1 :             double const c11(frand());
     909           1 :             double const c12(frand());
     910           1 :             double const c13(frand());
     911           1 :             double const c20(frand());
     912           1 :             double const c21(frand());
     913           1 :             double const c22(frand());
     914           1 :             double const c23(frand());
     915           1 :             double const c30(frand());
     916           1 :             double const c31(frand());
     917           1 :             double const c32(frand());
     918           1 :             double const c33(frand());
     919             : 
     920           1 :             c[0][0] = c00;
     921           1 :             c[0][1] = c01;
     922           1 :             c[0][2] = c02;
     923           1 :             c[0][3] = c03;
     924           1 :             c[1][0] = c10;
     925           1 :             c[1][1] = c11;
     926           1 :             c[1][2] = c12;
     927           1 :             c[1][3] = c13;
     928           1 :             c[2][0] = c20;
     929           1 :             c[2][1] = c21;
     930           1 :             c[2][2] = c22;
     931           1 :             c[2][3] = c23;
     932           1 :             c[3][0] = c30;
     933           1 :             c[3][1] = c31;
     934           1 :             c[3][2] = c32;
     935           1 :             c[3][3] = c33;
     936             : 
     937           1 :             CATCH_REQUIRE(c[0][0] == c00);
     938           1 :             CATCH_REQUIRE(c[0][1] == c01);
     939           1 :             CATCH_REQUIRE(c[0][2] == c02);
     940           1 :             CATCH_REQUIRE(c[0][3] == c03);
     941           1 :             CATCH_REQUIRE(c[1][0] == c10);
     942           1 :             CATCH_REQUIRE(c[1][1] == c11);
     943           1 :             CATCH_REQUIRE(c[1][2] == c12);
     944           1 :             CATCH_REQUIRE(c[1][3] == c13);
     945           1 :             CATCH_REQUIRE(c[2][0] == c20);
     946           1 :             CATCH_REQUIRE(c[2][1] == c21);
     947           1 :             CATCH_REQUIRE(c[2][2] == c22);
     948           1 :             CATCH_REQUIRE(c[2][3] == c23);
     949           1 :             CATCH_REQUIRE(c[3][0] == c30);
     950           1 :             CATCH_REQUIRE(c[3][1] == c31);
     951           1 :             CATCH_REQUIRE(c[3][2] == c32);
     952           1 :             CATCH_REQUIRE(c[3][3] == c33);
     953             : 
     954             :             // run operation C = A + B
     955             :             //
     956           1 :             c = a + b;
     957             : 
     958           1 :             CATCH_REQUIRE(a[0][0] == a00);
     959           1 :             CATCH_REQUIRE(a[0][1] == a01);
     960           1 :             CATCH_REQUIRE(a[0][2] == a02);
     961           1 :             CATCH_REQUIRE(a[0][3] == a03);
     962           1 :             CATCH_REQUIRE(a[1][0] == a10);
     963           1 :             CATCH_REQUIRE(a[1][1] == a11);
     964           1 :             CATCH_REQUIRE(a[1][2] == a12);
     965           1 :             CATCH_REQUIRE(a[1][3] == a13);
     966           1 :             CATCH_REQUIRE(a[2][0] == a20);
     967           1 :             CATCH_REQUIRE(a[2][1] == a21);
     968           1 :             CATCH_REQUIRE(a[2][2] == a22);
     969           1 :             CATCH_REQUIRE(a[2][3] == a23);
     970           1 :             CATCH_REQUIRE(a[3][0] == a30);
     971           1 :             CATCH_REQUIRE(a[3][1] == a31);
     972           1 :             CATCH_REQUIRE(a[3][2] == a32);
     973           1 :             CATCH_REQUIRE(a[3][3] == a33);
     974             : 
     975           1 :             CATCH_REQUIRE(b[0][0] == b00);
     976           1 :             CATCH_REQUIRE(b[0][1] == b01);
     977           1 :             CATCH_REQUIRE(b[0][2] == b02);
     978           1 :             CATCH_REQUIRE(b[0][3] == b03);
     979           1 :             CATCH_REQUIRE(b[1][0] == b10);
     980           1 :             CATCH_REQUIRE(b[1][1] == b11);
     981           1 :             CATCH_REQUIRE(b[1][2] == b12);
     982           1 :             CATCH_REQUIRE(b[1][3] == b13);
     983           1 :             CATCH_REQUIRE(b[2][0] == b20);
     984           1 :             CATCH_REQUIRE(b[2][1] == b21);
     985           1 :             CATCH_REQUIRE(b[2][2] == b22);
     986           1 :             CATCH_REQUIRE(b[2][3] == b23);
     987           1 :             CATCH_REQUIRE(b[3][0] == b30);
     988           1 :             CATCH_REQUIRE(b[3][1] == b31);
     989           1 :             CATCH_REQUIRE(b[3][2] == b32);
     990           1 :             CATCH_REQUIRE(b[3][3] == b33);
     991             : 
     992             :             // this could fail
     993             :             //CATCH_REQUIRE_FALSE(c[0][0] == c00);
     994             :             //CATCH_REQUIRE_FALSE(c[0][1] == c01);
     995             :             //CATCH_REQUIRE_FALSE(c[0][2] == c02);
     996             :             //CATCH_REQUIRE_FALSE(c[0][3] == c03);
     997             :             //CATCH_REQUIRE_FALSE(c[1][0] == c10);
     998             :             //CATCH_REQUIRE_FALSE(c[1][1] == c11);
     999             :             //CATCH_REQUIRE_FALSE(c[1][2] == c12);
    1000             :             //CATCH_REQUIRE_FALSE(c[1][3] == c13);
    1001             :             //CATCH_REQUIRE_FALSE(c[2][0] == c20);
    1002             :             //CATCH_REQUIRE_FALSE(c[2][1] == c21);
    1003             :             //CATCH_REQUIRE_FALSE(c[2][2] == c22);
    1004             :             //CATCH_REQUIRE_FALSE(c[2][3] == c23);
    1005             :             //CATCH_REQUIRE_FALSE(c[3][0] == c30);
    1006             :             //CATCH_REQUIRE_FALSE(c[3][1] == c31);
    1007             :             //CATCH_REQUIRE_FALSE(c[3][2] == c32);
    1008             :             //CATCH_REQUIRE_FALSE(c[3][3] == c33);
    1009             : 
    1010           1 :             CATCH_REQUIRE(c[0][0] == a00 + b00);
    1011           1 :             CATCH_REQUIRE(c[0][1] == a01 + b01);
    1012           1 :             CATCH_REQUIRE(c[0][2] == a02 + b02);
    1013           1 :             CATCH_REQUIRE(c[0][3] == a03 + b03);
    1014           1 :             CATCH_REQUIRE(c[1][0] == a10 + b10);
    1015           1 :             CATCH_REQUIRE(c[1][1] == a11 + b11);
    1016           1 :             CATCH_REQUIRE(c[1][2] == a12 + b12);
    1017           1 :             CATCH_REQUIRE(c[1][3] == a13 + b13);
    1018           1 :             CATCH_REQUIRE(c[2][0] == a20 + b20);
    1019           1 :             CATCH_REQUIRE(c[2][1] == a21 + b21);
    1020           1 :             CATCH_REQUIRE(c[2][2] == a22 + b22);
    1021           1 :             CATCH_REQUIRE(c[2][3] == a23 + b23);
    1022           1 :             CATCH_REQUIRE(c[3][0] == a30 + b30);
    1023           1 :             CATCH_REQUIRE(c[3][1] == a31 + b31);
    1024           1 :             CATCH_REQUIRE(c[3][2] == a32 + b32);
    1025           1 :             CATCH_REQUIRE(c[3][3] == a33 + b33);
    1026             :         }
    1027             :         CATCH_END_SECTION()
    1028             : 
    1029           8 :         CATCH_START_SECTION("matrix: a+=b")
    1030             :         {
    1031             :             // setup A
    1032             :             //
    1033           2 :             snapdev::matrix<double> a(4, 4);
    1034             : 
    1035           1 :             CATCH_REQUIRE(a.rows() == 4);
    1036           1 :             CATCH_REQUIRE(a.columns() == 4);
    1037             : 
    1038           1 :             double const a00(frand());
    1039           1 :             double const a01(frand());
    1040           1 :             double const a02(frand());
    1041           1 :             double const a03(frand());
    1042           1 :             double const a10(frand());
    1043           1 :             double const a11(frand());
    1044           1 :             double const a12(frand());
    1045           1 :             double const a13(frand());
    1046           1 :             double const a20(frand());
    1047           1 :             double const a21(frand());
    1048           1 :             double const a22(frand());
    1049           1 :             double const a23(frand());
    1050           1 :             double const a30(frand());
    1051           1 :             double const a31(frand());
    1052           1 :             double const a32(frand());
    1053           1 :             double const a33(frand());
    1054             : 
    1055           1 :             a[0][0] = a00;
    1056           1 :             a[0][1] = a01;
    1057           1 :             a[0][2] = a02;
    1058           1 :             a[0][3] = a03;
    1059           1 :             a[1][0] = a10;
    1060           1 :             a[1][1] = a11;
    1061           1 :             a[1][2] = a12;
    1062           1 :             a[1][3] = a13;
    1063           1 :             a[2][0] = a20;
    1064           1 :             a[2][1] = a21;
    1065           1 :             a[2][2] = a22;
    1066           1 :             a[2][3] = a23;
    1067           1 :             a[3][0] = a30;
    1068           1 :             a[3][1] = a31;
    1069           1 :             a[3][2] = a32;
    1070           1 :             a[3][3] = a33;
    1071             : 
    1072           1 :             CATCH_REQUIRE(a[0][0] == a00);
    1073           1 :             CATCH_REQUIRE(a[0][1] == a01);
    1074           1 :             CATCH_REQUIRE(a[0][2] == a02);
    1075           1 :             CATCH_REQUIRE(a[0][3] == a03);
    1076           1 :             CATCH_REQUIRE(a[1][0] == a10);
    1077           1 :             CATCH_REQUIRE(a[1][1] == a11);
    1078           1 :             CATCH_REQUIRE(a[1][2] == a12);
    1079           1 :             CATCH_REQUIRE(a[1][3] == a13);
    1080           1 :             CATCH_REQUIRE(a[2][0] == a20);
    1081           1 :             CATCH_REQUIRE(a[2][1] == a21);
    1082           1 :             CATCH_REQUIRE(a[2][2] == a22);
    1083           1 :             CATCH_REQUIRE(a[2][3] == a23);
    1084           1 :             CATCH_REQUIRE(a[3][0] == a30);
    1085           1 :             CATCH_REQUIRE(a[3][1] == a31);
    1086           1 :             CATCH_REQUIRE(a[3][2] == a32);
    1087           1 :             CATCH_REQUIRE(a[3][3] == a33);
    1088             : 
    1089             :             // setup B
    1090             :             //
    1091           2 :             snapdev::matrix<double> b(4, 4);
    1092             : 
    1093           1 :             CATCH_REQUIRE(b.rows() == 4);
    1094           1 :             CATCH_REQUIRE(b.columns() == 4);
    1095             : 
    1096           1 :             double const b00(frand());
    1097           1 :             double const b01(frand());
    1098           1 :             double const b02(frand());
    1099           1 :             double const b03(frand());
    1100           1 :             double const b10(frand());
    1101           1 :             double const b11(frand());
    1102           1 :             double const b12(frand());
    1103           1 :             double const b13(frand());
    1104           1 :             double const b20(frand());
    1105           1 :             double const b21(frand());
    1106           1 :             double const b22(frand());
    1107           1 :             double const b23(frand());
    1108           1 :             double const b30(frand());
    1109           1 :             double const b31(frand());
    1110           1 :             double const b32(frand());
    1111           1 :             double const b33(frand());
    1112             : 
    1113           1 :             b[0][0] = b00;
    1114           1 :             b[0][1] = b01;
    1115           1 :             b[0][2] = b02;
    1116           1 :             b[0][3] = b03;
    1117           1 :             b[1][0] = b10;
    1118           1 :             b[1][1] = b11;
    1119           1 :             b[1][2] = b12;
    1120           1 :             b[1][3] = b13;
    1121           1 :             b[2][0] = b20;
    1122           1 :             b[2][1] = b21;
    1123           1 :             b[2][2] = b22;
    1124           1 :             b[2][3] = b23;
    1125           1 :             b[3][0] = b30;
    1126           1 :             b[3][1] = b31;
    1127           1 :             b[3][2] = b32;
    1128           1 :             b[3][3] = b33;
    1129             : 
    1130           1 :             CATCH_REQUIRE(b[0][0] == b00);
    1131           1 :             CATCH_REQUIRE(b[0][1] == b01);
    1132           1 :             CATCH_REQUIRE(b[0][2] == b02);
    1133           1 :             CATCH_REQUIRE(b[0][3] == b03);
    1134           1 :             CATCH_REQUIRE(b[1][0] == b10);
    1135           1 :             CATCH_REQUIRE(b[1][1] == b11);
    1136           1 :             CATCH_REQUIRE(b[1][2] == b12);
    1137           1 :             CATCH_REQUIRE(b[1][3] == b13);
    1138           1 :             CATCH_REQUIRE(b[2][0] == b20);
    1139           1 :             CATCH_REQUIRE(b[2][1] == b21);
    1140           1 :             CATCH_REQUIRE(b[2][2] == b22);
    1141           1 :             CATCH_REQUIRE(b[2][3] == b23);
    1142           1 :             CATCH_REQUIRE(b[3][0] == b30);
    1143           1 :             CATCH_REQUIRE(b[3][1] == b31);
    1144           1 :             CATCH_REQUIRE(b[3][2] == b32);
    1145           1 :             CATCH_REQUIRE(b[3][3] == b33);
    1146             : 
    1147             :             // run operation A += B
    1148             :             //
    1149           1 :             a += b;
    1150             : 
    1151             :             // this could fail if any bXX is 0.0
    1152             :             //
    1153             :             //CATCH_REQUIRE_FALSE(a[0][0] == a00);
    1154             :             //CATCH_REQUIRE_FALSE(a[0][1] == a01);
    1155             :             //CATCH_REQUIRE_FALSE(a[0][2] == a02);
    1156             :             //CATCH_REQUIRE_FALSE(a[0][3] == a03);
    1157             :             //CATCH_REQUIRE_FALSE(a[1][0] == a10);
    1158             :             //CATCH_REQUIRE_FALSE(a[1][1] == a11);
    1159             :             //CATCH_REQUIRE_FALSE(a[1][2] == a12);
    1160             :             //CATCH_REQUIRE_FALSE(a[1][3] == a13);
    1161             :             //CATCH_REQUIRE_FALSE(a[2][0] == a20);
    1162             :             //CATCH_REQUIRE_FALSE(a[2][1] == a21);
    1163             :             //CATCH_REQUIRE_FALSE(a[2][2] == a22);
    1164             :             //CATCH_REQUIRE_FALSE(a[2][3] == a23);
    1165             :             //CATCH_REQUIRE_FALSE(a[3][0] == a30);
    1166             :             //CATCH_REQUIRE_FALSE(a[3][1] == a31);
    1167             :             //CATCH_REQUIRE_FALSE(a[3][2] == a32);
    1168             :             //CATCH_REQUIRE_FALSE(a[3][3] == a33);
    1169             : 
    1170           1 :             CATCH_REQUIRE(b[0][0] == b00);
    1171           1 :             CATCH_REQUIRE(b[0][1] == b01);
    1172           1 :             CATCH_REQUIRE(b[0][2] == b02);
    1173           1 :             CATCH_REQUIRE(b[0][3] == b03);
    1174           1 :             CATCH_REQUIRE(b[1][0] == b10);
    1175           1 :             CATCH_REQUIRE(b[1][1] == b11);
    1176           1 :             CATCH_REQUIRE(b[1][2] == b12);
    1177           1 :             CATCH_REQUIRE(b[1][3] == b13);
    1178           1 :             CATCH_REQUIRE(b[2][0] == b20);
    1179           1 :             CATCH_REQUIRE(b[2][1] == b21);
    1180           1 :             CATCH_REQUIRE(b[2][2] == b22);
    1181           1 :             CATCH_REQUIRE(b[2][3] == b23);
    1182           1 :             CATCH_REQUIRE(b[3][0] == b30);
    1183           1 :             CATCH_REQUIRE(b[3][1] == b31);
    1184           1 :             CATCH_REQUIRE(b[3][2] == b32);
    1185           1 :             CATCH_REQUIRE(b[3][3] == b33);
    1186             : 
    1187           1 :             CATCH_REQUIRE(a[0][0] == a00 + b00);
    1188           1 :             CATCH_REQUIRE(a[0][1] == a01 + b01);
    1189           1 :             CATCH_REQUIRE(a[0][2] == a02 + b02);
    1190           1 :             CATCH_REQUIRE(a[0][3] == a03 + b03);
    1191           1 :             CATCH_REQUIRE(a[1][0] == a10 + b10);
    1192           1 :             CATCH_REQUIRE(a[1][1] == a11 + b11);
    1193           1 :             CATCH_REQUIRE(a[1][2] == a12 + b12);
    1194           1 :             CATCH_REQUIRE(a[1][3] == a13 + b13);
    1195           1 :             CATCH_REQUIRE(a[2][0] == a20 + b20);
    1196           1 :             CATCH_REQUIRE(a[2][1] == a21 + b21);
    1197           1 :             CATCH_REQUIRE(a[2][2] == a22 + b22);
    1198           1 :             CATCH_REQUIRE(a[2][3] == a23 + b23);
    1199           1 :             CATCH_REQUIRE(a[3][0] == a30 + b30);
    1200           1 :             CATCH_REQUIRE(a[3][1] == a31 + b31);
    1201           1 :             CATCH_REQUIRE(a[3][2] == a32 + b32);
    1202           1 :             CATCH_REQUIRE(a[3][3] == a33 + b33);
    1203             :         }
    1204             :         CATCH_END_SECTION()
    1205             :     }
    1206             : 
    1207             :     // create two random 4x4 matrices and make sure the add works
    1208             :     //
    1209          16 :     CATCH_GIVEN("sub")
    1210             :     {
    1211           8 :         CATCH_START_SECTION("matrix: b=a-<scalar>")
    1212             :         {
    1213             :             // setup A
    1214             :             //
    1215           2 :             snapdev::matrix<double> a(4, 4);
    1216             : 
    1217           1 :             CATCH_REQUIRE(a.rows() == 4);
    1218           1 :             CATCH_REQUIRE(a.columns() == 4);
    1219             : 
    1220           1 :             double const a00(frand());
    1221           1 :             double const a01(frand());
    1222           1 :             double const a02(frand());
    1223           1 :             double const a03(frand());
    1224           1 :             double const a10(frand());
    1225           1 :             double const a11(frand());
    1226           1 :             double const a12(frand());
    1227           1 :             double const a13(frand());
    1228           1 :             double const a20(frand());
    1229           1 :             double const a21(frand());
    1230           1 :             double const a22(frand());
    1231           1 :             double const a23(frand());
    1232           1 :             double const a30(frand());
    1233           1 :             double const a31(frand());
    1234           1 :             double const a32(frand());
    1235           1 :             double const a33(frand());
    1236             : 
    1237           1 :             a[0][0] = a00;
    1238           1 :             a[0][1] = a01;
    1239           1 :             a[0][2] = a02;
    1240           1 :             a[0][3] = a03;
    1241           1 :             a[1][0] = a10;
    1242           1 :             a[1][1] = a11;
    1243           1 :             a[1][2] = a12;
    1244           1 :             a[1][3] = a13;
    1245           1 :             a[2][0] = a20;
    1246           1 :             a[2][1] = a21;
    1247           1 :             a[2][2] = a22;
    1248           1 :             a[2][3] = a23;
    1249           1 :             a[3][0] = a30;
    1250           1 :             a[3][1] = a31;
    1251           1 :             a[3][2] = a32;
    1252           1 :             a[3][3] = a33;
    1253             : 
    1254           1 :             CATCH_REQUIRE(a[0][0] == a00);
    1255           1 :             CATCH_REQUIRE(a[0][1] == a01);
    1256           1 :             CATCH_REQUIRE(a[0][2] == a02);
    1257           1 :             CATCH_REQUIRE(a[0][3] == a03);
    1258           1 :             CATCH_REQUIRE(a[1][0] == a10);
    1259           1 :             CATCH_REQUIRE(a[1][1] == a11);
    1260           1 :             CATCH_REQUIRE(a[1][2] == a12);
    1261           1 :             CATCH_REQUIRE(a[1][3] == a13);
    1262           1 :             CATCH_REQUIRE(a[2][0] == a20);
    1263           1 :             CATCH_REQUIRE(a[2][1] == a21);
    1264           1 :             CATCH_REQUIRE(a[2][2] == a22);
    1265           1 :             CATCH_REQUIRE(a[2][3] == a23);
    1266           1 :             CATCH_REQUIRE(a[3][0] == a30);
    1267           1 :             CATCH_REQUIRE(a[3][1] == a31);
    1268           1 :             CATCH_REQUIRE(a[3][2] == a32);
    1269           1 :             CATCH_REQUIRE(a[3][3] == a33);
    1270             : 
    1271             :             // setup B
    1272             :             //
    1273           2 :             snapdev::matrix<double> b(4, 4);
    1274             : 
    1275           1 :             CATCH_REQUIRE(b.rows() == 4);
    1276           1 :             CATCH_REQUIRE(b.columns() == 4);
    1277             : 
    1278           1 :             double const b00(frand());
    1279           1 :             double const b01(frand());
    1280           1 :             double const b02(frand());
    1281           1 :             double const b03(frand());
    1282           1 :             double const b10(frand());
    1283           1 :             double const b11(frand());
    1284           1 :             double const b12(frand());
    1285           1 :             double const b13(frand());
    1286           1 :             double const b20(frand());
    1287           1 :             double const b21(frand());
    1288           1 :             double const b22(frand());
    1289           1 :             double const b23(frand());
    1290           1 :             double const b30(frand());
    1291           1 :             double const b31(frand());
    1292           1 :             double const b32(frand());
    1293           1 :             double const b33(frand());
    1294             : 
    1295           1 :             b[0][0] = b00;
    1296           1 :             b[0][1] = b01;
    1297           1 :             b[0][2] = b02;
    1298           1 :             b[0][3] = b03;
    1299           1 :             b[1][0] = b10;
    1300           1 :             b[1][1] = b11;
    1301           1 :             b[1][2] = b12;
    1302           1 :             b[1][3] = b13;
    1303           1 :             b[2][0] = b20;
    1304           1 :             b[2][1] = b21;
    1305           1 :             b[2][2] = b22;
    1306           1 :             b[2][3] = b23;
    1307           1 :             b[3][0] = b30;
    1308           1 :             b[3][1] = b31;
    1309           1 :             b[3][2] = b32;
    1310           1 :             b[3][3] = b33;
    1311             : 
    1312           1 :             CATCH_REQUIRE(b[0][0] == b00);
    1313           1 :             CATCH_REQUIRE(b[0][1] == b01);
    1314           1 :             CATCH_REQUIRE(b[0][2] == b02);
    1315           1 :             CATCH_REQUIRE(b[0][3] == b03);
    1316           1 :             CATCH_REQUIRE(b[1][0] == b10);
    1317           1 :             CATCH_REQUIRE(b[1][1] == b11);
    1318           1 :             CATCH_REQUIRE(b[1][2] == b12);
    1319           1 :             CATCH_REQUIRE(b[1][3] == b13);
    1320           1 :             CATCH_REQUIRE(b[2][0] == b20);
    1321           1 :             CATCH_REQUIRE(b[2][1] == b21);
    1322           1 :             CATCH_REQUIRE(b[2][2] == b22);
    1323           1 :             CATCH_REQUIRE(b[2][3] == b23);
    1324           1 :             CATCH_REQUIRE(b[3][0] == b30);
    1325           1 :             CATCH_REQUIRE(b[3][1] == b31);
    1326           1 :             CATCH_REQUIRE(b[3][2] == b32);
    1327           1 :             CATCH_REQUIRE(b[3][3] == b33);
    1328             : 
    1329             :             // create a scalar for our test
    1330             :             //
    1331           1 :             double scalar(frand());
    1332           1 :             while(scalar == 0.0)
    1333             :             {
    1334           0 :                 scalar = frand();
    1335             :             }
    1336             : 
    1337             :             // run operation B = A - <scalar>
    1338             :             //
    1339           1 :             b = a - scalar;
    1340             : 
    1341           1 :             CATCH_REQUIRE(a[0][0] == a00);
    1342           1 :             CATCH_REQUIRE(a[0][1] == a01);
    1343           1 :             CATCH_REQUIRE(a[0][2] == a02);
    1344           1 :             CATCH_REQUIRE(a[0][3] == a03);
    1345           1 :             CATCH_REQUIRE(a[1][0] == a10);
    1346           1 :             CATCH_REQUIRE(a[1][1] == a11);
    1347           1 :             CATCH_REQUIRE(a[1][2] == a12);
    1348           1 :             CATCH_REQUIRE(a[1][3] == a13);
    1349           1 :             CATCH_REQUIRE(a[2][0] == a20);
    1350           1 :             CATCH_REQUIRE(a[2][1] == a21);
    1351           1 :             CATCH_REQUIRE(a[2][2] == a22);
    1352           1 :             CATCH_REQUIRE(a[2][3] == a23);
    1353           1 :             CATCH_REQUIRE(a[3][0] == a30);
    1354           1 :             CATCH_REQUIRE(a[3][1] == a31);
    1355           1 :             CATCH_REQUIRE(a[3][2] == a32);
    1356           1 :             CATCH_REQUIRE(a[3][3] == a33);
    1357             : 
    1358             :             // this can't fail because we ensure scalar != 0
    1359             :             //
    1360           1 :             CATCH_REQUIRE_FALSE(b[0][0] == b00);
    1361           1 :             CATCH_REQUIRE_FALSE(b[0][1] == b01);
    1362           1 :             CATCH_REQUIRE_FALSE(b[0][2] == b02);
    1363           1 :             CATCH_REQUIRE_FALSE(b[0][3] == b03);
    1364           1 :             CATCH_REQUIRE_FALSE(b[1][0] == b10);
    1365           1 :             CATCH_REQUIRE_FALSE(b[1][1] == b11);
    1366           1 :             CATCH_REQUIRE_FALSE(b[1][2] == b12);
    1367           1 :             CATCH_REQUIRE_FALSE(b[1][3] == b13);
    1368           1 :             CATCH_REQUIRE_FALSE(b[2][0] == b20);
    1369           1 :             CATCH_REQUIRE_FALSE(b[2][1] == b21);
    1370           1 :             CATCH_REQUIRE_FALSE(b[2][2] == b22);
    1371           1 :             CATCH_REQUIRE_FALSE(b[2][3] == b23);
    1372           1 :             CATCH_REQUIRE_FALSE(b[3][0] == b30);
    1373           1 :             CATCH_REQUIRE_FALSE(b[3][1] == b31);
    1374           1 :             CATCH_REQUIRE_FALSE(b[3][2] == b32);
    1375           1 :             CATCH_REQUIRE_FALSE(b[3][3] == b33);
    1376             : 
    1377           1 :             CATCH_REQUIRE(b[0][0] == a00 - scalar);
    1378           1 :             CATCH_REQUIRE(b[0][1] == a01 - scalar);
    1379           1 :             CATCH_REQUIRE(b[0][2] == a02 - scalar);
    1380           1 :             CATCH_REQUIRE(b[0][3] == a03 - scalar);
    1381           1 :             CATCH_REQUIRE(b[1][0] == a10 - scalar);
    1382           1 :             CATCH_REQUIRE(b[1][1] == a11 - scalar);
    1383           1 :             CATCH_REQUIRE(b[1][2] == a12 - scalar);
    1384           1 :             CATCH_REQUIRE(b[1][3] == a13 - scalar);
    1385           1 :             CATCH_REQUIRE(b[2][0] == a20 - scalar);
    1386           1 :             CATCH_REQUIRE(b[2][1] == a21 - scalar);
    1387           1 :             CATCH_REQUIRE(b[2][2] == a22 - scalar);
    1388           1 :             CATCH_REQUIRE(b[2][3] == a23 - scalar);
    1389           1 :             CATCH_REQUIRE(b[3][0] == a30 - scalar);
    1390           1 :             CATCH_REQUIRE(b[3][1] == a31 - scalar);
    1391           1 :             CATCH_REQUIRE(b[3][2] == a32 - scalar);
    1392           1 :             CATCH_REQUIRE(b[3][3] == a33 - scalar);
    1393             :         }
    1394             :         CATCH_END_SECTION()
    1395             : 
    1396           8 :         CATCH_START_SECTION("matrix: a-=<scalar>")
    1397             :         {
    1398             :             // setup A
    1399             :             //
    1400           2 :             snapdev::matrix<double> a(4, 4);
    1401             : 
    1402           1 :             CATCH_REQUIRE(a.rows() == 4);
    1403           1 :             CATCH_REQUIRE(a.columns() == 4);
    1404             : 
    1405           1 :             double const a00(frand());
    1406           1 :             double const a01(frand());
    1407           1 :             double const a02(frand());
    1408           1 :             double const a03(frand());
    1409           1 :             double const a10(frand());
    1410           1 :             double const a11(frand());
    1411           1 :             double const a12(frand());
    1412           1 :             double const a13(frand());
    1413           1 :             double const a20(frand());
    1414           1 :             double const a21(frand());
    1415           1 :             double const a22(frand());
    1416           1 :             double const a23(frand());
    1417           1 :             double const a30(frand());
    1418           1 :             double const a31(frand());
    1419           1 :             double const a32(frand());
    1420           1 :             double const a33(frand());
    1421             : 
    1422           1 :             a[0][0] = a00;
    1423           1 :             a[0][1] = a01;
    1424           1 :             a[0][2] = a02;
    1425           1 :             a[0][3] = a03;
    1426           1 :             a[1][0] = a10;
    1427           1 :             a[1][1] = a11;
    1428           1 :             a[1][2] = a12;
    1429           1 :             a[1][3] = a13;
    1430           1 :             a[2][0] = a20;
    1431           1 :             a[2][1] = a21;
    1432           1 :             a[2][2] = a22;
    1433           1 :             a[2][3] = a23;
    1434           1 :             a[3][0] = a30;
    1435           1 :             a[3][1] = a31;
    1436           1 :             a[3][2] = a32;
    1437           1 :             a[3][3] = a33;
    1438             : 
    1439           1 :             CATCH_REQUIRE(a[0][0] == a00);
    1440           1 :             CATCH_REQUIRE(a[0][1] == a01);
    1441           1 :             CATCH_REQUIRE(a[0][2] == a02);
    1442           1 :             CATCH_REQUIRE(a[0][3] == a03);
    1443           1 :             CATCH_REQUIRE(a[1][0] == a10);
    1444           1 :             CATCH_REQUIRE(a[1][1] == a11);
    1445           1 :             CATCH_REQUIRE(a[1][2] == a12);
    1446           1 :             CATCH_REQUIRE(a[1][3] == a13);
    1447           1 :             CATCH_REQUIRE(a[2][0] == a20);
    1448           1 :             CATCH_REQUIRE(a[2][1] == a21);
    1449           1 :             CATCH_REQUIRE(a[2][2] == a22);
    1450           1 :             CATCH_REQUIRE(a[2][3] == a23);
    1451           1 :             CATCH_REQUIRE(a[3][0] == a30);
    1452           1 :             CATCH_REQUIRE(a[3][1] == a31);
    1453           1 :             CATCH_REQUIRE(a[3][2] == a32);
    1454           1 :             CATCH_REQUIRE(a[3][3] == a33);
    1455             : 
    1456             :             // create a scalar for our test
    1457             :             //
    1458           1 :             double scalar(frand());
    1459           1 :             while(scalar == 0.0)
    1460             :             {
    1461           0 :                 scalar = frand();
    1462             :             }
    1463             : 
    1464             :             // run operation A -= <scalar>
    1465             :             //
    1466           1 :             a -= scalar;
    1467             : 
    1468             :             // this can't fail because we ensure scalar != 0
    1469             :             //
    1470           1 :             CATCH_REQUIRE_FALSE(a[0][0] == a00);
    1471           1 :             CATCH_REQUIRE_FALSE(a[0][1] == a01);
    1472           1 :             CATCH_REQUIRE_FALSE(a[0][2] == a02);
    1473           1 :             CATCH_REQUIRE_FALSE(a[0][3] == a03);
    1474           1 :             CATCH_REQUIRE_FALSE(a[1][0] == a10);
    1475           1 :             CATCH_REQUIRE_FALSE(a[1][1] == a11);
    1476           1 :             CATCH_REQUIRE_FALSE(a[1][2] == a12);
    1477           1 :             CATCH_REQUIRE_FALSE(a[1][3] == a13);
    1478           1 :             CATCH_REQUIRE_FALSE(a[2][0] == a20);
    1479           1 :             CATCH_REQUIRE_FALSE(a[2][1] == a21);
    1480           1 :             CATCH_REQUIRE_FALSE(a[2][2] == a22);
    1481           1 :             CATCH_REQUIRE_FALSE(a[2][3] == a23);
    1482           1 :             CATCH_REQUIRE_FALSE(a[3][0] == a30);
    1483           1 :             CATCH_REQUIRE_FALSE(a[3][1] == a31);
    1484           1 :             CATCH_REQUIRE_FALSE(a[3][2] == a32);
    1485           1 :             CATCH_REQUIRE_FALSE(a[3][3] == a33);
    1486             : 
    1487           1 :             CATCH_REQUIRE(a[0][0] == a00 - scalar);
    1488           1 :             CATCH_REQUIRE(a[0][1] == a01 - scalar);
    1489           1 :             CATCH_REQUIRE(a[0][2] == a02 - scalar);
    1490           1 :             CATCH_REQUIRE(a[0][3] == a03 - scalar);
    1491           1 :             CATCH_REQUIRE(a[1][0] == a10 - scalar);
    1492           1 :             CATCH_REQUIRE(a[1][1] == a11 - scalar);
    1493           1 :             CATCH_REQUIRE(a[1][2] == a12 - scalar);
    1494           1 :             CATCH_REQUIRE(a[1][3] == a13 - scalar);
    1495           1 :             CATCH_REQUIRE(a[2][0] == a20 - scalar);
    1496           1 :             CATCH_REQUIRE(a[2][1] == a21 - scalar);
    1497           1 :             CATCH_REQUIRE(a[2][2] == a22 - scalar);
    1498           1 :             CATCH_REQUIRE(a[2][3] == a23 - scalar);
    1499           1 :             CATCH_REQUIRE(a[3][0] == a30 - scalar);
    1500           1 :             CATCH_REQUIRE(a[3][1] == a31 - scalar);
    1501           1 :             CATCH_REQUIRE(a[3][2] == a32 - scalar);
    1502           1 :             CATCH_REQUIRE(a[3][3] == a33 - scalar);
    1503             :         }
    1504             :         CATCH_END_SECTION()
    1505             : 
    1506           8 :         CATCH_START_SECTION("matrix: c=a-b")
    1507             :         {
    1508             :             // setup A
    1509             :             //
    1510           2 :             snapdev::matrix<double> a(4, 4);
    1511             : 
    1512           1 :             CATCH_REQUIRE(a.rows() == 4);
    1513           1 :             CATCH_REQUIRE(a.columns() == 4);
    1514             : 
    1515           1 :             double const a00(frand());
    1516           1 :             double const a01(frand());
    1517           1 :             double const a02(frand());
    1518           1 :             double const a03(frand());
    1519           1 :             double const a10(frand());
    1520           1 :             double const a11(frand());
    1521           1 :             double const a12(frand());
    1522           1 :             double const a13(frand());
    1523           1 :             double const a20(frand());
    1524           1 :             double const a21(frand());
    1525           1 :             double const a22(frand());
    1526           1 :             double const a23(frand());
    1527           1 :             double const a30(frand());
    1528           1 :             double const a31(frand());
    1529           1 :             double const a32(frand());
    1530           1 :             double const a33(frand());
    1531             : 
    1532           1 :             a[0][0] = a00;
    1533           1 :             a[0][1] = a01;
    1534           1 :             a[0][2] = a02;
    1535           1 :             a[0][3] = a03;
    1536           1 :             a[1][0] = a10;
    1537           1 :             a[1][1] = a11;
    1538           1 :             a[1][2] = a12;
    1539           1 :             a[1][3] = a13;
    1540           1 :             a[2][0] = a20;
    1541           1 :             a[2][1] = a21;
    1542           1 :             a[2][2] = a22;
    1543           1 :             a[2][3] = a23;
    1544           1 :             a[3][0] = a30;
    1545           1 :             a[3][1] = a31;
    1546           1 :             a[3][2] = a32;
    1547           1 :             a[3][3] = a33;
    1548             : 
    1549           1 :             CATCH_REQUIRE(a[0][0] == a00);
    1550           1 :             CATCH_REQUIRE(a[0][1] == a01);
    1551           1 :             CATCH_REQUIRE(a[0][2] == a02);
    1552           1 :             CATCH_REQUIRE(a[0][3] == a03);
    1553           1 :             CATCH_REQUIRE(a[1][0] == a10);
    1554           1 :             CATCH_REQUIRE(a[1][1] == a11);
    1555           1 :             CATCH_REQUIRE(a[1][2] == a12);
    1556           1 :             CATCH_REQUIRE(a[1][3] == a13);
    1557           1 :             CATCH_REQUIRE(a[2][0] == a20);
    1558           1 :             CATCH_REQUIRE(a[2][1] == a21);
    1559           1 :             CATCH_REQUIRE(a[2][2] == a22);
    1560           1 :             CATCH_REQUIRE(a[2][3] == a23);
    1561           1 :             CATCH_REQUIRE(a[3][0] == a30);
    1562           1 :             CATCH_REQUIRE(a[3][1] == a31);
    1563           1 :             CATCH_REQUIRE(a[3][2] == a32);
    1564           1 :             CATCH_REQUIRE(a[3][3] == a33);
    1565             : 
    1566             :             // setup B
    1567             :             //
    1568           2 :             snapdev::matrix<double> b(4, 4);
    1569             : 
    1570           1 :             CATCH_REQUIRE(b.rows() == 4);
    1571           1 :             CATCH_REQUIRE(b.columns() == 4);
    1572             : 
    1573           1 :             double const b00(frand());
    1574           1 :             double const b01(frand());
    1575           1 :             double const b02(frand());
    1576           1 :             double const b03(frand());
    1577           1 :             double const b10(frand());
    1578           1 :             double const b11(frand());
    1579           1 :             double const b12(frand());
    1580           1 :             double const b13(frand());
    1581           1 :             double const b20(frand());
    1582           1 :             double const b21(frand());
    1583           1 :             double const b22(frand());
    1584           1 :             double const b23(frand());
    1585           1 :             double const b30(frand());
    1586           1 :             double const b31(frand());
    1587           1 :             double const b32(frand());
    1588           1 :             double const b33(frand());
    1589             : 
    1590           1 :             b[0][0] = b00;
    1591           1 :             b[0][1] = b01;
    1592           1 :             b[0][2] = b02;
    1593           1 :             b[0][3] = b03;
    1594           1 :             b[1][0] = b10;
    1595           1 :             b[1][1] = b11;
    1596           1 :             b[1][2] = b12;
    1597           1 :             b[1][3] = b13;
    1598           1 :             b[2][0] = b20;
    1599           1 :             b[2][1] = b21;
    1600           1 :             b[2][2] = b22;
    1601           1 :             b[2][3] = b23;
    1602           1 :             b[3][0] = b30;
    1603           1 :             b[3][1] = b31;
    1604           1 :             b[3][2] = b32;
    1605           1 :             b[3][3] = b33;
    1606             : 
    1607           1 :             CATCH_REQUIRE(b[0][0] == b00);
    1608           1 :             CATCH_REQUIRE(b[0][1] == b01);
    1609           1 :             CATCH_REQUIRE(b[0][2] == b02);
    1610           1 :             CATCH_REQUIRE(b[0][3] == b03);
    1611           1 :             CATCH_REQUIRE(b[1][0] == b10);
    1612           1 :             CATCH_REQUIRE(b[1][1] == b11);
    1613           1 :             CATCH_REQUIRE(b[1][2] == b12);
    1614           1 :             CATCH_REQUIRE(b[1][3] == b13);
    1615           1 :             CATCH_REQUIRE(b[2][0] == b20);
    1616           1 :             CATCH_REQUIRE(b[2][1] == b21);
    1617           1 :             CATCH_REQUIRE(b[2][2] == b22);
    1618           1 :             CATCH_REQUIRE(b[2][3] == b23);
    1619           1 :             CATCH_REQUIRE(b[3][0] == b30);
    1620           1 :             CATCH_REQUIRE(b[3][1] == b31);
    1621           1 :             CATCH_REQUIRE(b[3][2] == b32);
    1622           1 :             CATCH_REQUIRE(b[3][3] == b33);
    1623             : 
    1624             :             // setup C
    1625             :             //
    1626           2 :             snapdev::matrix<double> c(4, 4);
    1627             : 
    1628           1 :             CATCH_REQUIRE(c.rows() == 4);
    1629           1 :             CATCH_REQUIRE(c.columns() == 4);
    1630             : 
    1631           1 :             double const c00(frand());
    1632           1 :             double const c01(frand());
    1633           1 :             double const c02(frand());
    1634           1 :             double const c03(frand());
    1635           1 :             double const c10(frand());
    1636           1 :             double const c11(frand());
    1637           1 :             double const c12(frand());
    1638           1 :             double const c13(frand());
    1639           1 :             double const c20(frand());
    1640           1 :             double const c21(frand());
    1641           1 :             double const c22(frand());
    1642           1 :             double const c23(frand());
    1643           1 :             double const c30(frand());
    1644           1 :             double const c31(frand());
    1645           1 :             double const c32(frand());
    1646           1 :             double const c33(frand());
    1647             : 
    1648           1 :             c[0][0] = c00;
    1649           1 :             c[0][1] = c01;
    1650           1 :             c[0][2] = c02;
    1651           1 :             c[0][3] = c03;
    1652           1 :             c[1][0] = c10;
    1653           1 :             c[1][1] = c11;
    1654           1 :             c[1][2] = c12;
    1655           1 :             c[1][3] = c13;
    1656           1 :             c[2][0] = c20;
    1657           1 :             c[2][1] = c21;
    1658           1 :             c[2][2] = c22;
    1659           1 :             c[2][3] = c23;
    1660           1 :             c[3][0] = c30;
    1661           1 :             c[3][1] = c31;
    1662           1 :             c[3][2] = c32;
    1663           1 :             c[3][3] = c33;
    1664             : 
    1665           1 :             CATCH_REQUIRE(c[0][0] == c00);
    1666           1 :             CATCH_REQUIRE(c[0][1] == c01);
    1667           1 :             CATCH_REQUIRE(c[0][2] == c02);
    1668           1 :             CATCH_REQUIRE(c[0][3] == c03);
    1669           1 :             CATCH_REQUIRE(c[1][0] == c10);
    1670           1 :             CATCH_REQUIRE(c[1][1] == c11);
    1671           1 :             CATCH_REQUIRE(c[1][2] == c12);
    1672           1 :             CATCH_REQUIRE(c[1][3] == c13);
    1673           1 :             CATCH_REQUIRE(c[2][0] == c20);
    1674           1 :             CATCH_REQUIRE(c[2][1] == c21);
    1675           1 :             CATCH_REQUIRE(c[2][2] == c22);
    1676           1 :             CATCH_REQUIRE(c[2][3] == c23);
    1677           1 :             CATCH_REQUIRE(c[3][0] == c30);
    1678           1 :             CATCH_REQUIRE(c[3][1] == c31);
    1679           1 :             CATCH_REQUIRE(c[3][2] == c32);
    1680           1 :             CATCH_REQUIRE(c[3][3] == c33);
    1681             : 
    1682             :             // run operation C = A - B
    1683             :             //
    1684           1 :             c = a - b;
    1685             : 
    1686           1 :             CATCH_REQUIRE(a[0][0] == a00);
    1687           1 :             CATCH_REQUIRE(a[0][1] == a01);
    1688           1 :             CATCH_REQUIRE(a[0][2] == a02);
    1689           1 :             CATCH_REQUIRE(a[0][3] == a03);
    1690           1 :             CATCH_REQUIRE(a[1][0] == a10);
    1691           1 :             CATCH_REQUIRE(a[1][1] == a11);
    1692           1 :             CATCH_REQUIRE(a[1][2] == a12);
    1693           1 :             CATCH_REQUIRE(a[1][3] == a13);
    1694           1 :             CATCH_REQUIRE(a[2][0] == a20);
    1695           1 :             CATCH_REQUIRE(a[2][1] == a21);
    1696           1 :             CATCH_REQUIRE(a[2][2] == a22);
    1697           1 :             CATCH_REQUIRE(a[2][3] == a23);
    1698           1 :             CATCH_REQUIRE(a[3][0] == a30);
    1699           1 :             CATCH_REQUIRE(a[3][1] == a31);
    1700           1 :             CATCH_REQUIRE(a[3][2] == a32);
    1701           1 :             CATCH_REQUIRE(a[3][3] == a33);
    1702             : 
    1703           1 :             CATCH_REQUIRE(b[0][0] == b00);
    1704           1 :             CATCH_REQUIRE(b[0][1] == b01);
    1705           1 :             CATCH_REQUIRE(b[0][2] == b02);
    1706           1 :             CATCH_REQUIRE(b[0][3] == b03);
    1707           1 :             CATCH_REQUIRE(b[1][0] == b10);
    1708           1 :             CATCH_REQUIRE(b[1][1] == b11);
    1709           1 :             CATCH_REQUIRE(b[1][2] == b12);
    1710           1 :             CATCH_REQUIRE(b[1][3] == b13);
    1711           1 :             CATCH_REQUIRE(b[2][0] == b20);
    1712           1 :             CATCH_REQUIRE(b[2][1] == b21);
    1713           1 :             CATCH_REQUIRE(b[2][2] == b22);
    1714           1 :             CATCH_REQUIRE(b[2][3] == b23);
    1715           1 :             CATCH_REQUIRE(b[3][0] == b30);
    1716           1 :             CATCH_REQUIRE(b[3][1] == b31);
    1717           1 :             CATCH_REQUIRE(b[3][2] == b32);
    1718           1 :             CATCH_REQUIRE(b[3][3] == b33);
    1719             : 
    1720             :             // this could fail
    1721             :             //CATCH_REQUIRE_FALSE(c[0][0] == c00);
    1722             :             //CATCH_REQUIRE_FALSE(c[0][1] == c01);
    1723             :             //CATCH_REQUIRE_FALSE(c[0][2] == c02);
    1724             :             //CATCH_REQUIRE_FALSE(c[0][3] == c03);
    1725             :             //CATCH_REQUIRE_FALSE(c[1][0] == c10);
    1726             :             //CATCH_REQUIRE_FALSE(c[1][1] == c11);
    1727             :             //CATCH_REQUIRE_FALSE(c[1][2] == c12);
    1728             :             //CATCH_REQUIRE_FALSE(c[1][3] == c13);
    1729             :             //CATCH_REQUIRE_FALSE(c[2][0] == c20);
    1730             :             //CATCH_REQUIRE_FALSE(c[2][1] == c21);
    1731             :             //CATCH_REQUIRE_FALSE(c[2][2] == c22);
    1732             :             //CATCH_REQUIRE_FALSE(c[2][3] == c23);
    1733             :             //CATCH_REQUIRE_FALSE(c[3][0] == c30);
    1734             :             //CATCH_REQUIRE_FALSE(c[3][1] == c31);
    1735             :             //CATCH_REQUIRE_FALSE(c[3][2] == c32);
    1736             :             //CATCH_REQUIRE_FALSE(c[3][3] == c33);
    1737             : 
    1738           1 :             CATCH_REQUIRE(c[0][0] == a00 - b00);
    1739           1 :             CATCH_REQUIRE(c[0][1] == a01 - b01);
    1740           1 :             CATCH_REQUIRE(c[0][2] == a02 - b02);
    1741           1 :             CATCH_REQUIRE(c[0][3] == a03 - b03);
    1742           1 :             CATCH_REQUIRE(c[1][0] == a10 - b10);
    1743           1 :             CATCH_REQUIRE(c[1][1] == a11 - b11);
    1744           1 :             CATCH_REQUIRE(c[1][2] == a12 - b12);
    1745           1 :             CATCH_REQUIRE(c[1][3] == a13 - b13);
    1746           1 :             CATCH_REQUIRE(c[2][0] == a20 - b20);
    1747           1 :             CATCH_REQUIRE(c[2][1] == a21 - b21);
    1748           1 :             CATCH_REQUIRE(c[2][2] == a22 - b22);
    1749           1 :             CATCH_REQUIRE(c[2][3] == a23 - b23);
    1750           1 :             CATCH_REQUIRE(c[3][0] == a30 - b30);
    1751           1 :             CATCH_REQUIRE(c[3][1] == a31 - b31);
    1752           1 :             CATCH_REQUIRE(c[3][2] == a32 - b32);
    1753           1 :             CATCH_REQUIRE(c[3][3] == a33 - b33);
    1754             :         }
    1755             :         CATCH_END_SECTION()
    1756             : 
    1757           8 :         CATCH_START_SECTION("matrix: a-=b")
    1758             :         {
    1759             :             // setup A
    1760             :             //
    1761           2 :             snapdev::matrix<double> a(4, 4);
    1762             : 
    1763           1 :             CATCH_REQUIRE(a.rows() == 4);
    1764           1 :             CATCH_REQUIRE(a.columns() == 4);
    1765             : 
    1766           1 :             double const a00(frand());
    1767           1 :             double const a01(frand());
    1768           1 :             double const a02(frand());
    1769           1 :             double const a03(frand());
    1770           1 :             double const a10(frand());
    1771           1 :             double const a11(frand());
    1772           1 :             double const a12(frand());
    1773           1 :             double const a13(frand());
    1774           1 :             double const a20(frand());
    1775           1 :             double const a21(frand());
    1776           1 :             double const a22(frand());
    1777           1 :             double const a23(frand());
    1778           1 :             double const a30(frand());
    1779           1 :             double const a31(frand());
    1780           1 :             double const a32(frand());
    1781           1 :             double const a33(frand());
    1782             : 
    1783           1 :             a[0][0] = a00;
    1784           1 :             a[0][1] = a01;
    1785           1 :             a[0][2] = a02;
    1786           1 :             a[0][3] = a03;
    1787           1 :             a[1][0] = a10;
    1788           1 :             a[1][1] = a11;
    1789           1 :             a[1][2] = a12;
    1790           1 :             a[1][3] = a13;
    1791           1 :             a[2][0] = a20;
    1792           1 :             a[2][1] = a21;
    1793           1 :             a[2][2] = a22;
    1794           1 :             a[2][3] = a23;
    1795           1 :             a[3][0] = a30;
    1796           1 :             a[3][1] = a31;
    1797           1 :             a[3][2] = a32;
    1798           1 :             a[3][3] = a33;
    1799             : 
    1800           1 :             CATCH_REQUIRE(a[0][0] == a00);
    1801           1 :             CATCH_REQUIRE(a[0][1] == a01);
    1802           1 :             CATCH_REQUIRE(a[0][2] == a02);
    1803           1 :             CATCH_REQUIRE(a[0][3] == a03);
    1804           1 :             CATCH_REQUIRE(a[1][0] == a10);
    1805           1 :             CATCH_REQUIRE(a[1][1] == a11);
    1806           1 :             CATCH_REQUIRE(a[1][2] == a12);
    1807           1 :             CATCH_REQUIRE(a[1][3] == a13);
    1808           1 :             CATCH_REQUIRE(a[2][0] == a20);
    1809           1 :             CATCH_REQUIRE(a[2][1] == a21);
    1810           1 :             CATCH_REQUIRE(a[2][2] == a22);
    1811           1 :             CATCH_REQUIRE(a[2][3] == a23);
    1812           1 :             CATCH_REQUIRE(a[3][0] == a30);
    1813           1 :             CATCH_REQUIRE(a[3][1] == a31);
    1814           1 :             CATCH_REQUIRE(a[3][2] == a32);
    1815           1 :             CATCH_REQUIRE(a[3][3] == a33);
    1816             : 
    1817             :             // setup B
    1818             :             //
    1819           2 :             snapdev::matrix<double> b(4, 4);
    1820             : 
    1821           1 :             CATCH_REQUIRE(b.rows() == 4);
    1822           1 :             CATCH_REQUIRE(b.columns() == 4);
    1823             : 
    1824           1 :             double const b00(frand());
    1825           1 :             double const b01(frand());
    1826           1 :             double const b02(frand());
    1827           1 :             double const b03(frand());
    1828           1 :             double const b10(frand());
    1829           1 :             double const b11(frand());
    1830           1 :             double const b12(frand());
    1831           1 :             double const b13(frand());
    1832           1 :             double const b20(frand());
    1833           1 :             double const b21(frand());
    1834           1 :             double const b22(frand());
    1835           1 :             double const b23(frand());
    1836           1 :             double const b30(frand());
    1837           1 :             double const b31(frand());
    1838           1 :             double const b32(frand());
    1839           1 :             double const b33(frand());
    1840             : 
    1841           1 :             b[0][0] = b00;
    1842           1 :             b[0][1] = b01;
    1843           1 :             b[0][2] = b02;
    1844           1 :             b[0][3] = b03;
    1845           1 :             b[1][0] = b10;
    1846           1 :             b[1][1] = b11;
    1847           1 :             b[1][2] = b12;
    1848           1 :             b[1][3] = b13;
    1849           1 :             b[2][0] = b20;
    1850           1 :             b[2][1] = b21;
    1851           1 :             b[2][2] = b22;
    1852           1 :             b[2][3] = b23;
    1853           1 :             b[3][0] = b30;
    1854           1 :             b[3][1] = b31;
    1855           1 :             b[3][2] = b32;
    1856           1 :             b[3][3] = b33;
    1857             : 
    1858           1 :             CATCH_REQUIRE(b[0][0] == b00);
    1859           1 :             CATCH_REQUIRE(b[0][1] == b01);
    1860           1 :             CATCH_REQUIRE(b[0][2] == b02);
    1861           1 :             CATCH_REQUIRE(b[0][3] == b03);
    1862           1 :             CATCH_REQUIRE(b[1][0] == b10);
    1863           1 :             CATCH_REQUIRE(b[1][1] == b11);
    1864           1 :             CATCH_REQUIRE(b[1][2] == b12);
    1865           1 :             CATCH_REQUIRE(b[1][3] == b13);
    1866           1 :             CATCH_REQUIRE(b[2][0] == b20);
    1867           1 :             CATCH_REQUIRE(b[2][1] == b21);
    1868           1 :             CATCH_REQUIRE(b[2][2] == b22);
    1869           1 :             CATCH_REQUIRE(b[2][3] == b23);
    1870           1 :             CATCH_REQUIRE(b[3][0] == b30);
    1871           1 :             CATCH_REQUIRE(b[3][1] == b31);
    1872           1 :             CATCH_REQUIRE(b[3][2] == b32);
    1873           1 :             CATCH_REQUIRE(b[3][3] == b33);
    1874             : 
    1875             :             // run operation A -= B
    1876             :             //
    1877           1 :             a -= b;
    1878             : 
    1879             :             // this could fail if one of bXX is 0.0
    1880             :             //
    1881             :             //CATCH_REQUIRE(a[0][0] == a00);
    1882             :             //CATCH_REQUIRE(a[0][1] == a01);
    1883             :             //CATCH_REQUIRE(a[0][2] == a02);
    1884             :             //CATCH_REQUIRE(a[0][3] == a03);
    1885             :             //CATCH_REQUIRE(a[1][0] == a10);
    1886             :             //CATCH_REQUIRE(a[1][1] == a11);
    1887             :             //CATCH_REQUIRE(a[1][2] == a12);
    1888             :             //CATCH_REQUIRE(a[1][3] == a13);
    1889             :             //CATCH_REQUIRE(a[2][0] == a20);
    1890             :             //CATCH_REQUIRE(a[2][1] == a21);
    1891             :             //CATCH_REQUIRE(a[2][2] == a22);
    1892             :             //CATCH_REQUIRE(a[2][3] == a23);
    1893             :             //CATCH_REQUIRE(a[3][0] == a30);
    1894             :             //CATCH_REQUIRE(a[3][1] == a31);
    1895             :             //CATCH_REQUIRE(a[3][2] == a32);
    1896             :             //CATCH_REQUIRE(a[3][3] == a33);
    1897             : 
    1898           1 :             CATCH_REQUIRE(b[0][0] == b00);
    1899           1 :             CATCH_REQUIRE(b[0][1] == b01);
    1900           1 :             CATCH_REQUIRE(b[0][2] == b02);
    1901           1 :             CATCH_REQUIRE(b[0][3] == b03);
    1902           1 :             CATCH_REQUIRE(b[1][0] == b10);
    1903           1 :             CATCH_REQUIRE(b[1][1] == b11);
    1904           1 :             CATCH_REQUIRE(b[1][2] == b12);
    1905           1 :             CATCH_REQUIRE(b[1][3] == b13);
    1906           1 :             CATCH_REQUIRE(b[2][0] == b20);
    1907           1 :             CATCH_REQUIRE(b[2][1] == b21);
    1908           1 :             CATCH_REQUIRE(b[2][2] == b22);
    1909           1 :             CATCH_REQUIRE(b[2][3] == b23);
    1910           1 :             CATCH_REQUIRE(b[3][0] == b30);
    1911           1 :             CATCH_REQUIRE(b[3][1] == b31);
    1912           1 :             CATCH_REQUIRE(b[3][2] == b32);
    1913           1 :             CATCH_REQUIRE(b[3][3] == b33);
    1914             : 
    1915           1 :             CATCH_REQUIRE(a[0][0] == a00 - b00);
    1916           1 :             CATCH_REQUIRE(a[0][1] == a01 - b01);
    1917           1 :             CATCH_REQUIRE(a[0][2] == a02 - b02);
    1918           1 :             CATCH_REQUIRE(a[0][3] == a03 - b03);
    1919           1 :             CATCH_REQUIRE(a[1][0] == a10 - b10);
    1920           1 :             CATCH_REQUIRE(a[1][1] == a11 - b11);
    1921           1 :             CATCH_REQUIRE(a[1][2] == a12 - b12);
    1922           1 :             CATCH_REQUIRE(a[1][3] == a13 - b13);
    1923           1 :             CATCH_REQUIRE(a[2][0] == a20 - b20);
    1924           1 :             CATCH_REQUIRE(a[2][1] == a21 - b21);
    1925           1 :             CATCH_REQUIRE(a[2][2] == a22 - b22);
    1926           1 :             CATCH_REQUIRE(a[2][3] == a23 - b23);
    1927           1 :             CATCH_REQUIRE(a[3][0] == a30 - b30);
    1928           1 :             CATCH_REQUIRE(a[3][1] == a31 - b31);
    1929           1 :             CATCH_REQUIRE(a[3][2] == a32 - b32);
    1930           1 :             CATCH_REQUIRE(a[3][3] == a33 - b33);
    1931             :         }
    1932             :         CATCH_END_SECTION()
    1933             :     }
    1934           8 : }
    1935             : 
    1936             : 
    1937          10 : CATCH_TEST_CASE("matrix_util", "[matrix]")
    1938             : {
    1939             :     // various ways to change the data order, minor, enlarge
    1940             :     //
    1941          16 :     CATCH_GIVEN("move data")
    1942             :     {
    1943          10 :         CATCH_START_SECTION("matrix: minor")
    1944             :         {
    1945             :             // setup A
    1946             :             //
    1947           2 :             snapdev::matrix<double> a(4, 4);
    1948             : 
    1949           1 :             CATCH_REQUIRE(a.rows() == 4);
    1950           1 :             CATCH_REQUIRE(a.columns() == 4);
    1951             : 
    1952           1 :             double const a00(frand());
    1953           1 :             double const a01(frand());
    1954           1 :             double const a02(frand());
    1955           1 :             double const a03(frand());
    1956           1 :             double const a10(frand());
    1957           1 :             double const a11(frand());
    1958           1 :             double const a12(frand());
    1959           1 :             double const a13(frand());
    1960           1 :             double const a20(frand());
    1961           1 :             double const a21(frand());
    1962           1 :             double const a22(frand());
    1963           1 :             double const a23(frand());
    1964           1 :             double const a30(frand());
    1965           1 :             double const a31(frand());
    1966           1 :             double const a32(frand());
    1967           1 :             double const a33(frand());
    1968             : 
    1969           1 :             a[0][0] = a00;
    1970           1 :             a[0][1] = a01;
    1971           1 :             a[0][2] = a02;
    1972           1 :             a[0][3] = a03;
    1973           1 :             a[1][0] = a10;
    1974           1 :             a[1][1] = a11;
    1975           1 :             a[1][2] = a12;
    1976           1 :             a[1][3] = a13;
    1977           1 :             a[2][0] = a20;
    1978           1 :             a[2][1] = a21;
    1979           1 :             a[2][2] = a22;
    1980           1 :             a[2][3] = a23;
    1981           1 :             a[3][0] = a30;
    1982           1 :             a[3][1] = a31;
    1983           1 :             a[3][2] = a32;
    1984           1 :             a[3][3] = a33;
    1985             : 
    1986           1 :             CATCH_REQUIRE(a[0][0] == a00);
    1987           1 :             CATCH_REQUIRE(a[0][1] == a01);
    1988           1 :             CATCH_REQUIRE(a[0][2] == a02);
    1989           1 :             CATCH_REQUIRE(a[0][3] == a03);
    1990           1 :             CATCH_REQUIRE(a[1][0] == a10);
    1991           1 :             CATCH_REQUIRE(a[1][1] == a11);
    1992           1 :             CATCH_REQUIRE(a[1][2] == a12);
    1993           1 :             CATCH_REQUIRE(a[1][3] == a13);
    1994           1 :             CATCH_REQUIRE(a[2][0] == a20);
    1995           1 :             CATCH_REQUIRE(a[2][1] == a21);
    1996           1 :             CATCH_REQUIRE(a[2][2] == a22);
    1997           1 :             CATCH_REQUIRE(a[2][3] == a23);
    1998           1 :             CATCH_REQUIRE(a[3][0] == a30);
    1999           1 :             CATCH_REQUIRE(a[3][1] == a31);
    2000           1 :             CATCH_REQUIRE(a[3][2] == a32);
    2001           1 :             CATCH_REQUIRE(a[3][3] == a33);
    2002             : 
    2003             :             {
    2004           2 :                 snapdev::matrix<double> p = a.minor_matrix(0, 0);
    2005             : 
    2006           1 :                 CATCH_REQUIRE(p.rows() == 3);
    2007           1 :                 CATCH_REQUIRE(p.columns() == 3);
    2008             : 
    2009           1 :                 CATCH_REQUIRE(p[0][0] == a11);
    2010           1 :                 CATCH_REQUIRE(p[0][1] == a12);
    2011           1 :                 CATCH_REQUIRE(p[0][2] == a13);
    2012           1 :                 CATCH_REQUIRE(p[1][0] == a21);
    2013           1 :                 CATCH_REQUIRE(p[1][1] == a22);
    2014           1 :                 CATCH_REQUIRE(p[1][2] == a23);
    2015           1 :                 CATCH_REQUIRE(p[2][0] == a31);
    2016           1 :                 CATCH_REQUIRE(p[2][1] == a32);
    2017           1 :                 CATCH_REQUIRE(p[2][2] == a33);
    2018             :             }
    2019             : 
    2020             :             {
    2021           2 :                 snapdev::matrix<double> p = a.minor_matrix(0, 1);
    2022             : 
    2023           1 :                 CATCH_REQUIRE(p.rows() == 3);
    2024           1 :                 CATCH_REQUIRE(p.columns() == 3);
    2025             : 
    2026           1 :                 CATCH_REQUIRE(p[0][0] == a10);
    2027           1 :                 CATCH_REQUIRE(p[0][1] == a12);
    2028           1 :                 CATCH_REQUIRE(p[0][2] == a13);
    2029           1 :                 CATCH_REQUIRE(p[1][0] == a20);
    2030           1 :                 CATCH_REQUIRE(p[1][1] == a22);
    2031           1 :                 CATCH_REQUIRE(p[1][2] == a23);
    2032           1 :                 CATCH_REQUIRE(p[2][0] == a30);
    2033           1 :                 CATCH_REQUIRE(p[2][1] == a32);
    2034           1 :                 CATCH_REQUIRE(p[2][2] == a33);
    2035             :             }
    2036             : 
    2037             :             {
    2038           2 :                 snapdev::matrix<double> p = a.minor_matrix(0, 2);
    2039             : 
    2040           1 :                 CATCH_REQUIRE(p.rows() == 3);
    2041           1 :                 CATCH_REQUIRE(p.columns() == 3);
    2042             : 
    2043           1 :                 CATCH_REQUIRE(p[0][0] == a10);
    2044           1 :                 CATCH_REQUIRE(p[0][1] == a11);
    2045           1 :                 CATCH_REQUIRE(p[0][2] == a13);
    2046           1 :                 CATCH_REQUIRE(p[1][0] == a20);
    2047           1 :                 CATCH_REQUIRE(p[1][1] == a21);
    2048           1 :                 CATCH_REQUIRE(p[1][2] == a23);
    2049           1 :                 CATCH_REQUIRE(p[2][0] == a30);
    2050           1 :                 CATCH_REQUIRE(p[2][1] == a31);
    2051           1 :                 CATCH_REQUIRE(p[2][2] == a33);
    2052             :             }
    2053             : 
    2054             :             {
    2055           2 :                 snapdev::matrix<double> p = a.minor_matrix(0, 3);
    2056             : 
    2057           1 :                 CATCH_REQUIRE(p.rows() == 3);
    2058           1 :                 CATCH_REQUIRE(p.columns() == 3);
    2059             : 
    2060           1 :                 CATCH_REQUIRE(p[0][0] == a10);
    2061           1 :                 CATCH_REQUIRE(p[0][1] == a11);
    2062           1 :                 CATCH_REQUIRE(p[0][2] == a12);
    2063           1 :                 CATCH_REQUIRE(p[1][0] == a20);
    2064           1 :                 CATCH_REQUIRE(p[1][1] == a21);
    2065           1 :                 CATCH_REQUIRE(p[1][2] == a22);
    2066           1 :                 CATCH_REQUIRE(p[2][0] == a30);
    2067           1 :                 CATCH_REQUIRE(p[2][1] == a31);
    2068           1 :                 CATCH_REQUIRE(p[2][2] == a32);
    2069             :             }
    2070             : 
    2071             :             {
    2072           2 :                 snapdev::matrix<double> p = a.minor_matrix(2, 1);
    2073             : 
    2074           1 :                 CATCH_REQUIRE(p.rows() == 3);
    2075           1 :                 CATCH_REQUIRE(p.columns() == 3);
    2076             : 
    2077           1 :                 CATCH_REQUIRE(p[0][0] == a00);
    2078           1 :                 CATCH_REQUIRE(p[0][1] == a02);
    2079           1 :                 CATCH_REQUIRE(p[0][2] == a03);
    2080           1 :                 CATCH_REQUIRE(p[1][0] == a10);
    2081           1 :                 CATCH_REQUIRE(p[1][1] == a12);
    2082           1 :                 CATCH_REQUIRE(p[1][2] == a13);
    2083           1 :                 CATCH_REQUIRE(p[2][0] == a30);
    2084           1 :                 CATCH_REQUIRE(p[2][1] == a32);
    2085           1 :                 CATCH_REQUIRE(p[2][2] == a33);
    2086             :             }
    2087             :         }
    2088             :         CATCH_END_SECTION()
    2089             : 
    2090          10 :         CATCH_START_SECTION("matrix: transpose 4x4")
    2091             :         {
    2092             :             // setup A
    2093             :             //
    2094           2 :             snapdev::matrix<double> a(4, 4);
    2095             : 
    2096           1 :             CATCH_REQUIRE(a.rows() == 4);
    2097           1 :             CATCH_REQUIRE(a.columns() == 4);
    2098             : 
    2099           1 :             double const a00(frand());
    2100           1 :             double const a01(frand());
    2101           1 :             double const a02(frand());
    2102           1 :             double const a03(frand());
    2103           1 :             double const a10(frand());
    2104           1 :             double const a11(frand());
    2105           1 :             double const a12(frand());
    2106           1 :             double const a13(frand());
    2107           1 :             double const a20(frand());
    2108           1 :             double const a21(frand());
    2109           1 :             double const a22(frand());
    2110           1 :             double const a23(frand());
    2111           1 :             double const a30(frand());
    2112           1 :             double const a31(frand());
    2113           1 :             double const a32(frand());
    2114           1 :             double const a33(frand());
    2115             : 
    2116           1 :             a[0][0] = a00;
    2117           1 :             a[0][1] = a01;
    2118           1 :             a[0][2] = a02;
    2119           1 :             a[0][3] = a03;
    2120           1 :             a[1][0] = a10;
    2121           1 :             a[1][1] = a11;
    2122           1 :             a[1][2] = a12;
    2123           1 :             a[1][3] = a13;
    2124           1 :             a[2][0] = a20;
    2125           1 :             a[2][1] = a21;
    2126           1 :             a[2][2] = a22;
    2127           1 :             a[2][3] = a23;
    2128           1 :             a[3][0] = a30;
    2129           1 :             a[3][1] = a31;
    2130           1 :             a[3][2] = a32;
    2131           1 :             a[3][3] = a33;
    2132             : 
    2133           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2134           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2135           1 :             CATCH_REQUIRE(a[0][2] == a02);
    2136           1 :             CATCH_REQUIRE(a[0][3] == a03);
    2137           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2138           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2139           1 :             CATCH_REQUIRE(a[1][2] == a12);
    2140           1 :             CATCH_REQUIRE(a[1][3] == a13);
    2141           1 :             CATCH_REQUIRE(a[2][0] == a20);
    2142           1 :             CATCH_REQUIRE(a[2][1] == a21);
    2143           1 :             CATCH_REQUIRE(a[2][2] == a22);
    2144           1 :             CATCH_REQUIRE(a[2][3] == a23);
    2145           1 :             CATCH_REQUIRE(a[3][0] == a30);
    2146           1 :             CATCH_REQUIRE(a[3][1] == a31);
    2147           1 :             CATCH_REQUIRE(a[3][2] == a32);
    2148           1 :             CATCH_REQUIRE(a[3][3] == a33);
    2149             : 
    2150           2 :             snapdev::matrix<double> t = a.transpose();
    2151             : 
    2152           1 :             CATCH_REQUIRE(t.rows() == 4);
    2153           1 :             CATCH_REQUIRE(t.columns() == 4);
    2154             : 
    2155           1 :             CATCH_REQUIRE(t[0][0] == a00);
    2156           1 :             CATCH_REQUIRE(t[0][1] == a10);
    2157           1 :             CATCH_REQUIRE(t[0][2] == a20);
    2158           1 :             CATCH_REQUIRE(t[0][3] == a30);
    2159           1 :             CATCH_REQUIRE(t[1][0] == a01);
    2160           1 :             CATCH_REQUIRE(t[1][1] == a11);
    2161           1 :             CATCH_REQUIRE(t[1][2] == a21);
    2162           1 :             CATCH_REQUIRE(t[1][3] == a31);
    2163           1 :             CATCH_REQUIRE(t[2][0] == a02);
    2164           1 :             CATCH_REQUIRE(t[2][1] == a12);
    2165           1 :             CATCH_REQUIRE(t[2][2] == a22);
    2166           1 :             CATCH_REQUIRE(t[2][3] == a32);
    2167           1 :             CATCH_REQUIRE(t[3][0] == a03);
    2168           1 :             CATCH_REQUIRE(t[3][1] == a13);
    2169           1 :             CATCH_REQUIRE(t[3][2] == a23);
    2170           1 :             CATCH_REQUIRE(t[3][3] == a33);
    2171             :         }
    2172             :         CATCH_END_SECTION()
    2173             : 
    2174          10 :         CATCH_START_SECTION("matrix: transpose 6x2")
    2175             :         {
    2176             :             // setup A
    2177             :             //
    2178           2 :             snapdev::matrix<double> a(6, 2);
    2179             : 
    2180           1 :             CATCH_REQUIRE(a.rows() == 6);
    2181           1 :             CATCH_REQUIRE(a.columns() == 2);
    2182             : 
    2183           1 :             double const a00(frand());
    2184           1 :             double const a01(frand());
    2185           1 :             double const a10(frand());
    2186           1 :             double const a11(frand());
    2187           1 :             double const a20(frand());
    2188           1 :             double const a21(frand());
    2189           1 :             double const a30(frand());
    2190           1 :             double const a31(frand());
    2191           1 :             double const a40(frand());
    2192           1 :             double const a41(frand());
    2193           1 :             double const a50(frand());
    2194           1 :             double const a51(frand());
    2195             : 
    2196           1 :             a[0][0] = a00;
    2197           1 :             a[0][1] = a01;
    2198           1 :             a[1][0] = a10;
    2199           1 :             a[1][1] = a11;
    2200           1 :             a[2][0] = a20;
    2201           1 :             a[2][1] = a21;
    2202           1 :             a[3][0] = a30;
    2203           1 :             a[3][1] = a31;
    2204           1 :             a[4][0] = a40;
    2205           1 :             a[4][1] = a41;
    2206           1 :             a[5][0] = a50;
    2207           1 :             a[5][1] = a51;
    2208             : 
    2209           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2210           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2211           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2212           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2213           1 :             CATCH_REQUIRE(a[2][0] == a20);
    2214           1 :             CATCH_REQUIRE(a[2][1] == a21);
    2215           1 :             CATCH_REQUIRE(a[3][0] == a30);
    2216           1 :             CATCH_REQUIRE(a[3][1] == a31);
    2217           1 :             CATCH_REQUIRE(a[4][0] == a40);
    2218           1 :             CATCH_REQUIRE(a[4][1] == a41);
    2219           1 :             CATCH_REQUIRE(a[5][0] == a50);
    2220           1 :             CATCH_REQUIRE(a[5][1] == a51);
    2221             : 
    2222           2 :             snapdev::matrix<double> t = a.transpose();
    2223             : 
    2224           1 :             CATCH_REQUIRE(t.rows() == 2);
    2225           1 :             CATCH_REQUIRE(t.columns() == 6);
    2226             : 
    2227             :             // original
    2228             :             // a00 a01
    2229             :             // a10 a11
    2230             :             // a20 a21
    2231             :             // a30 a31
    2232             :             // a40 a41
    2233             :             // a50 a51
    2234             :             //
    2235             :             // transposed
    2236             :             // a00 a10 a20 a30 a40 a50
    2237             :             // a01 a11 a21 a31 a41 a51
    2238             :             //
    2239             : 
    2240           1 :             CATCH_REQUIRE(t[0][0] == a00);
    2241           1 :             CATCH_REQUIRE(t[0][1] == a10);
    2242           1 :             CATCH_REQUIRE(t[0][2] == a20);
    2243           1 :             CATCH_REQUIRE(t[0][3] == a30);
    2244           1 :             CATCH_REQUIRE(t[0][4] == a40);
    2245           1 :             CATCH_REQUIRE(t[0][5] == a50);
    2246           1 :             CATCH_REQUIRE(t[1][0] == a01);
    2247           1 :             CATCH_REQUIRE(t[1][1] == a11);
    2248           1 :             CATCH_REQUIRE(t[1][2] == a21);
    2249           1 :             CATCH_REQUIRE(t[1][3] == a31);
    2250           1 :             CATCH_REQUIRE(t[1][4] == a41);
    2251           1 :             CATCH_REQUIRE(t[1][5] == a51);
    2252             :         }
    2253             :         CATCH_END_SECTION()
    2254             : 
    2255          10 :         CATCH_START_SECTION("matrix: adjugate 2x2")
    2256             :         {
    2257             :             // setup A
    2258             :             //
    2259           2 :             snapdev::matrix<double> a(2, 2);
    2260             : 
    2261           1 :             CATCH_REQUIRE(a.rows() == 2);
    2262           1 :             CATCH_REQUIRE(a.columns() == 2);
    2263             : 
    2264           1 :             double const a00(frand());
    2265           1 :             double const a01(frand());
    2266           1 :             double const a10(frand());
    2267           1 :             double const a11(frand());
    2268             : 
    2269           1 :             a[0][0] = a00;
    2270           1 :             a[0][1] = a01;
    2271           1 :             a[1][0] = a10;
    2272           1 :             a[1][1] = a11;
    2273             : 
    2274           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2275           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2276           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2277           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2278             : 
    2279           2 :             snapdev::matrix<double> m = a.adjugate();
    2280             : 
    2281           1 :             CATCH_REQUIRE(m.rows() == 2);
    2282           1 :             CATCH_REQUIRE(m.columns() == 2);
    2283             : 
    2284           1 :             CATCH_REQUIRE(m[0][0] ==  a11);
    2285           1 :             CATCH_REQUIRE(m[0][1] == -a01);
    2286           1 :             CATCH_REQUIRE(m[1][0] == -a10);
    2287           1 :             CATCH_REQUIRE(m[1][1] ==  a00);
    2288             :         }
    2289             :         CATCH_END_SECTION()
    2290             : 
    2291          10 :         CATCH_START_SECTION("matrix: adjugate 3x3")
    2292             :         {
    2293             :             // setup A
    2294             :             //
    2295           2 :             snapdev::matrix<double> a(3, 3);
    2296             : 
    2297           1 :             CATCH_REQUIRE(a.rows() == 3);
    2298           1 :             CATCH_REQUIRE(a.columns() == 3);
    2299             : 
    2300           1 :             double const a00(frand());
    2301           1 :             double const a01(frand());
    2302           1 :             double const a02(frand());
    2303           1 :             double const a10(frand());
    2304           1 :             double const a11(frand());
    2305           1 :             double const a12(frand());
    2306           1 :             double const a20(frand());
    2307           1 :             double const a21(frand());
    2308           1 :             double const a22(frand());
    2309             : 
    2310           1 :             a[0][0] = a00;
    2311           1 :             a[0][1] = a01;
    2312           1 :             a[0][2] = a02;
    2313           1 :             a[1][0] = a10;
    2314           1 :             a[1][1] = a11;
    2315           1 :             a[1][2] = a12;
    2316           1 :             a[2][0] = a20;
    2317           1 :             a[2][1] = a21;
    2318           1 :             a[2][2] = a22;
    2319             : 
    2320           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2321           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2322           1 :             CATCH_REQUIRE(a[0][2] == a02);
    2323           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2324           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2325           1 :             CATCH_REQUIRE(a[1][2] == a12);
    2326           1 :             CATCH_REQUIRE(a[2][0] == a20);
    2327           1 :             CATCH_REQUIRE(a[2][1] == a21);
    2328           1 :             CATCH_REQUIRE(a[2][2] == a22);
    2329             : 
    2330           2 :             snapdev::matrix<double> m = a.adjugate();
    2331             : 
    2332           1 :             CATCH_REQUIRE(m.rows() == 3);
    2333           1 :             CATCH_REQUIRE(m.columns() == 3);
    2334             : 
    2335             :             {
    2336           2 :                 snapdev::matrix<double> p(a.minor_matrix(0, 0));
    2337           1 :                 double const e(m[0][0] - p.determinant());
    2338           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2339             :             }
    2340             : 
    2341             :             {
    2342           2 :                 snapdev::matrix<double> p(a.minor_matrix(0, 1));
    2343           1 :                 double const e(m[1][0] + p.determinant());
    2344           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2345             :             }
    2346             : 
    2347             :             {
    2348           2 :                 snapdev::matrix<double> p(a.minor_matrix(0, 2));
    2349           1 :                 double const e(m[2][0] - p.determinant());
    2350           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2351             :             }
    2352             : 
    2353             :             {
    2354           2 :                 snapdev::matrix<double> p(a.minor_matrix(1, 0));
    2355           1 :                 double const e(m[0][1] + p.determinant());
    2356           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2357             :             }
    2358             : 
    2359             :             {
    2360           2 :                 snapdev::matrix<double> p(a.minor_matrix(1, 1));
    2361           1 :                 double const e(m[1][1] - p.determinant());
    2362           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2363             :             }
    2364             : 
    2365             :             {
    2366           2 :                 snapdev::matrix<double> p(a.minor_matrix(1, 2));
    2367           1 :                 double const e(m[2][1] + p.determinant());
    2368           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2369             :             }
    2370             : 
    2371             :             {
    2372           2 :                 snapdev::matrix<double> p(a.minor_matrix(2, 0));
    2373           1 :                 double const e(m[0][2] - p.determinant());
    2374           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2375             :             }
    2376             : 
    2377             :             {
    2378           2 :                 snapdev::matrix<double> p(a.minor_matrix(2, 1));
    2379           1 :                 double const e(m[1][2] + p.determinant());
    2380           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2381             :             }
    2382             : 
    2383             :             {
    2384           2 :                 snapdev::matrix<double> p(a.minor_matrix(2, 2));
    2385           1 :                 double const e(m[2][2] - p.determinant());
    2386           1 :                 CATCH_REQUIRE(fabs(e) < 0.0001);
    2387             :             }
    2388             :         }
    2389             :         CATCH_END_SECTION()
    2390             :     }
    2391             : 
    2392             :     // calculations such as the determinant
    2393             :     //
    2394          16 :     CATCH_GIVEN("calculations")
    2395             :     {
    2396           6 :         CATCH_START_SECTION("matrix: 2x2 determinant")
    2397             :         {
    2398             :             // setup A
    2399             :             //
    2400           2 :             snapdev::matrix<double> a(2, 2);
    2401             : 
    2402           1 :             CATCH_REQUIRE(a.rows() == 2);
    2403           1 :             CATCH_REQUIRE(a.columns() == 2);
    2404             : 
    2405           1 :             double const a00(frand());
    2406           1 :             double const a01(frand());
    2407           1 :             double const a10(frand());
    2408           1 :             double const a11(frand());
    2409             : 
    2410           1 :             a[0][0] = a00;
    2411           1 :             a[0][1] = a01;
    2412           1 :             a[1][0] = a10;
    2413           1 :             a[1][1] = a11;
    2414             : 
    2415           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2416           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2417           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2418           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2419             : 
    2420           1 :             double d(a.determinant());
    2421             : 
    2422             :             //return f_vector[0 + 0 * 2] * f_vector[1 + 1 * 2]
    2423             :             //     - f_vector[1 + 0 * 2] * f_vector[0 + 1 * 2];
    2424           1 :             double e = a00 * a11 - a10 * a01;
    2425             : 
    2426           1 :             CATCH_REQUIRE(fabs(d - e) < 0.0001);
    2427             :         }
    2428             :         CATCH_END_SECTION()
    2429             : 
    2430           6 :         CATCH_START_SECTION("matrix: 3x3 determinant with specific data")
    2431             :         {
    2432             :             // setup A
    2433             :             //
    2434           2 :             snapdev::matrix<double> a(3, 3);
    2435             : 
    2436           1 :             CATCH_REQUIRE(a.rows() == 3);
    2437           1 :             CATCH_REQUIRE(a.columns() == 3);
    2438             : 
    2439           1 :             a[0][0] =  5.0;
    2440           1 :             a[0][1] = -2.0;
    2441           1 :             a[0][2] =  1.0;
    2442           1 :             a[1][0] =  0.0;
    2443           1 :             a[1][1] =  3.0;
    2444           1 :             a[1][2] = -1.0;
    2445           1 :             a[2][0] =  2.0;
    2446           1 :             a[2][1] =  0.0;
    2447           1 :             a[2][2] =  7.0;
    2448             : 
    2449           1 :             double d(a.determinant());
    2450             : 
    2451             :             // we know that the exact answer for this one is 103
    2452             :             //
    2453           1 :             CATCH_REQUIRE(fabs(d - 103.0) < 0.0001);
    2454             :         }
    2455             :         CATCH_END_SECTION()
    2456             : 
    2457           6 :         CATCH_START_SECTION("matrix: 3x3 determinant with random data")
    2458             :         {
    2459             :             // verify 10 times (with different values)
    2460             :             //
    2461          11 :             for(int repeat(0); repeat < 10; ++repeat)
    2462             :             {
    2463             :                 // setup A
    2464             :                 //
    2465          20 :                 snapdev::matrix<double> a(3, 3);
    2466             : 
    2467          10 :                 CATCH_REQUIRE(a.rows() == 3);
    2468          10 :                 CATCH_REQUIRE(a.columns() == 3);
    2469             : 
    2470          10 :                 double const a00(frand());
    2471          10 :                 double const a01(frand());
    2472          10 :                 double const a02(frand());
    2473          10 :                 double const a10(frand());
    2474          10 :                 double const a11(frand());
    2475          10 :                 double const a12(frand());
    2476          10 :                 double const a20(frand());
    2477          10 :                 double const a21(frand());
    2478          10 :                 double const a22(frand());
    2479             : 
    2480          10 :                 a[0][0] = a00;
    2481          10 :                 a[0][1] = a01;
    2482          10 :                 a[0][2] = a02;
    2483          10 :                 a[1][0] = a10;
    2484          10 :                 a[1][1] = a11;
    2485          10 :                 a[1][2] = a12;
    2486          10 :                 a[2][0] = a20;
    2487          10 :                 a[2][1] = a21;
    2488          10 :                 a[2][2] = a22;
    2489             : 
    2490          10 :                 CATCH_REQUIRE(a[0][0] == a00);
    2491          10 :                 CATCH_REQUIRE(a[0][1] == a01);
    2492          10 :                 CATCH_REQUIRE(a[0][2] == a02);
    2493          10 :                 CATCH_REQUIRE(a[1][0] == a10);
    2494          10 :                 CATCH_REQUIRE(a[1][1] == a11);
    2495          10 :                 CATCH_REQUIRE(a[1][2] == a12);
    2496          10 :                 CATCH_REQUIRE(a[2][0] == a20);
    2497          10 :                 CATCH_REQUIRE(a[2][1] == a21);
    2498          10 :                 CATCH_REQUIRE(a[2][2] == a22);
    2499             : 
    2500          10 :                 double d(a.determinant());
    2501             : 
    2502             :                 //return f_vector[0 + 0 * 2] * f_vector[1 + 1 * 2]
    2503             :                 //     - f_vector[1 + 0 * 2] * f_vector[0 + 1 * 2];
    2504          20 :                 double e = a00 * a11 * a22 + a01 * a12 * a20 + a02 * a10 * a21
    2505          10 :                          - a00 * a12 * a21 - a01 * a10 * a22 - a02 * a11 * a20;
    2506             : 
    2507          10 :                 CATCH_REQUIRE(fabs(d - e) < 0.0001);
    2508             :             }
    2509             :         }
    2510             :         CATCH_END_SECTION()
    2511             :     }
    2512           8 : }
    2513             : 
    2514             : 
    2515          10 : CATCH_TEST_CASE("matrix_multiplicative", "[matrix]")
    2516             : {
    2517             :     // create two random 4x4 matrices and make sure the add works
    2518             :     //
    2519          16 :     CATCH_GIVEN("mul")
    2520             :     {
    2521           8 :         CATCH_START_SECTION("matrix: b=a*<scalar>")
    2522             :         {
    2523             :             // setup A
    2524             :             //
    2525           2 :             snapdev::matrix<double> a(4, 4);
    2526             : 
    2527           1 :             CATCH_REQUIRE(a.rows() == 4);
    2528           1 :             CATCH_REQUIRE(a.columns() == 4);
    2529             : 
    2530           1 :             double const a00(frand());
    2531           1 :             double const a01(frand());
    2532           1 :             double const a02(frand());
    2533           1 :             double const a03(frand());
    2534           1 :             double const a10(frand());
    2535           1 :             double const a11(frand());
    2536           1 :             double const a12(frand());
    2537           1 :             double const a13(frand());
    2538           1 :             double const a20(frand());
    2539           1 :             double const a21(frand());
    2540           1 :             double const a22(frand());
    2541           1 :             double const a23(frand());
    2542           1 :             double const a30(frand());
    2543           1 :             double const a31(frand());
    2544           1 :             double const a32(frand());
    2545           1 :             double const a33(frand());
    2546             : 
    2547           1 :             a[0][0] = a00;
    2548           1 :             a[0][1] = a01;
    2549           1 :             a[0][2] = a02;
    2550           1 :             a[0][3] = a03;
    2551           1 :             a[1][0] = a10;
    2552           1 :             a[1][1] = a11;
    2553           1 :             a[1][2] = a12;
    2554           1 :             a[1][3] = a13;
    2555           1 :             a[2][0] = a20;
    2556           1 :             a[2][1] = a21;
    2557           1 :             a[2][2] = a22;
    2558           1 :             a[2][3] = a23;
    2559           1 :             a[3][0] = a30;
    2560           1 :             a[3][1] = a31;
    2561           1 :             a[3][2] = a32;
    2562           1 :             a[3][3] = a33;
    2563             : 
    2564           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2565           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2566           1 :             CATCH_REQUIRE(a[0][2] == a02);
    2567           1 :             CATCH_REQUIRE(a[0][3] == a03);
    2568           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2569           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2570           1 :             CATCH_REQUIRE(a[1][2] == a12);
    2571           1 :             CATCH_REQUIRE(a[1][3] == a13);
    2572           1 :             CATCH_REQUIRE(a[2][0] == a20);
    2573           1 :             CATCH_REQUIRE(a[2][1] == a21);
    2574           1 :             CATCH_REQUIRE(a[2][2] == a22);
    2575           1 :             CATCH_REQUIRE(a[2][3] == a23);
    2576           1 :             CATCH_REQUIRE(a[3][0] == a30);
    2577           1 :             CATCH_REQUIRE(a[3][1] == a31);
    2578           1 :             CATCH_REQUIRE(a[3][2] == a32);
    2579           1 :             CATCH_REQUIRE(a[3][3] == a33);
    2580             : 
    2581             :             // setup B
    2582             :             //
    2583           2 :             snapdev::matrix<double> b(4, 4);
    2584             : 
    2585           1 :             CATCH_REQUIRE(b.rows() == 4);
    2586           1 :             CATCH_REQUIRE(b.columns() == 4);
    2587             : 
    2588           1 :             double const b00(frand());
    2589           1 :             double const b01(frand());
    2590           1 :             double const b02(frand());
    2591           1 :             double const b03(frand());
    2592           1 :             double const b10(frand());
    2593           1 :             double const b11(frand());
    2594           1 :             double const b12(frand());
    2595           1 :             double const b13(frand());
    2596           1 :             double const b20(frand());
    2597           1 :             double const b21(frand());
    2598           1 :             double const b22(frand());
    2599           1 :             double const b23(frand());
    2600           1 :             double const b30(frand());
    2601           1 :             double const b31(frand());
    2602           1 :             double const b32(frand());
    2603           1 :             double const b33(frand());
    2604             : 
    2605           1 :             b[0][0] = b00;
    2606           1 :             b[0][1] = b01;
    2607           1 :             b[0][2] = b02;
    2608           1 :             b[0][3] = b03;
    2609           1 :             b[1][0] = b10;
    2610           1 :             b[1][1] = b11;
    2611           1 :             b[1][2] = b12;
    2612           1 :             b[1][3] = b13;
    2613           1 :             b[2][0] = b20;
    2614           1 :             b[2][1] = b21;
    2615           1 :             b[2][2] = b22;
    2616           1 :             b[2][3] = b23;
    2617           1 :             b[3][0] = b30;
    2618           1 :             b[3][1] = b31;
    2619           1 :             b[3][2] = b32;
    2620           1 :             b[3][3] = b33;
    2621             : 
    2622           1 :             CATCH_REQUIRE(b[0][0] == b00);
    2623           1 :             CATCH_REQUIRE(b[0][1] == b01);
    2624           1 :             CATCH_REQUIRE(b[0][2] == b02);
    2625           1 :             CATCH_REQUIRE(b[0][3] == b03);
    2626           1 :             CATCH_REQUIRE(b[1][0] == b10);
    2627           1 :             CATCH_REQUIRE(b[1][1] == b11);
    2628           1 :             CATCH_REQUIRE(b[1][2] == b12);
    2629           1 :             CATCH_REQUIRE(b[1][3] == b13);
    2630           1 :             CATCH_REQUIRE(b[2][0] == b20);
    2631           1 :             CATCH_REQUIRE(b[2][1] == b21);
    2632           1 :             CATCH_REQUIRE(b[2][2] == b22);
    2633           1 :             CATCH_REQUIRE(b[2][3] == b23);
    2634           1 :             CATCH_REQUIRE(b[3][0] == b30);
    2635           1 :             CATCH_REQUIRE(b[3][1] == b31);
    2636           1 :             CATCH_REQUIRE(b[3][2] == b32);
    2637           1 :             CATCH_REQUIRE(b[3][3] == b33);
    2638             : 
    2639             :             // create a scalar for our test
    2640             :             //
    2641           1 :             double scalar(frand());
    2642           1 :             while(scalar == 0.0 || scalar == 1.0)
    2643             :             {
    2644           0 :                 scalar = frand();
    2645             :             }
    2646             : 
    2647             :             // run operation B = A + <scalar>
    2648             :             //
    2649           1 :             b = a * scalar;
    2650             : 
    2651           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2652           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2653           1 :             CATCH_REQUIRE(a[0][2] == a02);
    2654           1 :             CATCH_REQUIRE(a[0][3] == a03);
    2655           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2656           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2657           1 :             CATCH_REQUIRE(a[1][2] == a12);
    2658           1 :             CATCH_REQUIRE(a[1][3] == a13);
    2659           1 :             CATCH_REQUIRE(a[2][0] == a20);
    2660           1 :             CATCH_REQUIRE(a[2][1] == a21);
    2661           1 :             CATCH_REQUIRE(a[2][2] == a22);
    2662           1 :             CATCH_REQUIRE(a[2][3] == a23);
    2663           1 :             CATCH_REQUIRE(a[3][0] == a30);
    2664           1 :             CATCH_REQUIRE(a[3][1] == a31);
    2665           1 :             CATCH_REQUIRE(a[3][2] == a32);
    2666           1 :             CATCH_REQUIRE(a[3][3] == a33);
    2667             : 
    2668             :             // this can't fail because we ensure scalar != 0 or 1
    2669             :             //
    2670           1 :             CATCH_REQUIRE_FALSE(b[0][0] == b00);
    2671           1 :             CATCH_REQUIRE_FALSE(b[0][1] == b01);
    2672           1 :             CATCH_REQUIRE_FALSE(b[0][2] == b02);
    2673           1 :             CATCH_REQUIRE_FALSE(b[0][3] == b03);
    2674           1 :             CATCH_REQUIRE_FALSE(b[1][0] == b10);
    2675           1 :             CATCH_REQUIRE_FALSE(b[1][1] == b11);
    2676           1 :             CATCH_REQUIRE_FALSE(b[1][2] == b12);
    2677           1 :             CATCH_REQUIRE_FALSE(b[1][3] == b13);
    2678           1 :             CATCH_REQUIRE_FALSE(b[2][0] == b20);
    2679           1 :             CATCH_REQUIRE_FALSE(b[2][1] == b21);
    2680           1 :             CATCH_REQUIRE_FALSE(b[2][2] == b22);
    2681           1 :             CATCH_REQUIRE_FALSE(b[2][3] == b23);
    2682           1 :             CATCH_REQUIRE_FALSE(b[3][0] == b30);
    2683           1 :             CATCH_REQUIRE_FALSE(b[3][1] == b31);
    2684           1 :             CATCH_REQUIRE_FALSE(b[3][2] == b32);
    2685           1 :             CATCH_REQUIRE_FALSE(b[3][3] == b33);
    2686             : 
    2687           1 :             CATCH_REQUIRE(b[0][0] == a00 * scalar);
    2688           1 :             CATCH_REQUIRE(b[0][1] == a01 * scalar);
    2689           1 :             CATCH_REQUIRE(b[0][2] == a02 * scalar);
    2690           1 :             CATCH_REQUIRE(b[0][3] == a03 * scalar);
    2691           1 :             CATCH_REQUIRE(b[1][0] == a10 * scalar);
    2692           1 :             CATCH_REQUIRE(b[1][1] == a11 * scalar);
    2693           1 :             CATCH_REQUIRE(b[1][2] == a12 * scalar);
    2694           1 :             CATCH_REQUIRE(b[1][3] == a13 * scalar);
    2695           1 :             CATCH_REQUIRE(b[2][0] == a20 * scalar);
    2696           1 :             CATCH_REQUIRE(b[2][1] == a21 * scalar);
    2697           1 :             CATCH_REQUIRE(b[2][2] == a22 * scalar);
    2698           1 :             CATCH_REQUIRE(b[2][3] == a23 * scalar);
    2699           1 :             CATCH_REQUIRE(b[3][0] == a30 * scalar);
    2700           1 :             CATCH_REQUIRE(b[3][1] == a31 * scalar);
    2701           1 :             CATCH_REQUIRE(b[3][2] == a32 * scalar);
    2702           1 :             CATCH_REQUIRE(b[3][3] == a33 * scalar);
    2703             :         }
    2704             :         CATCH_END_SECTION()
    2705             : 
    2706           8 :         CATCH_START_SECTION("matrix: a*=<scalar>")
    2707             :         {
    2708             :             // setup A
    2709             :             //
    2710           2 :             snapdev::matrix<double> a(4, 4);
    2711             : 
    2712           1 :             CATCH_REQUIRE(a.rows() == 4);
    2713           1 :             CATCH_REQUIRE(a.columns() == 4);
    2714             : 
    2715           1 :             double const a00(frand());
    2716           1 :             double const a01(frand());
    2717           1 :             double const a02(frand());
    2718           1 :             double const a03(frand());
    2719           1 :             double const a10(frand());
    2720           1 :             double const a11(frand());
    2721           1 :             double const a12(frand());
    2722           1 :             double const a13(frand());
    2723           1 :             double const a20(frand());
    2724           1 :             double const a21(frand());
    2725           1 :             double const a22(frand());
    2726           1 :             double const a23(frand());
    2727           1 :             double const a30(frand());
    2728           1 :             double const a31(frand());
    2729           1 :             double const a32(frand());
    2730           1 :             double const a33(frand());
    2731             : 
    2732           1 :             a[0][0] = a00;
    2733           1 :             a[0][1] = a01;
    2734           1 :             a[0][2] = a02;
    2735           1 :             a[0][3] = a03;
    2736           1 :             a[1][0] = a10;
    2737           1 :             a[1][1] = a11;
    2738           1 :             a[1][2] = a12;
    2739           1 :             a[1][3] = a13;
    2740           1 :             a[2][0] = a20;
    2741           1 :             a[2][1] = a21;
    2742           1 :             a[2][2] = a22;
    2743           1 :             a[2][3] = a23;
    2744           1 :             a[3][0] = a30;
    2745           1 :             a[3][1] = a31;
    2746           1 :             a[3][2] = a32;
    2747           1 :             a[3][3] = a33;
    2748             : 
    2749           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2750           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2751           1 :             CATCH_REQUIRE(a[0][2] == a02);
    2752           1 :             CATCH_REQUIRE(a[0][3] == a03);
    2753           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2754           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2755           1 :             CATCH_REQUIRE(a[1][2] == a12);
    2756           1 :             CATCH_REQUIRE(a[1][3] == a13);
    2757           1 :             CATCH_REQUIRE(a[2][0] == a20);
    2758           1 :             CATCH_REQUIRE(a[2][1] == a21);
    2759           1 :             CATCH_REQUIRE(a[2][2] == a22);
    2760           1 :             CATCH_REQUIRE(a[2][3] == a23);
    2761           1 :             CATCH_REQUIRE(a[3][0] == a30);
    2762           1 :             CATCH_REQUIRE(a[3][1] == a31);
    2763           1 :             CATCH_REQUIRE(a[3][2] == a32);
    2764           1 :             CATCH_REQUIRE(a[3][3] == a33);
    2765             : 
    2766             :             // create a scalar for our test
    2767             :             //
    2768           1 :             double scalar(frand());
    2769           1 :             while(scalar == 0.0 || scalar == 1.0)
    2770             :             {
    2771           0 :                 scalar = frand();
    2772             :             }
    2773             : 
    2774             :             // run operation A *= <scalar>
    2775             :             //
    2776           1 :             a *= scalar;
    2777             : 
    2778             :             // this can't fail because we ensure scalar != 0 or 1
    2779             :             //
    2780           1 :             CATCH_REQUIRE_FALSE(a[0][0] == a00);
    2781           1 :             CATCH_REQUIRE_FALSE(a[0][1] == a01);
    2782           1 :             CATCH_REQUIRE_FALSE(a[0][2] == a02);
    2783           1 :             CATCH_REQUIRE_FALSE(a[0][3] == a03);
    2784           1 :             CATCH_REQUIRE_FALSE(a[1][0] == a10);
    2785           1 :             CATCH_REQUIRE_FALSE(a[1][1] == a11);
    2786           1 :             CATCH_REQUIRE_FALSE(a[1][2] == a12);
    2787           1 :             CATCH_REQUIRE_FALSE(a[1][3] == a13);
    2788           1 :             CATCH_REQUIRE_FALSE(a[2][0] == a20);
    2789           1 :             CATCH_REQUIRE_FALSE(a[2][1] == a21);
    2790           1 :             CATCH_REQUIRE_FALSE(a[2][2] == a22);
    2791           1 :             CATCH_REQUIRE_FALSE(a[2][3] == a23);
    2792           1 :             CATCH_REQUIRE_FALSE(a[3][0] == a30);
    2793           1 :             CATCH_REQUIRE_FALSE(a[3][1] == a31);
    2794           1 :             CATCH_REQUIRE_FALSE(a[3][2] == a32);
    2795           1 :             CATCH_REQUIRE_FALSE(a[3][3] == a33);
    2796             : 
    2797           1 :             CATCH_REQUIRE(a[0][0] == a00 * scalar);
    2798           1 :             CATCH_REQUIRE(a[0][1] == a01 * scalar);
    2799           1 :             CATCH_REQUIRE(a[0][2] == a02 * scalar);
    2800           1 :             CATCH_REQUIRE(a[0][3] == a03 * scalar);
    2801           1 :             CATCH_REQUIRE(a[1][0] == a10 * scalar);
    2802           1 :             CATCH_REQUIRE(a[1][1] == a11 * scalar);
    2803           1 :             CATCH_REQUIRE(a[1][2] == a12 * scalar);
    2804           1 :             CATCH_REQUIRE(a[1][3] == a13 * scalar);
    2805           1 :             CATCH_REQUIRE(a[2][0] == a20 * scalar);
    2806           1 :             CATCH_REQUIRE(a[2][1] == a21 * scalar);
    2807           1 :             CATCH_REQUIRE(a[2][2] == a22 * scalar);
    2808           1 :             CATCH_REQUIRE(a[2][3] == a23 * scalar);
    2809           1 :             CATCH_REQUIRE(a[3][0] == a30 * scalar);
    2810           1 :             CATCH_REQUIRE(a[3][1] == a31 * scalar);
    2811           1 :             CATCH_REQUIRE(a[3][2] == a32 * scalar);
    2812           1 :             CATCH_REQUIRE(a[3][3] == a33 * scalar);
    2813             :         }
    2814             :         CATCH_END_SECTION()
    2815             : 
    2816           8 :         CATCH_START_SECTION("matrix: c=a*b")
    2817             :         {
    2818             :             // setup A
    2819             :             //
    2820           2 :             snapdev::matrix<double> a(4, 4);
    2821             : 
    2822           1 :             CATCH_REQUIRE(a.rows() == 4);
    2823           1 :             CATCH_REQUIRE(a.columns() == 4);
    2824             : 
    2825           1 :             double const a00(frand());
    2826           1 :             double const a01(frand());
    2827           1 :             double const a02(frand());
    2828           1 :             double const a03(frand());
    2829           1 :             double const a10(frand());
    2830           1 :             double const a11(frand());
    2831           1 :             double const a12(frand());
    2832           1 :             double const a13(frand());
    2833           1 :             double const a20(frand());
    2834           1 :             double const a21(frand());
    2835           1 :             double const a22(frand());
    2836           1 :             double const a23(frand());
    2837           1 :             double const a30(frand());
    2838           1 :             double const a31(frand());
    2839           1 :             double const a32(frand());
    2840           1 :             double const a33(frand());
    2841             : 
    2842           1 :             a[0][0] = a00;
    2843           1 :             a[0][1] = a01;
    2844           1 :             a[0][2] = a02;
    2845           1 :             a[0][3] = a03;
    2846           1 :             a[1][0] = a10;
    2847           1 :             a[1][1] = a11;
    2848           1 :             a[1][2] = a12;
    2849           1 :             a[1][3] = a13;
    2850           1 :             a[2][0] = a20;
    2851           1 :             a[2][1] = a21;
    2852           1 :             a[2][2] = a22;
    2853           1 :             a[2][3] = a23;
    2854           1 :             a[3][0] = a30;
    2855           1 :             a[3][1] = a31;
    2856           1 :             a[3][2] = a32;
    2857           1 :             a[3][3] = a33;
    2858             : 
    2859           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2860           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2861           1 :             CATCH_REQUIRE(a[0][2] == a02);
    2862           1 :             CATCH_REQUIRE(a[0][3] == a03);
    2863           1 :             CATCH_REQUIRE(a[1][0] == a10);
    2864           1 :             CATCH_REQUIRE(a[1][1] == a11);
    2865           1 :             CATCH_REQUIRE(a[1][2] == a12);
    2866           1 :             CATCH_REQUIRE(a[1][3] == a13);
    2867           1 :             CATCH_REQUIRE(a[2][0] == a20);
    2868           1 :             CATCH_REQUIRE(a[2][1] == a21);
    2869           1 :             CATCH_REQUIRE(a[2][2] == a22);
    2870           1 :             CATCH_REQUIRE(a[2][3] == a23);
    2871           1 :             CATCH_REQUIRE(a[3][0] == a30);
    2872           1 :             CATCH_REQUIRE(a[3][1] == a31);
    2873           1 :             CATCH_REQUIRE(a[3][2] == a32);
    2874           1 :             CATCH_REQUIRE(a[3][3] == a33);
    2875             : 
    2876             :             // setup B
    2877             :             //
    2878           2 :             snapdev::matrix<double> b(4, 4);
    2879             : 
    2880           1 :             CATCH_REQUIRE(b.rows() == 4);
    2881           1 :             CATCH_REQUIRE(b.columns() == 4);
    2882             : 
    2883           1 :             double const b00(frand());
    2884           1 :             double const b01(frand());
    2885           1 :             double const b02(frand());
    2886           1 :             double const b03(frand());
    2887           1 :             double const b10(frand());
    2888           1 :             double const b11(frand());
    2889           1 :             double const b12(frand());
    2890           1 :             double const b13(frand());
    2891           1 :             double const b20(frand());
    2892           1 :             double const b21(frand());
    2893           1 :             double const b22(frand());
    2894           1 :             double const b23(frand());
    2895           1 :             double const b30(frand());
    2896           1 :             double const b31(frand());
    2897           1 :             double const b32(frand());
    2898           1 :             double const b33(frand());
    2899             : 
    2900           1 :             b[0][0] = b00;
    2901           1 :             b[0][1] = b01;
    2902           1 :             b[0][2] = b02;
    2903           1 :             b[0][3] = b03;
    2904           1 :             b[1][0] = b10;
    2905           1 :             b[1][1] = b11;
    2906           1 :             b[1][2] = b12;
    2907           1 :             b[1][3] = b13;
    2908           1 :             b[2][0] = b20;
    2909           1 :             b[2][1] = b21;
    2910           1 :             b[2][2] = b22;
    2911           1 :             b[2][3] = b23;
    2912           1 :             b[3][0] = b30;
    2913           1 :             b[3][1] = b31;
    2914           1 :             b[3][2] = b32;
    2915           1 :             b[3][3] = b33;
    2916             : 
    2917           1 :             CATCH_REQUIRE(b[0][0] == b00);
    2918           1 :             CATCH_REQUIRE(b[0][1] == b01);
    2919           1 :             CATCH_REQUIRE(b[0][2] == b02);
    2920           1 :             CATCH_REQUIRE(b[0][3] == b03);
    2921           1 :             CATCH_REQUIRE(b[1][0] == b10);
    2922           1 :             CATCH_REQUIRE(b[1][1] == b11);
    2923           1 :             CATCH_REQUIRE(b[1][2] == b12);
    2924           1 :             CATCH_REQUIRE(b[1][3] == b13);
    2925           1 :             CATCH_REQUIRE(b[2][0] == b20);
    2926           1 :             CATCH_REQUIRE(b[2][1] == b21);
    2927           1 :             CATCH_REQUIRE(b[2][2] == b22);
    2928           1 :             CATCH_REQUIRE(b[2][3] == b23);
    2929           1 :             CATCH_REQUIRE(b[3][0] == b30);
    2930           1 :             CATCH_REQUIRE(b[3][1] == b31);
    2931           1 :             CATCH_REQUIRE(b[3][2] == b32);
    2932           1 :             CATCH_REQUIRE(b[3][3] == b33);
    2933             : 
    2934             :             // setup C
    2935             :             //
    2936           2 :             snapdev::matrix<double> c(4, 4);
    2937             : 
    2938           1 :             CATCH_REQUIRE(c.rows() == 4);
    2939           1 :             CATCH_REQUIRE(c.columns() == 4);
    2940             : 
    2941           1 :             double const c00(frand());
    2942           1 :             double const c01(frand());
    2943           1 :             double const c02(frand());
    2944           1 :             double const c03(frand());
    2945           1 :             double const c10(frand());
    2946           1 :             double const c11(frand());
    2947           1 :             double const c12(frand());
    2948           1 :             double const c13(frand());
    2949           1 :             double const c20(frand());
    2950           1 :             double const c21(frand());
    2951           1 :             double const c22(frand());
    2952           1 :             double const c23(frand());
    2953           1 :             double const c30(frand());
    2954           1 :             double const c31(frand());
    2955           1 :             double const c32(frand());
    2956           1 :             double const c33(frand());
    2957             : 
    2958           1 :             c[0][0] = c00;
    2959           1 :             c[0][1] = c01;
    2960           1 :             c[0][2] = c02;
    2961           1 :             c[0][3] = c03;
    2962           1 :             c[1][0] = c10;
    2963           1 :             c[1][1] = c11;
    2964           1 :             c[1][2] = c12;
    2965           1 :             c[1][3] = c13;
    2966           1 :             c[2][0] = c20;
    2967           1 :             c[2][1] = c21;
    2968           1 :             c[2][2] = c22;
    2969           1 :             c[2][3] = c23;
    2970           1 :             c[3][0] = c30;
    2971           1 :             c[3][1] = c31;
    2972           1 :             c[3][2] = c32;
    2973           1 :             c[3][3] = c33;
    2974             : 
    2975           1 :             CATCH_REQUIRE(c[0][0] == c00);
    2976           1 :             CATCH_REQUIRE(c[0][1] == c01);
    2977           1 :             CATCH_REQUIRE(c[0][2] == c02);
    2978           1 :             CATCH_REQUIRE(c[0][3] == c03);
    2979           1 :             CATCH_REQUIRE(c[1][0] == c10);
    2980           1 :             CATCH_REQUIRE(c[1][1] == c11);
    2981           1 :             CATCH_REQUIRE(c[1][2] == c12);
    2982           1 :             CATCH_REQUIRE(c[1][3] == c13);
    2983           1 :             CATCH_REQUIRE(c[2][0] == c20);
    2984           1 :             CATCH_REQUIRE(c[2][1] == c21);
    2985           1 :             CATCH_REQUIRE(c[2][2] == c22);
    2986           1 :             CATCH_REQUIRE(c[2][3] == c23);
    2987           1 :             CATCH_REQUIRE(c[3][0] == c30);
    2988           1 :             CATCH_REQUIRE(c[3][1] == c31);
    2989           1 :             CATCH_REQUIRE(c[3][2] == c32);
    2990           1 :             CATCH_REQUIRE(c[3][3] == c33);
    2991             : 
    2992             :             // run operation C = A * B
    2993             :             //
    2994           1 :             c = a * b;
    2995             : 
    2996           1 :             CATCH_REQUIRE(a[0][0] == a00);
    2997           1 :             CATCH_REQUIRE(a[0][1] == a01);
    2998           1 :             CATCH_REQUIRE(a[0][2] == a02);
    2999           1 :             CATCH_REQUIRE(a[0][3] == a03);
    3000           1 :             CATCH_REQUIRE(a[1][0] == a10);
    3001           1 :             CATCH_REQUIRE(a[1][1] == a11);
    3002           1 :             CATCH_REQUIRE(a[1][2] == a12);
    3003           1 :             CATCH_REQUIRE(a[1][3] == a13);
    3004           1 :             CATCH_REQUIRE(a[2][0] == a20);
    3005           1 :             CATCH_REQUIRE(a[2][1] == a21);
    3006           1 :             CATCH_REQUIRE(a[2][2] == a22);
    3007           1 :             CATCH_REQUIRE(a[2][3] == a23);
    3008           1 :             CATCH_REQUIRE(a[3][0] == a30);
    3009           1 :             CATCH_REQUIRE(a[3][1] == a31);
    3010           1 :             CATCH_REQUIRE(a[3][2] == a32);
    3011           1 :             CATCH_REQUIRE(a[3][3] == a33);
    3012             : 
    3013           1 :             CATCH_REQUIRE(b[0][0] == b00);
    3014           1 :             CATCH_REQUIRE(b[0][1] == b01);
    3015           1 :             CATCH_REQUIRE(b[0][2] == b02);
    3016           1 :             CATCH_REQUIRE(b[0][3] == b03);
    3017           1 :             CATCH_REQUIRE(b[1][0] == b10);
    3018           1 :             CATCH_REQUIRE(b[1][1] == b11);
    3019           1 :             CATCH_REQUIRE(b[1][2] == b12);
    3020           1 :             CATCH_REQUIRE(b[1][3] == b13);
    3021           1 :             CATCH_REQUIRE(b[2][0] == b20);
    3022           1 :             CATCH_REQUIRE(b[2][1] == b21);
    3023           1 :             CATCH_REQUIRE(b[2][2] == b22);
    3024           1 :             CATCH_REQUIRE(b[2][3] == b23);
    3025           1 :             CATCH_REQUIRE(b[3][0] == b30);
    3026           1 :             CATCH_REQUIRE(b[3][1] == b31);
    3027           1 :             CATCH_REQUIRE(b[3][2] == b32);
    3028           1 :             CATCH_REQUIRE(b[3][3] == b33);
    3029             : 
    3030             :             // this could fail (albeit rather unlikely)
    3031             :             //CATCH_REQUIRE_FALSE(c[0][0] == c00);
    3032             :             //CATCH_REQUIRE_FALSE(c[0][1] == c01);
    3033             :             //CATCH_REQUIRE_FALSE(c[0][2] == c02);
    3034             :             //CATCH_REQUIRE_FALSE(c[0][3] == c03);
    3035             :             //CATCH_REQUIRE_FALSE(c[1][0] == c10);
    3036             :             //CATCH_REQUIRE_FALSE(c[1][1] == c11);
    3037             :             //CATCH_REQUIRE_FALSE(c[1][2] == c12);
    3038             :             //CATCH_REQUIRE_FALSE(c[1][3] == c13);
    3039             :             //CATCH_REQUIRE_FALSE(c[2][0] == c20);
    3040             :             //CATCH_REQUIRE_FALSE(c[2][1] == c21);
    3041             :             //CATCH_REQUIRE_FALSE(c[2][2] == c22);
    3042             :             //CATCH_REQUIRE_FALSE(c[2][3] == c23);
    3043             :             //CATCH_REQUIRE_FALSE(c[3][0] == c30);
    3044             :             //CATCH_REQUIRE_FALSE(c[3][1] == c31);
    3045             :             //CATCH_REQUIRE_FALSE(c[3][2] == c32);
    3046             :             //CATCH_REQUIRE_FALSE(c[3][3] == c33);
    3047             : 
    3048             : //std::cout << "A = " << a << std::endl;
    3049             : //std::cout << "B = " << b << std::endl;
    3050             : //std::cout << "C = " << c << std::endl;
    3051             : 
    3052           1 :             CATCH_REQUIRE(c[0][0] == a00 * b00 + a01 * b10 + a02 * b20 + a03 * b30);
    3053           1 :             CATCH_REQUIRE(c[0][1] == a00 * b01 + a01 * b11 + a02 * b21 + a03 * b31);
    3054           1 :             CATCH_REQUIRE(c[0][2] == a00 * b02 + a01 * b12 + a02 * b22 + a03 * b32);
    3055           1 :             CATCH_REQUIRE(c[0][3] == a00 * b03 + a01 * b13 + a02 * b23 + a03 * b33);
    3056             : 
    3057           1 :             CATCH_REQUIRE(c[1][0] == a10 * b00 + a11 * b10 + a12 * b20 + a13 * b30);
    3058           1 :             CATCH_REQUIRE(c[1][1] == a10 * b01 + a11 * b11 + a12 * b21 + a13 * b31);
    3059           1 :             CATCH_REQUIRE(c[1][2] == a10 * b02 + a11 * b12 + a12 * b22 + a13 * b32);
    3060           1 :             CATCH_REQUIRE(c[1][3] == a10 * b03 + a11 * b13 + a12 * b23 + a13 * b33);
    3061             : 
    3062           1 :             CATCH_REQUIRE(c[2][0] == a20 * b00 + a21 * b10 + a22 * b20 + a23 * b30);
    3063           1 :             CATCH_REQUIRE(c[2][1] == a20 * b01 + a21 * b11 + a22 * b21 + a23 * b31);
    3064           1 :             CATCH_REQUIRE(c[2][2] == a20 * b02 + a21 * b12 + a22 * b22 + a23 * b32);
    3065           1 :             CATCH_REQUIRE(c[2][3] == a20 * b03 + a21 * b13 + a22 * b23 + a23 * b33);
    3066             : 
    3067           1 :             CATCH_REQUIRE(c[3][0] == a30 * b00 + a31 * b10 + a32 * b20 + a33 * b30);
    3068           1 :             CATCH_REQUIRE(c[3][1] == a30 * b01 + a31 * b11 + a32 * b21 + a33 * b31);
    3069           1 :             CATCH_REQUIRE(c[3][2] == a30 * b02 + a31 * b12 + a32 * b22 + a33 * b32);
    3070           1 :             CATCH_REQUIRE(c[3][3] == a30 * b03 + a31 * b13 + a32 * b23 + a33 * b33);
    3071             : 
    3072             :             // if swapped it would equal this instead
    3073             :             //CATCH_REQUIRE(c[0][0] == a00 * b00 + a10 * b01 + a20 * b02 + a30 * b03);
    3074             :             //CATCH_REQUIRE(c[0][1] == a01 * b00 + a11 * b01 + a21 * b02 + a31 * b03);
    3075             :             //CATCH_REQUIRE(c[0][2] == a02 * b00 + a12 * b01 + a22 * b02 + a32 * b03);
    3076             :             //CATCH_REQUIRE(c[0][3] == a03 * b00 + a13 * b01 + a23 * b02 + a33 * b03);
    3077             : 
    3078             :             //CATCH_REQUIRE(c[1][0] == a00 * b10 + a10 * b11 + a20 * b12 + a30 * b13);
    3079             :             //CATCH_REQUIRE(c[1][1] == a01 * b10 + a11 * b11 + a21 * b12 + a31 * b13);
    3080             :             //CATCH_REQUIRE(c[1][2] == a02 * b10 + a12 * b11 + a22 * b12 + a32 * b13);
    3081             :             //CATCH_REQUIRE(c[1][3] == a03 * b10 + a13 * b11 + a23 * b12 + a33 * b13);
    3082             : 
    3083             :             //CATCH_REQUIRE(c[2][0] == a00 * b20 + a10 * b21 + a20 * b22 + a30 * b23);
    3084             :             //CATCH_REQUIRE(c[2][1] == a01 * b20 + a11 * b21 + a21 * b22 + a31 * b23);
    3085             :             //CATCH_REQUIRE(c[2][2] == a02 * b20 + a12 * b21 + a22 * b22 + a32 * b23);
    3086             :             //CATCH_REQUIRE(c[2][3] == a03 * b20 + a13 * b21 + a23 * b22 + a33 * b23);
    3087             : 
    3088             :             //CATCH_REQUIRE(c[3][0] == a00 * b30 + a10 * b31 + a20 * b32 + a30 * b33);
    3089             :             //CATCH_REQUIRE(c[3][1] == a01 * b30 + a11 * b31 + a21 * b32 + a31 * b33);
    3090             :             //CATCH_REQUIRE(c[3][2] == a02 * b30 + a12 * b31 + a22 * b32 + a32 * b33);
    3091             :             //CATCH_REQUIRE(c[3][3] == a03 * b30 + a13 * b31 + a23 * b32 + a33 * b33);
    3092             :         }
    3093             :         CATCH_END_SECTION()
    3094             : 
    3095           8 :         CATCH_START_SECTION("matrix: a*=b")
    3096             :         {
    3097             :             // setup A
    3098             :             //
    3099           2 :             snapdev::matrix<double> a(4, 4);
    3100             : 
    3101           1 :             CATCH_REQUIRE(a.rows() == 4);
    3102           1 :             CATCH_REQUIRE(a.columns() == 4);
    3103             : 
    3104           1 :             double const a00(frand());
    3105           1 :             double const a01(frand());
    3106           1 :             double const a02(frand());
    3107           1 :             double const a03(frand());
    3108           1 :             double const a10(frand());
    3109           1 :             double const a11(frand());
    3110           1 :             double const a12(frand());
    3111           1 :             double const a13(frand());
    3112           1 :             double const a20(frand());
    3113           1 :             double const a21(frand());
    3114           1 :             double const a22(frand());
    3115           1 :             double const a23(frand());
    3116           1 :             double const a30(frand());
    3117           1 :             double const a31(frand());
    3118           1 :             double const a32(frand());
    3119           1 :             double const a33(frand());
    3120             : 
    3121           1 :             a[0][0] = a00;
    3122           1 :             a[0][1] = a01;
    3123           1 :             a[0][2] = a02;
    3124           1 :             a[0][3] = a03;
    3125           1 :             a[1][0] = a10;
    3126           1 :             a[1][1] = a11;
    3127           1 :             a[1][2] = a12;
    3128           1 :             a[1][3] = a13;
    3129           1 :             a[2][0] = a20;
    3130           1 :             a[2][1] = a21;
    3131           1 :             a[2][2] = a22;
    3132           1 :             a[2][3] = a23;
    3133           1 :             a[3][0] = a30;
    3134           1 :             a[3][1] = a31;
    3135           1 :             a[3][2] = a32;
    3136           1 :             a[3][3] = a33;
    3137             : 
    3138           1 :             CATCH_REQUIRE(a[0][0] == a00);
    3139           1 :             CATCH_REQUIRE(a[0][1] == a01);
    3140           1 :             CATCH_REQUIRE(a[0][2] == a02);
    3141           1 :             CATCH_REQUIRE(a[0][3] == a03);
    3142           1 :             CATCH_REQUIRE(a[1][0] == a10);
    3143           1 :             CATCH_REQUIRE(a[1][1] == a11);
    3144           1 :             CATCH_REQUIRE(a[1][2] == a12);
    3145           1 :             CATCH_REQUIRE(a[1][3] == a13);
    3146           1 :             CATCH_REQUIRE(a[2][0] == a20);
    3147           1 :             CATCH_REQUIRE(a[2][1] == a21);
    3148           1 :             CATCH_REQUIRE(a[2][2] == a22);
    3149           1 :             CATCH_REQUIRE(a[2][3] == a23);
    3150           1 :             CATCH_REQUIRE(a[3][0] == a30);
    3151           1 :             CATCH_REQUIRE(a[3][1] == a31);
    3152           1 :             CATCH_REQUIRE(a[3][2] == a32);
    3153           1 :             CATCH_REQUIRE(a[3][3] == a33);
    3154             : 
    3155             :             // setup B
    3156             :             //
    3157           2 :             snapdev::matrix<double> b(4, 4);
    3158             : 
    3159           1 :             CATCH_REQUIRE(b.rows() == 4);
    3160           1 :             CATCH_REQUIRE(b.columns() == 4);
    3161             : 
    3162           1 :             double const b00(frand());
    3163           1 :             double const b01(frand());
    3164           1 :             double const b02(frand());
    3165           1 :             double const b03(frand());
    3166           1 :             double const b10(frand());
    3167           1 :             double const b11(frand());
    3168           1 :             double const b12(frand());
    3169           1 :             double const b13(frand());
    3170           1 :             double const b20(frand());
    3171           1 :             double const b21(frand());
    3172           1 :             double const b22(frand());
    3173           1 :             double const b23(frand());
    3174           1 :             double const b30(frand());
    3175           1 :             double const b31(frand());
    3176           1 :             double const b32(frand());
    3177           1 :             double const b33(frand());
    3178             : 
    3179           1 :             b[0][0] = b00;
    3180           1 :             b[0][1] = b01;
    3181           1 :             b[0][2] = b02;
    3182           1 :             b[0][3] = b03;
    3183           1 :             b[1][0] = b10;
    3184           1 :             b[1][1] = b11;
    3185           1 :             b[1][2] = b12;
    3186           1 :             b[1][3] = b13;
    3187           1 :             b[2][0] = b20;
    3188           1 :             b[2][1] = b21;
    3189           1 :             b[2][2] = b22;
    3190           1 :             b[2][3] = b23;
    3191           1 :             b[3][0] = b30;
    3192           1 :             b[3][1] = b31;
    3193           1 :             b[3][2] = b32;
    3194           1 :             b[3][3] = b33;
    3195             : 
    3196           1 :             CATCH_REQUIRE(b[0][0] == b00);
    3197           1 :             CATCH_REQUIRE(b[0][1] == b01);
    3198           1 :             CATCH_REQUIRE(b[0][2] == b02);
    3199           1 :             CATCH_REQUIRE(b[0][3] == b03);
    3200           1 :             CATCH_REQUIRE(b[1][0] == b10);
    3201           1 :             CATCH_REQUIRE(b[1][1] == b11);
    3202           1 :             CATCH_REQUIRE(b[1][2] == b12);
    3203           1 :             CATCH_REQUIRE(b[1][3] == b13);
    3204           1 :             CATCH_REQUIRE(b[2][0] == b20);
    3205           1 :             CATCH_REQUIRE(b[2][1] == b21);
    3206           1 :             CATCH_REQUIRE(b[2][2] == b22);
    3207           1 :             CATCH_REQUIRE(b[2][3] == b23);
    3208           1 :             CATCH_REQUIRE(b[3][0] == b30);
    3209           1 :             CATCH_REQUIRE(b[3][1] == b31);
    3210           1 :             CATCH_REQUIRE(b[3][2] == b32);
    3211           1 :             CATCH_REQUIRE(b[3][3] == b33);
    3212             : 
    3213             :             // run operation A *= B
    3214             :             //
    3215           1 :             a *= b;
    3216             : 
    3217             :             // this could fail
    3218             :             //
    3219             :             //CATCH_REQUIRE_FALSE(a[0][0] == a00);
    3220             :             //CATCH_REQUIRE_FALSE(a[0][1] == a01);
    3221             :             //CATCH_REQUIRE_FALSE(a[0][2] == a02);
    3222             :             //CATCH_REQUIRE_FALSE(a[0][3] == a03);
    3223             :             //CATCH_REQUIRE_FALSE(a[1][0] == a10);
    3224             :             //CATCH_REQUIRE_FALSE(a[1][1] == a11);
    3225             :             //CATCH_REQUIRE_FALSE(a[1][2] == a12);
    3226             :             //CATCH_REQUIRE_FALSE(a[1][3] == a13);
    3227             :             //CATCH_REQUIRE_FALSE(a[2][0] == a20);
    3228             :             //CATCH_REQUIRE_FALSE(a[2][1] == a21);
    3229             :             //CATCH_REQUIRE_FALSE(a[2][2] == a22);
    3230             :             //CATCH_REQUIRE_FALSE(a[2][3] == a23);
    3231             :             //CATCH_REQUIRE_FALSE(a[3][0] == a30);
    3232             :             //CATCH_REQUIRE_FALSE(a[3][1] == a31);
    3233             :             //CATCH_REQUIRE_FALSE(a[3][2] == a32);
    3234             :             //CATCH_REQUIRE_FALSE(a[3][3] == a33);
    3235             : 
    3236           1 :             CATCH_REQUIRE(b[0][0] == b00);
    3237           1 :             CATCH_REQUIRE(b[0][1] == b01);
    3238           1 :             CATCH_REQUIRE(b[0][2] == b02);
    3239           1 :             CATCH_REQUIRE(b[0][3] == b03);
    3240           1 :             CATCH_REQUIRE(b[1][0] == b10);
    3241           1 :             CATCH_REQUIRE(b[1][1] == b11);
    3242           1 :             CATCH_REQUIRE(b[1][2] == b12);
    3243           1 :             CATCH_REQUIRE(b[1][3] == b13);
    3244           1 :             CATCH_REQUIRE(b[2][0] == b20);
    3245           1 :             CATCH_REQUIRE(b[2][1] == b21);
    3246           1 :             CATCH_REQUIRE(b[2][2] == b22);
    3247           1 :             CATCH_REQUIRE(b[2][3] == b23);
    3248           1 :             CATCH_REQUIRE(b[3][0] == b30);
    3249           1 :             CATCH_REQUIRE(b[3][1] == b31);
    3250           1 :             CATCH_REQUIRE(b[3][2] == b32);
    3251           1 :             CATCH_REQUIRE(b[3][3] == b33);
    3252             : 
    3253           1 :             CATCH_REQUIRE(a[0][0] == a00 * b00 + a01 * b10 + a02 * b20 + a03 * b30);
    3254           1 :             CATCH_REQUIRE(a[0][1] == a00 * b01 + a01 * b11 + a02 * b21 + a03 * b31);
    3255           1 :             CATCH_REQUIRE(a[0][2] == a00 * b02 + a01 * b12 + a02 * b22 + a03 * b32);
    3256           1 :             CATCH_REQUIRE(a[0][3] == a00 * b03 + a01 * b13 + a02 * b23 + a03 * b33);
    3257             : 
    3258           1 :             CATCH_REQUIRE(a[1][0] == a10 * b00 + a11 * b10 + a12 * b20 + a13 * b30);
    3259           1 :             CATCH_REQUIRE(a[1][1] == a10 * b01 + a11 * b11 + a12 * b21 + a13 * b31);
    3260           1 :             CATCH_REQUIRE(a[1][2] == a10 * b02 + a11 * b12 + a12 * b22 + a13 * b32);
    3261           1 :             CATCH_REQUIRE(a[1][3] == a10 * b03 + a11 * b13 + a12 * b23 + a13 * b33);
    3262             : 
    3263           1 :             CATCH_REQUIRE(a[2][0] == a20 * b00 + a21 * b10 + a22 * b20 + a23 * b30);
    3264           1 :             CATCH_REQUIRE(a[2][1] == a20 * b01 + a21 * b11 + a22 * b21 + a23 * b31);
    3265           1 :             CATCH_REQUIRE(a[2][2] == a20 * b02 + a21 * b12 + a22 * b22 + a23 * b32);
    3266           1 :             CATCH_REQUIRE(a[2][3] == a20 * b03 + a21 * b13 + a22 * b23 + a23 * b33);
    3267             : 
    3268           1 :             CATCH_REQUIRE(a[3][0] == a30 * b00 + a31 * b10 + a32 * b20 + a33 * b30);
    3269           1 :             CATCH_REQUIRE(a[3][1] == a30 * b01 + a31 * b11 + a32 * b21 + a33 * b31);
    3270           1 :             CATCH_REQUIRE(a[3][2] == a30 * b02 + a31 * b12 + a32 * b22 + a33 * b32);
    3271           1 :             CATCH_REQUIRE(a[3][3] == a30 * b03 + a31 * b13 + a32 * b23 + a33 * b33);
    3272             : 
    3273             :             //CATCH_REQUIRE(a[0][0] == a00 * b00 + a10 * b01 + a20 * b02 + a30 * b03);
    3274             :             //CATCH_REQUIRE(a[0][1] == a01 * b00 + a11 * b01 + a21 * b02 + a31 * b03);
    3275             :             //CATCH_REQUIRE(a[0][2] == a02 * b00 + a12 * b01 + a22 * b02 + a32 * b03);
    3276             :             //CATCH_REQUIRE(a[0][3] == a03 * b00 + a13 * b01 + a23 * b02 + a33 * b03);
    3277             : 
    3278             :             //CATCH_REQUIRE(a[1][0] == a00 * b10 + a10 * b11 + a20 * b12 + a30 * b13);
    3279             :             //CATCH_REQUIRE(a[1][1] == a01 * b10 + a11 * b11 + a21 * b12 + a31 * b13);
    3280             :             //CATCH_REQUIRE(a[1][2] == a02 * b10 + a12 * b11 + a22 * b12 + a32 * b13);
    3281             :             //CATCH_REQUIRE(a[1][3] == a03 * b10 + a13 * b11 + a23 * b12 + a33 * b13);
    3282             : 
    3283             :             //CATCH_REQUIRE(a[2][0] == a00 * b20 + a10 * b21 + a20 * b22 + a30 * b23);
    3284             :             //CATCH_REQUIRE(a[2][1] == a01 * b20 + a11 * b21 + a21 * b22 + a31 * b23);
    3285             :             //CATCH_REQUIRE(a[2][2] == a02 * b20 + a12 * b21 + a22 * b22 + a32 * b23);
    3286             :             //CATCH_REQUIRE(a[2][3] == a03 * b20 + a13 * b21 + a23 * b22 + a33 * b23);
    3287             : 
    3288             :             //CATCH_REQUIRE(a[3][0] == a00 * b30 + a10 * b31 + a20 * b32 + a30 * b33);
    3289             :             //CATCH_REQUIRE(a[3][1] == a01 * b30 + a11 * b31 + a21 * b32 + a31 * b33);
    3290             :             //CATCH_REQUIRE(a[3][2] == a02 * b30 + a12 * b31 + a22 * b32 + a32 * b33);
    3291             :             //CATCH_REQUIRE(a[3][3] == a03 * b30 + a13 * b31 + a23 * b32 + a33 * b33);
    3292             :         }
    3293             :         CATCH_END_SECTION()
    3294             :     }
    3295             : 
    3296             :     // create two random 4x4 matrices and make sure the add works
    3297             :     //
    3298          16 :     CATCH_GIVEN("div")
    3299             :     {
    3300           8 :         CATCH_START_SECTION("matrix: b=a/<scalar>")
    3301             :         {
    3302             :             // setup A
    3303             :             //
    3304           2 :             snapdev::matrix<double> a(4, 4);
    3305             : 
    3306           1 :             CATCH_REQUIRE(a.rows() == 4);
    3307           1 :             CATCH_REQUIRE(a.columns() == 4);
    3308             : 
    3309           1 :             double const a00(frand());
    3310           1 :             double const a01(frand());
    3311           1 :             double const a02(frand());
    3312           1 :             double const a03(frand());
    3313           1 :             double const a10(frand());
    3314           1 :             double const a11(frand());
    3315           1 :             double const a12(frand());
    3316           1 :             double const a13(frand());
    3317           1 :             double const a20(frand());
    3318           1 :             double const a21(frand());
    3319           1 :             double const a22(frand());
    3320           1 :             double const a23(frand());
    3321           1 :             double const a30(frand());
    3322           1 :             double const a31(frand());
    3323           1 :             double const a32(frand());
    3324           1 :             double const a33(frand());
    3325             : 
    3326           1 :             a[0][0] = a00;
    3327           1 :             a[0][1] = a01;
    3328           1 :             a[0][2] = a02;
    3329           1 :             a[0][3] = a03;
    3330           1 :             a[1][0] = a10;
    3331           1 :             a[1][1] = a11;
    3332           1 :             a[1][2] = a12;
    3333           1 :             a[1][3] = a13;
    3334           1 :             a[2][0] = a20;
    3335           1 :             a[2][1] = a21;
    3336           1 :             a[2][2] = a22;
    3337           1 :             a[2][3] = a23;
    3338           1 :             a[3][0] = a30;
    3339           1 :             a[3][1] = a31;
    3340           1 :             a[3][2] = a32;
    3341           1 :             a[3][3] = a33;
    3342             : 
    3343           1 :             CATCH_REQUIRE(a[0][0] == a00);
    3344           1 :             CATCH_REQUIRE(a[0][1] == a01);
    3345           1 :             CATCH_REQUIRE(a[0][2] == a02);
    3346           1 :             CATCH_REQUIRE(a[0][3] == a03);
    3347           1 :             CATCH_REQUIRE(a[1][0] == a10);
    3348           1 :             CATCH_REQUIRE(a[1][1] == a11);
    3349           1 :             CATCH_REQUIRE(a[1][2] == a12);
    3350           1 :             CATCH_REQUIRE(a[1][3] == a13);
    3351           1 :             CATCH_REQUIRE(a[2][0] == a20);
    3352           1 :             CATCH_REQUIRE(a[2][1] == a21);
    3353           1 :             CATCH_REQUIRE(a[2][2] == a22);
    3354           1 :             CATCH_REQUIRE(a[2][3] == a23);
    3355           1 :             CATCH_REQUIRE(a[3][0] == a30);
    3356           1 :             CATCH_REQUIRE(a[3][1] == a31);
    3357           1 :             CATCH_REQUIRE(a[3][2] == a32);
    3358           1 :             CATCH_REQUIRE(a[3][3] == a33);
    3359             : 
    3360             :             // setup B
    3361             :             //
    3362           2 :             snapdev::matrix<double> b(4, 4);
    3363             : 
    3364           1 :             CATCH_REQUIRE(b.rows() == 4);
    3365           1 :             CATCH_REQUIRE(b.columns() == 4);
    3366             : 
    3367           1 :             double const b00(frand());
    3368           1 :             double const b01(frand());
    3369           1 :             double const b02(frand());
    3370           1 :             double const b03(frand());
    3371           1 :             double const b10(frand());
    3372           1 :             double const b11(frand());
    3373           1 :             double const b12(frand());
    3374           1 :             double const b13(frand());
    3375           1 :             double const b20(frand());
    3376           1 :             double const b21(frand());
    3377           1 :             double const b22(frand());
    3378           1 :             double const b23(frand());
    3379           1 :             double const b30(frand());
    3380           1 :             double const b31(frand());
    3381           1 :             double const b32(frand());
    3382           1 :             double const b33(frand());
    3383             : 
    3384           1 :             b[0][0] = b00;
    3385           1 :             b[0][1] = b01;
    3386           1 :             b[0][2] = b02;
    3387           1 :             b[0][3] = b03;
    3388           1 :             b[1][0] = b10;
    3389           1 :             b[1][1] = b11;
    3390           1 :             b[1][2] = b12;
    3391           1 :             b[1][3] = b13;
    3392           1 :             b[2][0] = b20;
    3393           1 :             b[2][1] = b21;
    3394           1 :             b[2][2] = b22;
    3395           1 :             b[2][3] = b23;
    3396           1 :             b[3][0] = b30;
    3397           1 :             b[3][1] = b31;
    3398           1 :             b[3][2] = b32;
    3399           1 :             b[3][3] = b33;
    3400             : 
    3401           1 :             CATCH_REQUIRE(b[0][0] == b00);
    3402           1 :             CATCH_REQUIRE(b[0][1] == b01);
    3403           1 :             CATCH_REQUIRE(b[0][2] == b02);
    3404           1 :             CATCH_REQUIRE(b[0][3] == b03);
    3405           1 :             CATCH_REQUIRE(b[1][0] == b10);
    3406           1 :             CATCH_REQUIRE(b[1][1] == b11);
    3407           1 :             CATCH_REQUIRE(b[1][2] == b12);
    3408           1 :             CATCH_REQUIRE(b[1][3] == b13);
    3409           1 :             CATCH_REQUIRE(b[2][0] == b20);
    3410           1 :             CATCH_REQUIRE(b[2][1] == b21);
    3411           1 :             CATCH_REQUIRE(b[2][2] == b22);
    3412           1 :             CATCH_REQUIRE(b[2][3] == b23);
    3413           1 :             CATCH_REQUIRE(b[3][0] == b30);
    3414           1 :             CATCH_REQUIRE(b[3][1] == b31);
    3415           1 :             CATCH_REQUIRE(b[3][2] == b32);
    3416           1 :             CATCH_REQUIRE(b[3][3] == b33);
    3417             : 
    3418             :             // create a scalar for our test
    3419             :             //
    3420           1 :             double scalar(frand());
    3421           1 :             while(scalar == 0.0 || scalar == 1.0)
    3422             :             {
    3423           0 :                 scalar = frand();
    3424             :             }
    3425             : 
    3426             :             // run operation B = A / <scalar>
    3427             :             //
    3428           1 :             b = a / scalar;
    3429             : 
    3430           1 :             CATCH_REQUIRE(a[0][0] == a00);
    3431           1 :             CATCH_REQUIRE(a[0][1] == a01);
    3432           1 :             CATCH_REQUIRE(a[0][2] == a02);
    3433           1 :             CATCH_REQUIRE(a[0][3] == a03);
    3434           1 :             CATCH_REQUIRE(a[1][0] == a10);
    3435           1 :             CATCH_REQUIRE(a[1][1] == a11);
    3436           1 :             CATCH_REQUIRE(a[1][2] == a12);
    3437           1 :             CATCH_REQUIRE(a[1][3] == a13);
    3438           1 :             CATCH_REQUIRE(a[2][0] == a20);
    3439           1 :             CATCH_REQUIRE(a[2][1] == a21);
    3440           1 :             CATCH_REQUIRE(a[2][2] == a22);
    3441           1 :             CATCH_REQUIRE(a[2][3] == a23);
    3442           1 :             CATCH_REQUIRE(a[3][0] == a30);
    3443           1 :             CATCH_REQUIRE(a[3][1] == a31);
    3444           1 :             CATCH_REQUIRE(a[3][2] == a32);
    3445           1 :             CATCH_REQUIRE(a[3][3] == a33);
    3446             : 
    3447             :             // this can't fail because we ensure scalar != 0 or 1
    3448             :             //
    3449           1 :             CATCH_REQUIRE_FALSE(b[0][0] == b00);
    3450           1 :             CATCH_REQUIRE_FALSE(b[0][1] == b01);
    3451           1 :             CATCH_REQUIRE_FALSE(b[0][2] == b02);
    3452           1 :             CATCH_REQUIRE_FALSE(b[0][3] == b03);
    3453           1 :             CATCH_REQUIRE_FALSE(b[1][0] == b10);
    3454           1 :             CATCH_REQUIRE_FALSE(b[1][1] == b11);
    3455           1 :             CATCH_REQUIRE_FALSE(b[1][2] == b12);
    3456           1 :             CATCH_REQUIRE_FALSE(b[1][3] == b13);
    3457           1 :             CATCH_REQUIRE_FALSE(b[2][0] == b20);
    3458           1 :             CATCH_REQUIRE_FALSE(b[2][1] == b21);
    3459           1 :             CATCH_REQUIRE_FALSE(b[2][2] == b22);
    3460           1 :             CATCH_REQUIRE_FALSE(b[2][3] == b23);
    3461           1 :             CATCH_REQUIRE_FALSE(b[3][0] == b30);
    3462           1 :             CATCH_REQUIRE_FALSE(b[3][1] == b31);
    3463           1 :             CATCH_REQUIRE_FALSE(b[3][2] == b32);
    3464           1 :             CATCH_REQUIRE_FALSE(b[3][3] == b33);
    3465             : 
    3466           1 :             CATCH_REQUIRE(b[0][0] == a00 / scalar);
    3467           1 :             CATCH_REQUIRE(b[0][1] == a01 / scalar);
    3468           1 :             CATCH_REQUIRE(b[0][2] == a02 / scalar);
    3469           1 :             CATCH_REQUIRE(b[0][3] == a03 / scalar);
    3470           1 :             CATCH_REQUIRE(b[1][0] == a10 / scalar);
    3471           1 :             CATCH_REQUIRE(b[1][1] == a11 / scalar);
    3472           1 :             CATCH_REQUIRE(b[1][2] == a12 / scalar);
    3473           1 :             CATCH_REQUIRE(b[1][3] == a13 / scalar);
    3474           1 :             CATCH_REQUIRE(b[2][0] == a20 / scalar);
    3475           1 :             CATCH_REQUIRE(b[2][1] == a21 / scalar);
    3476           1 :             CATCH_REQUIRE(b[2][2] == a22 / scalar);
    3477           1 :             CATCH_REQUIRE(b[2][3] == a23 / scalar);
    3478           1 :             CATCH_REQUIRE(b[3][0] == a30 / scalar);
    3479           1 :             CATCH_REQUIRE(b[3][1] == a31 / scalar);
    3480           1 :             CATCH_REQUIRE(b[3][2] == a32 / scalar);
    3481           1 :             CATCH_REQUIRE(b[3][3] == a33 / scalar);
    3482             :         }
    3483             :         CATCH_END_SECTION()
    3484             : 
    3485           8 :         CATCH_START_SECTION("matrix: a/=<scalar>")
    3486             :         {
    3487             :             // setup A
    3488             :             //
    3489           2 :             snapdev::matrix<double> a(4, 4);
    3490             : 
    3491           1 :             CATCH_REQUIRE(a.rows() == 4);
    3492           1 :             CATCH_REQUIRE(a.columns() == 4);
    3493             : 
    3494           1 :             double const a00(frand());
    3495           1 :             double const a01(frand());
    3496           1 :             double const a02(frand());
    3497           1 :             double const a03(frand());
    3498           1 :             double const a10(frand());
    3499           1 :             double const a11(frand());
    3500           1 :             double const a12(frand());
    3501           1 :             double const a13(frand());
    3502           1 :             double const a20(frand());
    3503           1 :             double const a21(frand());
    3504           1 :             double const a22(frand());
    3505           1 :             double const a23(frand());
    3506           1 :             double const a30(frand());
    3507           1 :             double const a31(frand());
    3508           1 :             double const a32(frand());
    3509           1 :             double const a33(frand());
    3510             : 
    3511           1 :             a[0][0] = a00;
    3512           1 :             a[0][1] = a01;
    3513           1 :             a[0][2] = a02;
    3514           1 :             a[0][3] = a03;
    3515           1 :             a[1][0] = a10;
    3516           1 :             a[1][1] = a11;
    3517           1 :             a[1][2] = a12;
    3518           1 :             a[1][3] = a13;
    3519           1 :             a[2][0] = a20;
    3520           1 :             a[2][1] = a21;
    3521           1 :             a[2][2] = a22;
    3522           1 :             a[2][3] = a23;
    3523           1 :             a[3][0] = a30;
    3524           1 :             a[3][1] = a31;
    3525           1 :             a[3][2] = a32;
    3526           1 :             a[3][3] = a33;
    3527             : 
    3528           1 :             CATCH_REQUIRE(a[0][0] == a00);
    3529           1 :             CATCH_REQUIRE(a[0][1] == a01);
    3530           1 :             CATCH_REQUIRE(a[0][2] == a02);
    3531           1 :             CATCH_REQUIRE(a[0][3] == a03);
    3532           1 :             CATCH_REQUIRE(a[1][0] == a10);
    3533           1 :             CATCH_REQUIRE(a[1][1] == a11);
    3534           1 :             CATCH_REQUIRE(a[1][2] == a12);
    3535           1 :             CATCH_REQUIRE(a[1][3] == a13);
    3536           1 :             CATCH_REQUIRE(a[2][0] == a20);
    3537           1 :             CATCH_REQUIRE(a[2][1] == a21);
    3538           1 :             CATCH_REQUIRE(a[2][2] == a22);
    3539           1 :             CATCH_REQUIRE(a[2][3] == a23);
    3540           1 :             CATCH_REQUIRE(a[3][0] == a30);
    3541           1 :             CATCH_REQUIRE(a[3][1] == a31);
    3542           1 :             CATCH_REQUIRE(a[3][2] == a32);
    3543           1 :             CATCH_REQUIRE(a[3][3] == a33);
    3544             : 
    3545             :             // create a scalar for our test
    3546             :             //
    3547           1 :             double scalar(frand());
    3548           1 :             while(scalar == 0.0)
    3549             :             {
    3550           0 :                 scalar = frand();
    3551             :             }
    3552             : 
    3553             :             // run operation A /= <scalar>
    3554             :             //
    3555           1 :             a /= scalar;
    3556             : 
    3557             :             // this can't fail because we ensure scalar != 0 or 1
    3558             :             //
    3559           1 :             CATCH_REQUIRE_FALSE(a[0][0] == a00);
    3560           1 :             CATCH_REQUIRE_FALSE(a[0][1] == a01);
    3561           1 :             CATCH_REQUIRE_FALSE(a[0][2] == a02);
    3562           1 :             CATCH_REQUIRE_FALSE(a[0][3] == a03);
    3563           1 :             CATCH_REQUIRE_FALSE(a[1][0] == a10);
    3564           1 :             CATCH_REQUIRE_FALSE(a[1][1] == a11);
    3565           1 :             CATCH_REQUIRE_FALSE(a[1][2] == a12);
    3566           1 :             CATCH_REQUIRE_FALSE(a[1][3] == a13);
    3567           1 :             CATCH_REQUIRE_FALSE(a[2][0] == a20);
    3568           1 :             CATCH_REQUIRE_FALSE(a[2][1] == a21);
    3569           1 :             CATCH_REQUIRE_FALSE(a[2][2] == a22);
    3570           1 :             CATCH_REQUIRE_FALSE(a[2][3] == a23);
    3571           1 :             CATCH_REQUIRE_FALSE(a[3][0] == a30);
    3572           1 :             CATCH_REQUIRE_FALSE(a[3][1] == a31);
    3573           1 :             CATCH_REQUIRE_FALSE(a[3][2] == a32);
    3574           1 :             CATCH_REQUIRE_FALSE(a[3][3] == a33);
    3575             : 
    3576           1 :             CATCH_REQUIRE(a[0][0] == a00 / scalar);
    3577           1 :             CATCH_REQUIRE(a[0][1] == a01 / scalar);
    3578           1 :             CATCH_REQUIRE(a[0][2] == a02 / scalar);
    3579           1 :             CATCH_REQUIRE(a[0][3] == a03 / scalar);
    3580           1 :             CATCH_REQUIRE(a[1][0] == a10 / scalar);
    3581           1 :             CATCH_REQUIRE(a[1][1] == a11 / scalar);
    3582           1 :             CATCH_REQUIRE(a[1][2] == a12 / scalar);
    3583           1 :             CATCH_REQUIRE(a[1][3] == a13 / scalar);
    3584           1 :             CATCH_REQUIRE(a[2][0] == a20 / scalar);
    3585           1 :             CATCH_REQUIRE(a[2][1] == a21 / scalar);
    3586           1 :             CATCH_REQUIRE(a[2][2] == a22 / scalar);
    3587           1 :             CATCH_REQUIRE(a[2][3] == a23 / scalar);
    3588           1 :             CATCH_REQUIRE(a[3][0] == a30 / scalar);
    3589           1 :             CATCH_REQUIRE(a[3][1] == a31 / scalar);
    3590           1 :             CATCH_REQUIRE(a[3][2] == a32 / scalar);
    3591           1 :             CATCH_REQUIRE(a[3][3] == a33 / scalar);
    3592             :         }
    3593             :         CATCH_END_SECTION()
    3594             : 
    3595           8 :         CATCH_START_SECTION("matrix: c=a/b")
    3596             :         {
    3597             :             // setup A
    3598             :             //
    3599           2 :             snapdev::matrix<double> a(4, 4);
    3600             : 
    3601           1 :             CATCH_REQUIRE(a.rows() == 4);
    3602           1 :             CATCH_REQUIRE(a.columns() == 4);
    3603             : 
    3604           1 :             double const a00(frand());
    3605           1 :             double const a01(frand());
    3606           1 :             double const a02(frand());
    3607           1 :             double const a03(frand());
    3608           1 :             double const a10(frand());
    3609           1 :             double const a11(frand());
    3610           1 :             double const a12(frand());
    3611           1 :             double const a13(frand());
    3612           1 :             double const a20(frand());
    3613           1 :             double const a21(frand());
    3614           1 :             double const a22(frand());
    3615           1 :             double const a23(frand());
    3616           1 :             double const a30(frand());
    3617           1 :             double const a31(frand());
    3618           1 :             double const a32(frand());
    3619           1 :             double const a33(frand());
    3620             : 
    3621           1 :             a[0][0] = a00;
    3622           1 :             a[0][1] = a01;
    3623           1 :             a[0][2] = a02;
    3624           1 :             a[0][3] = a03;
    3625           1 :             a[1][0] = a10;
    3626           1 :             a[1][1] = a11;
    3627           1 :             a[1][2] = a12;
    3628           1 :             a[1][3] = a13;
    3629           1 :             a[2][0] = a20;
    3630           1 :             a[2][1] = a21;
    3631           1 :             a[2][2] = a22;
    3632           1 :             a[2][3] = a23;
    3633           1 :             a[3][0] = a30;
    3634           1 :             a[3][1] = a31;
    3635           1 :             a[3][2] = a32;
    3636           1 :             a[3][3] = a33;
    3637             : 
    3638           1 :             CATCH_REQUIRE(a[0][0] == a00);
    3639           1 :             CATCH_REQUIRE(a[0][1] == a01);
    3640           1 :             CATCH_REQUIRE(a[0][2] == a02);
    3641           1 :             CATCH_REQUIRE(a[0][3] == a03);
    3642           1 :             CATCH_REQUIRE(a[1][0] == a10);
    3643           1 :             CATCH_REQUIRE(a[1][1] == a11);
    3644           1 :             CATCH_REQUIRE(a[1][2] == a12);
    3645           1 :             CATCH_REQUIRE(a[1][3] == a13);
    3646           1 :             CATCH_REQUIRE(a[2][0] == a20);
    3647           1 :             CATCH_REQUIRE(a[2][1] == a21);
    3648           1 :             CATCH_REQUIRE(a[2][2] == a22);
    3649           1 :             CATCH_REQUIRE(a[2][3] == a23);
    3650           1 :             CATCH_REQUIRE(a[3][0] == a30);
    3651           1 :             CATCH_REQUIRE(a[3][1] == a31);
    3652           1 :             CATCH_REQUIRE(a[3][2] == a32);
    3653           1 :             CATCH_REQUIRE(a[3][3] == a33);
    3654             : 
    3655             :             // setup B
    3656             :             //
    3657           2 :             snapdev::matrix<double> b(4, 4);
    3658             : 
    3659           1 :             CATCH_REQUIRE(b.rows() == 4);
    3660           1 :             CATCH_REQUIRE(b.columns() == 4);
    3661             : 
    3662             :             double b00;
    3663             :             double b01;
    3664             :             double b02;
    3665             :             double b03;
    3666             :             double b10;
    3667             :             double b11;
    3668             :             double b12;
    3669             :             double b13;
    3670             :             double b20;
    3671             :             double b21;
    3672             :             double b22;
    3673             :             double b23;
    3674             :             double b30;
    3675             :             double b31;
    3676             :             double b32;
    3677             :             double b33;
    3678             : 
    3679             :             // create an inversible matrix (most are with random numbers)
    3680           2 :             snapdev::matrix<double> t(4, 4);
    3681           0 :             do
    3682             :             {
    3683           1 :                 b00 = frand();
    3684           1 :                 b01 = frand();
    3685           1 :                 b02 = frand();
    3686           1 :                 b03 = frand();
    3687           1 :                 b10 = frand();
    3688           1 :                 b11 = frand();
    3689           1 :                 b12 = frand();
    3690           1 :                 b13 = frand();
    3691           1 :                 b20 = frand();
    3692           1 :                 b21 = frand();
    3693           1 :                 b22 = frand();
    3694           1 :                 b23 = frand();
    3695           1 :                 b30 = frand();
    3696           1 :                 b31 = frand();
    3697           1 :                 b32 = frand();
    3698           1 :                 b33 = frand();
    3699             : 
    3700           1 :                 b[0][0] = b00;
    3701           1 :                 b[0][1] = b01;
    3702           1 :                 b[0][2] = b02;
    3703           1 :                 b[0][3] = b03;
    3704           1 :                 b[1][0] = b10;
    3705           1 :                 b[1][1] = b11;
    3706           1 :                 b[1][2] = b12;
    3707           1 :                 b[1][3] = b13;
    3708           1 :                 b[2][0] = b20;
    3709           1 :                 b[2][1] = b21;
    3710           1 :                 b[2][2] = b22;
    3711           1 :                 b[2][3] = b23;
    3712           1 :                 b[3][0] = b30;
    3713           1 :                 b[3][1] = b31;
    3714           1 :                 b[3][2] = b32;
    3715           1 :                 b[3][3] = b33;
    3716             : 
    3717           1 :                 CATCH_REQUIRE(b[0][0] == b00);
    3718           1 :                 CATCH_REQUIRE(b[0][1] == b01);
    3719           1 :                 CATCH_REQUIRE(b[0][2] == b02);
    3720           1 :                 CATCH_REQUIRE(b[0][3] == b03);
    3721           1 :                 CATCH_REQUIRE(b[1][0] == b10);
    3722           1 :                 CATCH_REQUIRE(b[1][1] == b11);
    3723           1 :                 CATCH_REQUIRE(b[1][2] == b12);
    3724           1 :                 CATCH_REQUIRE(b[1][3] == b13);
    3725           1 :                 CATCH_REQUIRE(b[2][0] == b20);
    3726           1 :                 CATCH_REQUIRE(b[2][1] == b21);
    3727           1 :                 CATCH_REQUIRE(b[2][2] == b22);
    3728           1 :                 CATCH_REQUIRE(b[2][3] == b23);
    3729           1 :                 CATCH_REQUIRE(b[3][0] == b30);
    3730           1 :                 CATCH_REQUIRE(b[3][1] == b31);
    3731           1 :                 CATCH_REQUIRE(b[3][2] == b32);
    3732           1 :                 CATCH_REQUIRE(b[3][3] == b33);
    3733             : 
    3734           1 :                 t = b;
    3735             :             }
    3736           1 :             while(!t.inverse());
    3737             : 
    3738             :             // setup C
    3739             :             //
    3740           2 :             snapdev::matrix<double> c(4, 4);
    3741             : 
    3742           1 :             CATCH_REQUIRE(c.rows() == 4);
    3743           1 :             CATCH_REQUIRE(c.columns() == 4);
    3744             : 
    3745           1 :             double const c00(frand());
    3746           1 :             double const c01(frand());
    3747           1 :             double const c02(frand());
    3748           1 :             double const c03(frand());
    3749           1 :             double const c10(frand());
    3750           1 :             double const c11(frand());
    3751           1 :             double const c12(frand());
    3752           1 :             double const c13(frand());
    3753           1 :             double const c20(frand());
    3754           1 :             double const c21(frand());
    3755           1 :             double const c22(frand());
    3756           1 :             double const c23(frand());
    3757           1 :             double const c30(frand());
    3758           1 :             double const c31(frand());
    3759           1 :             double const c32(frand());
    3760           1 :             double const c33(frand());
    3761             : 
    3762           1 :             c[0][0] = c00;
    3763           1 :             c[0][1] = c01;
    3764           1 :             c[0][2] = c02;
    3765           1 :             c[0][3] = c03;
    3766           1 :             c[1][0] = c10;
    3767           1 :             c[1][1] = c11;
    3768           1 :             c[1][2] = c12;
    3769           1 :             c[1][3] = c13;
    3770           1 :             c[2][0] = c20;
    3771           1 :             c[2][1] = c21;
    3772           1 :             c[2][2] = c22;
    3773           1 :             c[2][3] = c23;
    3774           1 :             c[3][0] = c30;
    3775           1 :             c[3][1] = c31;
    3776           1 :             c[3][2] = c32;
    3777           1 :             c[3][3] = c33;
    3778             : 
    3779           1 :             CATCH_REQUIRE(c[0][0] == c00);
    3780           1 :             CATCH_REQUIRE(c[0][1] == c01);
    3781           1 :             CATCH_REQUIRE(c[0][2] == c02);
    3782           1 :             CATCH_REQUIRE(c[0][3] == c03);
    3783           1 :             CATCH_REQUIRE(c[1][0] == c10);
    3784           1 :             CATCH_REQUIRE(c[1][1] == c11);
    3785           1 :             CATCH_REQUIRE(c[1][2] == c12);
    3786           1 :             CATCH_REQUIRE(c[1][3] == c13);
    3787           1 :             CATCH_REQUIRE(c[2][0] == c20);
    3788           1 :             CATCH_REQUIRE(c[2][1] == c21);
    3789           1 :             CATCH_REQUIRE(c[2][2] == c22);
    3790           1 :             CATCH_REQUIRE(c[2][3] == c23);
    3791           1 :             CATCH_REQUIRE(c[3][0] == c30);
    3792           1 :             CATCH_REQUIRE(c[3][1] == c31);
    3793           1 :             CATCH_REQUIRE(c[3][2] == c32);
    3794           1 :             CATCH_REQUIRE(c[3][3] == c33);
    3795             : 
    3796             :             // run operation C = A / B
    3797             :             //
    3798           1 :             c = a / b;
    3799             : 
    3800           1 :             CATCH_REQUIRE(a[0][0] == a00);
    3801           1 :             CATCH_REQUIRE(a[0][1] == a01);
    3802           1 :             CATCH_REQUIRE(a[0][2] == a02);
    3803           1 :             CATCH_REQUIRE(a[0][3] == a03);
    3804           1 :             CATCH_REQUIRE(a[1][0] == a10);
    3805           1 :             CATCH_REQUIRE(a[1][1] == a11);
    3806           1 :             CATCH_REQUIRE(a[1][2] == a12);
    3807           1 :             CATCH_REQUIRE(a[1][3] == a13);
    3808           1 :             CATCH_REQUIRE(a[2][0] == a20);
    3809           1 :             CATCH_REQUIRE(a[2][1] == a21);
    3810           1 :             CATCH_REQUIRE(a[2][2] == a22);
    3811           1 :             CATCH_REQUIRE(a[2][3] == a23);
    3812           1 :             CATCH_REQUIRE(a[3][0] == a30);
    3813           1 :             CATCH_REQUIRE(a[3][1] == a31);
    3814           1 :             CATCH_REQUIRE(a[3][2] == a32);
    3815           1 :             CATCH_REQUIRE(a[3][3] == a33);
    3816             : 
    3817           1 :             CATCH_REQUIRE(b[0][0] == b00);
    3818           1 :             CATCH_REQUIRE(b[0][1] == b01);
    3819           1 :             CATCH_REQUIRE(b[0][2] == b02);
    3820           1 :             CATCH_REQUIRE(b[0][3] == b03);
    3821           1 :             CATCH_REQUIRE(b[1][0] == b10);
    3822           1 :             CATCH_REQUIRE(b[1][1] == b11);
    3823           1 :             CATCH_REQUIRE(b[1][2] == b12);
    3824           1 :             CATCH_REQUIRE(b[1][3] == b13);
    3825           1 :             CATCH_REQUIRE(b[2][0] == b20);
    3826           1 :             CATCH_REQUIRE(b[2][1] == b21);
    3827           1 :             CATCH_REQUIRE(b[2][2] == b22);
    3828           1 :             CATCH_REQUIRE(b[2][3] == b23);
    3829           1 :             CATCH_REQUIRE(b[3][0] == b30);
    3830           1 :             CATCH_REQUIRE(b[3][1] == b31);
    3831           1 :             CATCH_REQUIRE(b[3][2] == b32);
    3832           1 :             CATCH_REQUIRE(b[3][3] == b33);
    3833             : 
    3834             :             // this could fail
    3835             :             //CATCH_REQUIRE_FALSE(c[0][0] == c00);
    3836             :             //CATCH_REQUIRE_FALSE(c[0][1] == c01);
    3837             :             //CATCH_REQUIRE_FALSE(c[0][2] == c02);
    3838             :             //CATCH_REQUIRE_FALSE(c[0][3] == c03);
    3839             :             //CATCH_REQUIRE_FALSE(c[1][0] == c10);
    3840             :             //CATCH_REQUIRE_FALSE(c[1][1] == c11);
    3841             :             //CATCH_REQUIRE_FALSE(c[1][2] == c12);
    3842             :             //CATCH_REQUIRE_FALSE(c[1][3] == c13);
    3843             :             //CATCH_REQUIRE_FALSE(c[2][0] == c20);
    3844             :             //CATCH_REQUIRE_FALSE(c[2][1] == c21);
    3845             :             //CATCH_REQUIRE_FALSE(c[2][2] == c22);
    3846             :             //CATCH_REQUIRE_FALSE(c[2][3] == c23);
    3847             :             //CATCH_REQUIRE_FALSE(c[3][0] == c30);
    3848             :             //CATCH_REQUIRE_FALSE(c[3][1] == c31);
    3849             :             //CATCH_REQUIRE_FALSE(c[3][2] == c32);
    3850             :             //CATCH_REQUIRE_FALSE(c[3][3] == c33);
    3851             : 
    3852             :             // first level verification, which is exactly what the
    3853             :             // division function does, it does not mean that the
    3854             :             // division works properly, though
    3855             :             //
    3856           2 :             snapdev::matrix<double> r(4, 4);
    3857           1 :             r = a * t;
    3858             : 
    3859           1 :             CATCH_REQUIRE(c[0][0] == r[0][0]);
    3860           1 :             CATCH_REQUIRE(c[0][1] == r[0][1]);
    3861           1 :             CATCH_REQUIRE(c[0][2] == r[0][2]);
    3862           1 :             CATCH_REQUIRE(c[0][3] == r[0][3]);
    3863           1 :             CATCH_REQUIRE(c[1][0] == r[1][0]);
    3864           1 :             CATCH_REQUIRE(c[1][1] == r[1][1]);
    3865           1 :             CATCH_REQUIRE(c[1][2] == r[1][2]);
    3866           1 :             CATCH_REQUIRE(c[1][3] == r[1][3]);
    3867           1 :             CATCH_REQUIRE(c[2][0] == r[2][0]);
    3868           1 :             CATCH_REQUIRE(c[2][1] == r[2][1]);
    3869           1 :             CATCH_REQUIRE(c[2][2] == r[2][2]);
    3870           1 :             CATCH_REQUIRE(c[2][3] == r[2][3]);
    3871           1 :             CATCH_REQUIRE(c[3][0] == r[3][0]);
    3872           1 :             CATCH_REQUIRE(c[3][1] == r[3][1]);
    3873           1 :             CATCH_REQUIRE(c[3][2] == r[3][2]);
    3874           1 :             CATCH_REQUIRE(c[3][3] == r[3][3]);
    3875             : 
    3876           1 :             double const determinant(b.determinant());
    3877           2 :             snapdev::matrix<double> adjugate(b.adjugate());
    3878             : 
    3879           2 :             snapdev::matrix<double> inv(adjugate * (1.0 / determinant));
    3880           2 :             snapdev::matrix<double> div(a * inv);
    3881             : 
    3882           1 :             CATCH_REQUIRE(fabs(c[0][0] - div[0][0]) < 0.0001);
    3883           1 :             CATCH_REQUIRE(fabs(c[0][1] - div[0][1]) < 0.0001);
    3884           1 :             CATCH_REQUIRE(fabs(c[0][2] - div[0][2]) < 0.0001);
    3885           1 :             CATCH_REQUIRE(fabs(c[0][3] - div[0][3]) < 0.0001);
    3886           1 :             CATCH_REQUIRE(fabs(c[1][0] - div[1][0]) < 0.0001);
    3887           1 :             CATCH_REQUIRE(fabs(c[1][1] - div[1][1]) < 0.0001);
    3888           1 :             CATCH_REQUIRE(fabs(c[1][2] - div[1][2]) < 0.0001);
    3889           1 :             CATCH_REQUIRE(fabs(c[1][3] - div[1][3]) < 0.0001);
    3890           1 :             CATCH_REQUIRE(fabs(c[2][0] - div[2][0]) < 0.0001);
    3891           1 :             CATCH_REQUIRE(fabs(c[2][1] - div[2][1]) < 0.0001);
    3892           1 :             CATCH_REQUIRE(fabs(c[2][2] - div[2][2]) < 0.0001);
    3893           1 :             CATCH_REQUIRE(fabs(c[2][3] - div[2][3]) < 0.0001);
    3894           1 :             CATCH_REQUIRE(fabs(c[3][0] - div[3][0]) < 0.0001);
    3895           1 :             CATCH_REQUIRE(fabs(c[3][1] - div[3][1]) < 0.0001);
    3896           1 :             CATCH_REQUIRE(fabs(c[3][2] - div[3][2]) < 0.0001);
    3897           1 :             CATCH_REQUIRE(fabs(c[3][3] - div[3][3]) < 0.0001);
    3898             :         }
    3899             :         CATCH_END_SECTION()
    3900             : 
    3901           8 :         CATCH_START_SECTION("matrix: a/=b")
    3902             :         {
    3903             :             // setup A
    3904             :             //
    3905           2 :             snapdev::matrix<double> a(4, 4);
    3906             : 
    3907           1 :             CATCH_REQUIRE(a.rows() == 4);
    3908           1 :             CATCH_REQUIRE(a.columns() == 4);
    3909             : 
    3910           1 :             double const a00(frand());
    3911           1 :             double const a01(frand());
    3912           1 :             double const a02(frand());
    3913           1 :             double const a03(frand());
    3914           1 :             double const a10(frand());
    3915           1 :             double const a11(frand());
    3916           1 :             double const a12(frand());
    3917           1 :             double const a13(frand());
    3918           1 :             double const a20(frand());
    3919           1 :             double const a21(frand());
    3920           1 :             double const a22(frand());
    3921           1 :             double const a23(frand());
    3922           1 :             double const a30(frand());
    3923           1 :             double const a31(frand());
    3924           1 :             double const a32(frand());
    3925           1 :             double const a33(frand());
    3926             : 
    3927           1 :             a[0][0] = a00;
    3928           1 :             a[0][1] = a01;
    3929           1 :             a[0][2] = a02;
    3930           1 :             a[0][3] = a03;
    3931           1 :             a[1][0] = a10;
    3932           1 :             a[1][1] = a11;
    3933           1 :             a[1][2] = a12;
    3934           1 :             a[1][3] = a13;
    3935           1 :             a[2][0] = a20;
    3936           1 :             a[2][1] = a21;
    3937           1 :             a[2][2] = a22;
    3938           1 :             a[2][3] = a23;
    3939           1 :             a[3][0] = a30;
    3940           1 :             a[3][1] = a31;
    3941           1 :             a[3][2] = a32;
    3942           1 :             a[3][3] = a33;
    3943             : 
    3944           1 :             CATCH_REQUIRE(a[0][0] == a00);
    3945           1 :             CATCH_REQUIRE(a[0][1] == a01);
    3946           1 :             CATCH_REQUIRE(a[0][2] == a02);
    3947           1 :             CATCH_REQUIRE(a[0][3] == a03);
    3948           1 :             CATCH_REQUIRE(a[1][0] == a10);
    3949           1 :             CATCH_REQUIRE(a[1][1] == a11);
    3950           1 :             CATCH_REQUIRE(a[1][2] == a12);
    3951           1 :             CATCH_REQUIRE(a[1][3] == a13);
    3952           1 :             CATCH_REQUIRE(a[2][0] == a20);
    3953           1 :             CATCH_REQUIRE(a[2][1] == a21);
    3954           1 :             CATCH_REQUIRE(a[2][2] == a22);
    3955           1 :             CATCH_REQUIRE(a[2][3] == a23);
    3956           1 :             CATCH_REQUIRE(a[3][0] == a30);
    3957           1 :             CATCH_REQUIRE(a[3][1] == a31);
    3958           1 :             CATCH_REQUIRE(a[3][2] == a32);
    3959           1 :             CATCH_REQUIRE(a[3][3] == a33);
    3960             : 
    3961             :             // setup B
    3962             :             //
    3963           2 :             snapdev::matrix<double> b(4, 4);
    3964             : 
    3965           1 :             CATCH_REQUIRE(b.rows() == 4);
    3966           1 :             CATCH_REQUIRE(b.columns() == 4);
    3967             : 
    3968           1 :             double const b00(frand());
    3969           1 :             double const b01(frand());
    3970           1 :             double const b02(frand());
    3971           1 :             double const b03(frand());
    3972           1 :             double const b10(frand());
    3973           1 :             double const b11(frand());
    3974           1 :             double const b12(frand());
    3975           1 :             double const b13(frand());
    3976           1 :             double const b20(frand());
    3977           1 :             double const b21(frand());
    3978           1 :             double const b22(frand());
    3979           1 :             double const b23(frand());
    3980           1 :             double const b30(frand());
    3981           1 :             double const b31(frand());
    3982           1 :             double const b32(frand());
    3983           1 :             double const b33(frand());
    3984             : 
    3985           1 :             b[0][0] = b00;
    3986           1 :             b[0][1] = b01;
    3987           1 :             b[0][2] = b02;
    3988           1 :             b[0][3] = b03;
    3989           1 :             b[1][0] = b10;
    3990           1 :             b[1][1] = b11;
    3991           1 :             b[1][2] = b12;
    3992           1 :             b[1][3] = b13;
    3993           1 :             b[2][0] = b20;
    3994           1 :             b[2][1] = b21;
    3995           1 :             b[2][2] = b22;
    3996           1 :             b[2][3] = b23;
    3997           1 :             b[3][0] = b30;
    3998           1 :             b[3][1] = b31;
    3999           1 :             b[3][2] = b32;
    4000           1 :             b[3][3] = b33;
    4001             : 
    4002           1 :             CATCH_REQUIRE(b[0][0] == b00);
    4003           1 :             CATCH_REQUIRE(b[0][1] == b01);
    4004           1 :             CATCH_REQUIRE(b[0][2] == b02);
    4005           1 :             CATCH_REQUIRE(b[0][3] == b03);
    4006           1 :             CATCH_REQUIRE(b[1][0] == b10);
    4007           1 :             CATCH_REQUIRE(b[1][1] == b11);
    4008           1 :             CATCH_REQUIRE(b[1][2] == b12);
    4009           1 :             CATCH_REQUIRE(b[1][3] == b13);
    4010           1 :             CATCH_REQUIRE(b[2][0] == b20);
    4011           1 :             CATCH_REQUIRE(b[2][1] == b21);
    4012           1 :             CATCH_REQUIRE(b[2][2] == b22);
    4013           1 :             CATCH_REQUIRE(b[2][3] == b23);
    4014           1 :             CATCH_REQUIRE(b[3][0] == b30);
    4015           1 :             CATCH_REQUIRE(b[3][1] == b31);
    4016           1 :             CATCH_REQUIRE(b[3][2] == b32);
    4017           1 :             CATCH_REQUIRE(b[3][3] == b33);
    4018             : 
    4019           1 :             double const determinant(b.determinant());
    4020           2 :             snapdev::matrix<double> adjugate(b.adjugate());
    4021             : 
    4022           2 :             snapdev::matrix<double> inv(adjugate * (1.0 / determinant));
    4023           2 :             snapdev::matrix<double> div(a * inv);
    4024             : 
    4025             :             // run operation A /= B
    4026             :             //
    4027           1 :             a /= b;
    4028             : 
    4029             :             // this could fail if one of bXX is 0.0
    4030             :             //
    4031             :             //CATCH_REQUIRE(a[0][0] == a00);
    4032             :             //CATCH_REQUIRE(a[0][1] == a01);
    4033             :             //CATCH_REQUIRE(a[0][2] == a02);
    4034             :             //CATCH_REQUIRE(a[0][3] == a03);
    4035             :             //CATCH_REQUIRE(a[1][0] == a10);
    4036             :             //CATCH_REQUIRE(a[1][1] == a11);
    4037             :             //CATCH_REQUIRE(a[1][2] == a12);
    4038             :             //CATCH_REQUIRE(a[1][3] == a13);
    4039             :             //CATCH_REQUIRE(a[2][0] == a20);
    4040             :             //CATCH_REQUIRE(a[2][1] == a21);
    4041             :             //CATCH_REQUIRE(a[2][2] == a22);
    4042             :             //CATCH_REQUIRE(a[2][3] == a23);
    4043             :             //CATCH_REQUIRE(a[3][0] == a30);
    4044             :             //CATCH_REQUIRE(a[3][1] == a31);
    4045             :             //CATCH_REQUIRE(a[3][2] == a32);
    4046             :             //CATCH_REQUIRE(a[3][3] == a33);
    4047             : 
    4048           1 :             CATCH_REQUIRE(b[0][0] == b00);
    4049           1 :             CATCH_REQUIRE(b[0][1] == b01);
    4050           1 :             CATCH_REQUIRE(b[0][2] == b02);
    4051           1 :             CATCH_REQUIRE(b[0][3] == b03);
    4052           1 :             CATCH_REQUIRE(b[1][0] == b10);
    4053           1 :             CATCH_REQUIRE(b[1][1] == b11);
    4054           1 :             CATCH_REQUIRE(b[1][2] == b12);
    4055           1 :             CATCH_REQUIRE(b[1][3] == b13);
    4056           1 :             CATCH_REQUIRE(b[2][0] == b20);
    4057           1 :             CATCH_REQUIRE(b[2][1] == b21);
    4058           1 :             CATCH_REQUIRE(b[2][2] == b22);
    4059           1 :             CATCH_REQUIRE(b[2][3] == b23);
    4060           1 :             CATCH_REQUIRE(b[3][0] == b30);
    4061           1 :             CATCH_REQUIRE(b[3][1] == b31);
    4062           1 :             CATCH_REQUIRE(b[3][2] == b32);
    4063           1 :             CATCH_REQUIRE(b[3][3] == b33);
    4064             : 
    4065           1 :             CATCH_REQUIRE(fabs(a[0][0] - div[0][0]) < 0.0001);
    4066           1 :             CATCH_REQUIRE(fabs(a[0][1] - div[0][1]) < 0.0001);
    4067           1 :             CATCH_REQUIRE(fabs(a[0][2] - div[0][2]) < 0.0001);
    4068           1 :             CATCH_REQUIRE(fabs(a[0][3] - div[0][3]) < 0.0001);
    4069           1 :             CATCH_REQUIRE(fabs(a[1][0] - div[1][0]) < 0.0001);
    4070           1 :             CATCH_REQUIRE(fabs(a[1][1] - div[1][1]) < 0.0001);
    4071           1 :             CATCH_REQUIRE(fabs(a[1][2] - div[1][2]) < 0.0001);
    4072           1 :             CATCH_REQUIRE(fabs(a[1][3] - div[1][3]) < 0.0001);
    4073           1 :             CATCH_REQUIRE(fabs(a[2][0] - div[2][0]) < 0.0001);
    4074           1 :             CATCH_REQUIRE(fabs(a[2][1] - div[2][1]) < 0.0001);
    4075           1 :             CATCH_REQUIRE(fabs(a[2][2] - div[2][2]) < 0.0001);
    4076           1 :             CATCH_REQUIRE(fabs(a[2][3] - div[2][3]) < 0.0001);
    4077           1 :             CATCH_REQUIRE(fabs(a[3][0] - div[3][0]) < 0.0001);
    4078           1 :             CATCH_REQUIRE(fabs(a[3][1] - div[3][1]) < 0.0001);
    4079           1 :             CATCH_REQUIRE(fabs(a[3][2] - div[3][2]) < 0.0001);
    4080           1 :             CATCH_REQUIRE(fabs(a[3][3] - div[3][3]) < 0.0001);
    4081             :         }
    4082             :         CATCH_END_SECTION()
    4083             :     }
    4084           8 : }
    4085             : 
    4086             : 
    4087           5 : CATCH_TEST_CASE("matrix_color", "[matrix]")
    4088             : {
    4089             :     // generate a few brightness matrices
    4090             :     //
    4091           6 :     CATCH_GIVEN("brightness")
    4092             :     {
    4093           2 :         CATCH_START_SECTION("matrix: b=a*<brightness> (a is identity)")
    4094             :         {
    4095             :             // setup A
    4096             :             //
    4097           2 :             snapdev::matrix<double> a(4, 4);
    4098             : 
    4099           1 :             CATCH_REQUIRE(a.rows() == 4);
    4100           1 :             CATCH_REQUIRE(a.columns() == 4);
    4101             : 
    4102             :             // setup B
    4103             :             //
    4104           2 :             snapdev::matrix<double> b(4, 4);
    4105             : 
    4106           1 :             CATCH_REQUIRE(b.rows() == 4);
    4107           1 :             CATCH_REQUIRE(b.columns() == 4);
    4108             : 
    4109        1001 :             for(int idx(0); idx < 1000; ++idx)
    4110             :             {
    4111        1000 :                 double const brightness(static_cast<double>(idx) / 1000.0);
    4112        1000 :                 b = a.brightness(brightness);
    4113             : 
    4114        1000 :                 CATCH_REQUIRE(fabs(b[0][0] - brightness) < 0.0001);
    4115        1000 :                 CATCH_REQUIRE(fabs(b[0][1] - 0.0) < 0.0001);
    4116        1000 :                 CATCH_REQUIRE(fabs(b[0][2] - 0.0) < 0.0001);
    4117        1000 :                 CATCH_REQUIRE(fabs(b[0][3] - 0.0) < 0.0001);
    4118        1000 :                 CATCH_REQUIRE(fabs(b[1][0] - 0.0) < 0.0001);
    4119        1000 :                 CATCH_REQUIRE(fabs(b[1][1] - brightness) < 0.0001);
    4120        1000 :                 CATCH_REQUIRE(fabs(b[1][2] - 0.0) < 0.0001);
    4121        1000 :                 CATCH_REQUIRE(fabs(b[1][3] - 0.0) < 0.0001);
    4122        1000 :                 CATCH_REQUIRE(fabs(b[2][0] - 0.0) < 0.0001);
    4123        1000 :                 CATCH_REQUIRE(fabs(b[2][1] - 0.0) < 0.0001);
    4124        1000 :                 CATCH_REQUIRE(fabs(b[2][2] - brightness) < 0.0001);
    4125        1000 :                 CATCH_REQUIRE(fabs(b[2][3] - 0.0) < 0.0001);
    4126        1000 :                 CATCH_REQUIRE(fabs(b[3][0] - 0.0) < 0.0001);
    4127        1000 :                 CATCH_REQUIRE(fabs(b[3][1] - 0.0) < 0.0001);
    4128        1000 :                 CATCH_REQUIRE(fabs(b[3][2] - 0.0) < 0.0001);
    4129        1000 :                 CATCH_REQUIRE(fabs(b[3][3] - 1.0) < 0.0001);
    4130             : 
    4131             : //     *       brightness_r & 0 & 0 & 0
    4132             : //     *    \\ 0 & brightness_g & 0 & 0
    4133             : //     *    \\ 0 & 0 & brightness_b & 0
    4134             : //     *    \\ 0 & 0 & 0 & 1
    4135             : 
    4136             :             }
    4137             :         }
    4138             :         CATCH_END_SECTION()
    4139             :     }
    4140             : 
    4141             :     // generate a few saturation matrices
    4142             :     //
    4143           6 :     CATCH_GIVEN("saturation")
    4144             :     {
    4145           2 :         CATCH_START_SECTION("matrix: b=a*<saturation> (a is identity)")
    4146             :         {
    4147             :             // setup A
    4148             :             //
    4149           2 :             snapdev::matrix<double> a(4, 4);
    4150             : 
    4151           1 :             CATCH_REQUIRE(a.rows() == 4);
    4152           1 :             CATCH_REQUIRE(a.columns() == 4);
    4153             : 
    4154             :             // setup B
    4155             :             //
    4156           2 :             snapdev::matrix<double> b(4, 4);
    4157             : 
    4158           1 :             CATCH_REQUIRE(b.rows() == 4);
    4159           1 :             CATCH_REQUIRE(b.columns() == 4);
    4160             : 
    4161           1 :             double luma[6][3];
    4162             : 
    4163           1 :             luma[0][0] = a.HDTV_LUMA_RED  ;
    4164           1 :             luma[0][1] = a.HDTV_LUMA_GREEN;
    4165           1 :             luma[0][2] = a.HDTV_LUMA_BLUE ;
    4166             : 
    4167           1 :             luma[1][0] = a.LED_LUMA_RED  ;
    4168           1 :             luma[1][1] = a.LED_LUMA_GREEN;
    4169           1 :             luma[1][2] = a.LED_LUMA_BLUE ;
    4170             : 
    4171           1 :             luma[2][0] = a.CRT_LUMA_RED  ;
    4172           1 :             luma[2][1] = a.CRT_LUMA_GREEN;
    4173           1 :             luma[2][2] = a.CRT_LUMA_BLUE ;
    4174             : 
    4175           1 :             luma[3][0] = a.NTSC_LUMA_RED  ;
    4176           1 :             luma[3][1] = a.NTSC_LUMA_GREEN;
    4177           1 :             luma[3][2] = a.NTSC_LUMA_BLUE ;
    4178             : 
    4179           1 :             luma[4][0] = a.AVERAGE_LUMA_RED  ;
    4180           1 :             luma[4][1] = a.AVERAGE_LUMA_GREEN;
    4181           1 :             luma[4][2] = a.AVERAGE_LUMA_BLUE ;
    4182             : 
    4183           1 :             luma[5][0] = 0.2; // to test any dynamic version if such exists for this algorithm
    4184           1 :             luma[5][1] = 0.7;
    4185           1 :             luma[5][2] = 0.1;
    4186             : 
    4187           7 :             for(int l(0); l < 6; ++l)
    4188             :             {
    4189             :                 // set the luma weights
    4190             :                 //
    4191           6 :                 a.set_luma_vector(luma[l][0], luma[l][1], luma[l][2]);
    4192             : 
    4193        6006 :                 for(int idx(0); idx < 1000; ++idx)
    4194             :                 {
    4195        6000 :                     double const s(static_cast<double>(idx) / 1000.0);
    4196        6000 :                     b = a.saturation(s);
    4197             : 
    4198        6000 :                     CATCH_REQUIRE(fabs(b[0][0] - (luma[l][0] * (1.0 - s) + s)) < 0.0001);
    4199        6000 :                     CATCH_REQUIRE(fabs(b[0][1] - (luma[l][0] * (1.0 - s))) < 0.0001);
    4200        6000 :                     CATCH_REQUIRE(fabs(b[0][2] - (luma[l][0] * (1.0 - s))) < 0.0001);
    4201        6000 :                     CATCH_REQUIRE(fabs(b[0][3] - (0.0)) < 0.0001);
    4202        6000 :                     CATCH_REQUIRE(fabs(b[1][0] - (luma[l][1] * (1.0 - s))) < 0.0001);
    4203        6000 :                     CATCH_REQUIRE(fabs(b[1][1] - (luma[l][1] * (1.0 - s) + s)) < 0.0001);
    4204        6000 :                     CATCH_REQUIRE(fabs(b[1][2] - (luma[l][1] * (1.0 - s))) < 0.0001);
    4205        6000 :                     CATCH_REQUIRE(fabs(b[1][3] - (0.0)) < 0.0001);
    4206        6000 :                     CATCH_REQUIRE(fabs(b[2][0] - (luma[l][2] * (1.0 - s))) < 0.0001);
    4207        6000 :                     CATCH_REQUIRE(fabs(b[2][1] - (luma[l][2] * (1.0 - s))) < 0.0001);
    4208        6000 :                     CATCH_REQUIRE(fabs(b[2][2] - (luma[l][2] * (1.0 - s) + s)) < 0.0001);
    4209        6000 :                     CATCH_REQUIRE(fabs(b[2][3] - (0.0)) < 0.0001);
    4210        6000 :                     CATCH_REQUIRE(fabs(b[3][0] - (0.0)) < 0.0001);
    4211        6000 :                     CATCH_REQUIRE(fabs(b[3][1] - (0.0)) < 0.0001);
    4212        6000 :                     CATCH_REQUIRE(fabs(b[3][2] - (0.0)) < 0.0001);
    4213        6000 :                     CATCH_REQUIRE(fabs(b[3][3] - (1.0)) < 0.0001);
    4214             :                 }
    4215             :             }
    4216             :         }
    4217             :         CATCH_END_SECTION()
    4218             :     }
    4219             : 
    4220             :     // generate a few brightness matrices
    4221             :     //
    4222           6 :     CATCH_GIVEN("hue")
    4223             :     {
    4224           2 :         CATCH_START_SECTION("matrix: b=a*<hue> (a is identity)")
    4225             :         {
    4226             :             // setup A
    4227             :             //
    4228           2 :             snapdev::matrix<double> a(4, 4);
    4229             : 
    4230           1 :             CATCH_REQUIRE(a.rows() == 4);
    4231           1 :             CATCH_REQUIRE(a.columns() == 4);
    4232             : 
    4233             :             // setup B
    4234             :             //
    4235           2 :             snapdev::matrix<double> b(4, 4);
    4236             : 
    4237           1 :             CATCH_REQUIRE(b.rows() == 4);
    4238           1 :             CATCH_REQUIRE(b.columns() == 4);
    4239             : 
    4240           1 :             double luma[6][3];
    4241             : 
    4242           1 :             luma[0][0] = a.HDTV_LUMA_RED  ;
    4243           1 :             luma[0][1] = a.HDTV_LUMA_GREEN;
    4244           1 :             luma[0][2] = a.HDTV_LUMA_BLUE ;
    4245             : 
    4246           1 :             luma[1][0] = a.LED_LUMA_RED  ;
    4247           1 :             luma[1][1] = a.LED_LUMA_GREEN;
    4248           1 :             luma[1][2] = a.LED_LUMA_BLUE ;
    4249             : 
    4250           1 :             luma[2][0] = a.CRT_LUMA_RED  ;
    4251           1 :             luma[2][1] = a.CRT_LUMA_GREEN;
    4252           1 :             luma[2][2] = a.CRT_LUMA_BLUE ;
    4253             : 
    4254           1 :             luma[3][0] = a.NTSC_LUMA_RED  ;
    4255           1 :             luma[3][1] = a.NTSC_LUMA_GREEN;
    4256           1 :             luma[3][2] = a.NTSC_LUMA_BLUE ;
    4257             : 
    4258           1 :             luma[4][0] = a.AVERAGE_LUMA_RED  ;
    4259           1 :             luma[4][1] = a.AVERAGE_LUMA_GREEN;
    4260           1 :             luma[4][2] = a.AVERAGE_LUMA_BLUE ;
    4261             : 
    4262           1 :             luma[5][0] = 0.2; // to test any dynamic version if such exists for this algorithm
    4263           1 :             luma[5][1] = 0.7;
    4264           1 :             luma[5][2] = 0.1;
    4265             : 
    4266           1 :             std::vector<std::string> last_hue_matrix{
    4267             :                 "hdtv",
    4268             :                 "led",
    4269             :                 "crt",
    4270             :                 "ntsc",
    4271             :                 "average",
    4272             :                 "dynamic"
    4273           2 :             };
    4274             : 
    4275           7 :             for(int l(0); l < 6; ++l)
    4276             :             {
    4277             :                 // set the luma weights
    4278             :                 //
    4279           6 :                 a.set_luma_vector(luma[l][0], luma[l][1], luma[l][2]);
    4280             : 
    4281             : //std::cerr << "*** working on luma " << l << "\n";
    4282        6006 :                 for(int idx(0); idx < 1000; ++idx)
    4283             :                 {
    4284             :                     // hue is an angle from 0 to 360 (albeit in radian)
    4285             :                     //
    4286        6000 :                     double const hue(static_cast<double>(idx) / 1000.0 * M_PI * 2.0);
    4287        6000 :                     b = a.hue(hue);
    4288             : 
    4289             : #ifdef _DEBUG
    4290        6000 :                     CATCH_REQUIRE(a.get_last_hue_matrix() == last_hue_matrix[l]);
    4291             : #endif
    4292             : 
    4293             :                     // this one requires us to work!
    4294             :                     // here we calculate the hue matrix by hand instead of using
    4295             :                     // pre-calculated numbers; that way we know we have it right
    4296             :                     //
    4297             : 
    4298             :                     // $R_r$ (red rotation)
    4299             :                     //
    4300       12000 :                     snapdev::matrix<double> r_r(4, 4);
    4301        6000 :                     r_r[1][1] =  1.0 / sqrt(2.0);
    4302        6000 :                     r_r[1][2] =  1.0 / sqrt(2.0);
    4303        6000 :                     r_r[2][1] = -1.0 / sqrt(2.0);
    4304        6000 :                     r_r[2][2] =  1.0 / sqrt(2.0);
    4305             : 
    4306             :                     // $R_g$ (green rotation)
    4307             :                     //
    4308       12000 :                     snapdev::matrix<double> r_g(4, 4);
    4309        6000 :                     r_g[0][0] =  sqrt(2.0) / sqrt(3.0);
    4310        6000 :                     r_g[0][2] =  1.0 / sqrt(3.0);
    4311        6000 :                     r_g[2][0] = -1.0 / sqrt(3.0);
    4312        6000 :                     r_g[2][2] =  sqrt(2.0) / sqrt(3.0);
    4313             : 
    4314             :                     // $R_rg = R_r R_g$ (red and green rotations combined)
    4315             :                     //
    4316       12000 :                     snapdev::matrix<double> r_rg(r_r * r_g);
    4317             : 
    4318             :                     // $w$ (luma vector, a.k.a. color weights to calcalute the perfect grayscale for your monitor)
    4319             :                     //
    4320       12000 :                     snapdev::matrix<double> w(a.get_luma_vector());
    4321             :                     //w[0][0] = a.RED_WEIGHT;
    4322             :                     //w[1][0] = a.GREEN_WEIGHT;
    4323             :                     //w[2][0] = a.BLUE_WEIGHT;
    4324             :                     //w[3][0] = 0;
    4325             : 
    4326             :                     // sm (skew matrix)
    4327             :                     //
    4328       12000 :                     snapdev::matrix<double> sm(r_rg * w);
    4329             : 
    4330             :                     // s (scaling with skew matrix)
    4331             :                     //
    4332       12000 :                     snapdev::matrix<double> s(4, 4);
    4333        6000 :                     s[0][2] = sm[0][0] / sm[2][0];
    4334        6000 :                     s[1][2] = sm[1][0] / sm[2][0];
    4335             : 
    4336             :                     // p (skewed red & green rotation)
    4337             :                     //
    4338       12000 :                     snapdev::matrix<double> p(r_rg);
    4339        6000 :                     p *= s;
    4340             : 
    4341             :                     // r_b (blue rotation)
    4342             :                     //
    4343       12000 :                     snapdev::matrix<double> r_b(4, 4);
    4344        6000 :                     double const rot_cos(cos(hue));
    4345        6000 :                     double const rot_sin(sin(hue));
    4346        6000 :                     r_b[0][0] =  rot_cos;
    4347        6000 :                     r_b[0][1] =  rot_sin;
    4348        6000 :                     r_b[1][0] = -rot_sin;
    4349        6000 :                     r_b[1][1] =  rot_cos;
    4350             : 
    4351             :                     // hue_matrix (the resulting matrix depending on angle)
    4352             :                     //
    4353       12000 :                     snapdev::matrix<double> hue_matrix(p * r_b / p);
    4354             : 
    4355             : //std::cout << "Got b = " << b << "\n";
    4356             : //std::cout << "Got our hue matrix = " << hue_matrix << "\n";
    4357             : 
    4358        6000 :                     CATCH_REQUIRE(fabs(b[0][0] - hue_matrix[0][0]) < 0.0001);
    4359        6000 :                     CATCH_REQUIRE(fabs(b[0][1] - hue_matrix[0][1]) < 0.0001);
    4360        6000 :                     CATCH_REQUIRE(fabs(b[0][2] - hue_matrix[0][2]) < 0.0001);
    4361        6000 :                     CATCH_REQUIRE(fabs(b[0][3] - hue_matrix[0][3]) < 0.0001);
    4362        6000 :                     CATCH_REQUIRE(fabs(b[1][0] - hue_matrix[1][0]) < 0.0001);
    4363        6000 :                     CATCH_REQUIRE(fabs(b[1][1] - hue_matrix[1][1]) < 0.0001);
    4364        6000 :                     CATCH_REQUIRE(fabs(b[1][2] - hue_matrix[1][2]) < 0.0001);
    4365        6000 :                     CATCH_REQUIRE(fabs(b[1][3] - hue_matrix[1][3]) < 0.0001);
    4366        6000 :                     CATCH_REQUIRE(fabs(b[2][0] - hue_matrix[2][0]) < 0.0001);
    4367        6000 :                     CATCH_REQUIRE(fabs(b[2][1] - hue_matrix[2][1]) < 0.0001);
    4368        6000 :                     CATCH_REQUIRE(fabs(b[2][2] - hue_matrix[2][2]) < 0.0001);
    4369        6000 :                     CATCH_REQUIRE(fabs(b[2][3] - hue_matrix[2][3]) < 0.0001);
    4370        6000 :                     CATCH_REQUIRE(fabs(b[3][0] - hue_matrix[3][0]) < 0.0001);
    4371        6000 :                     CATCH_REQUIRE(fabs(b[3][1] - hue_matrix[3][1]) < 0.0001);
    4372        6000 :                     CATCH_REQUIRE(fabs(b[3][2] - hue_matrix[3][2]) < 0.0001);
    4373        6000 :                     CATCH_REQUIRE(fabs(b[3][3] - hue_matrix[3][3]) < 0.0001);
    4374             : 
    4375             : //     * $$
    4376             : //     * R_r =
    4377             : //     * \begin{bmatrix}
    4378             : //     *       1 &  0                  & 0                  & 0
    4379             : //     *    \\ 0 &  {1 \over \sqrt 2 } & {1 \over \sqrt 2 } & 0
    4380             : //     *    \\ 0 & -{1 \over \sqrt 2 } & {1 \over \sqrt 2 } & 0
    4381             : //     *    \\ 0 &  0                  & 0                  & 1
    4382             : //     * \end{bmatrix}
    4383             : //     * $$
    4384             : //     *
    4385             : //     * $$
    4386             : //     * R_g =
    4387             : //     * \begin{bmatrix}
    4388             : //     *       1 & 0 & 0 & 0
    4389             : //     *    \\ 0 & {\sqrt 2 \over \sqrt 3} & -{1 \over \sqrt 3} & 0
    4390             : //     *    \\ 0 & {1 \over \sqrt 3} & {\sqrt 2 \over \sqrt 3} & 0
    4391             : //     *    \\ 0 & 0 & 0 & 1
    4392             : //     * \end{bmatrix}
    4393             : //     * $$
    4394             : //     *
    4395             : //     * $$
    4396             : //     * R_{rg} = R_r R_g
    4397             : //     * $$
    4398             : //     *
    4399             : //     * $$
    4400             : //     * \begin{bmatrix}
    4401             : //     *        L_r
    4402             : //     *     \\ L_g
    4403             : //     *     \\ L_b
    4404             : //     *     \\ 0
    4405             : //     * \end{bmatrix}
    4406             : //     * =
    4407             : //     * R_{rg}
    4408             : //     * \begin{bmatrix}
    4409             : //     *        W_r
    4410             : //     *     \\ W_g
    4411             : //     *     \\ W_b
    4412             : //     *     \\ 0
    4413             : //     * \end{bmatrix}
    4414             : //     * $$
    4415             : //     *
    4416             : //     * Now we define a rotation matrix for the blue axis. This one uses
    4417             : //     * a variable as the angle. This is the actual rotation angle which
    4418             : //     * can go from 0 to $2 \pi$:
    4419             : //     *
    4420             : //     * $$
    4421             : //     * R_b =
    4422             : //     * \begin{bmatrix}
    4423             : //     *       cos( \theta )  & sin( \theta ) & 0 & 0
    4424             : //     *    \\ -sin( \theta ) & cos( \theta ) & 0 & 0
    4425             : //     *    \\ 0              & 0             & 1 & 0
    4426             : //     *    \\ 0              & 0             & 0 & 1
    4427             : //     * \end{bmatrix}
    4428             : //     * $$
    4429             : //     *
    4430             : //     * Now we have all the matrices to caculate the hue rotation
    4431             : //     * of all the components of an image:
    4432             : //     *
    4433             : //     * $$
    4434             : //     * H = R_{rg} S R_b S^{-1} R_{rg}^{-1}
    4435             : //     * $$
    4436             : //     *
    4437             : //     * $H$ can then be used as in:
    4438             : //     *
    4439             : //     * $$
    4440             : //     * \begin{bmatrix}
    4441             : //     *      R'
    4442             : //     *   \\ G'
    4443             : //     *   \\ B'
    4444             : //     * \end{bmatrix}
    4445             : //     * =
    4446             : //     * H
    4447             : //     * \begin{bmatrix}
    4448             : //     *      R
    4449             : //     *   \\ G
    4450             : //     *   \\ B
    4451             : //     * \end{bmatrix}
    4452             : //     * $$
    4453             : //     *
    4454             :                 }
    4455             :             }
    4456             :         }
    4457             :         CATCH_END_SECTION()
    4458             :     }
    4459           9 : }
    4460             : 
    4461             : 
    4462             : // vim: ts=4 sw=4 et

Generated by: LCOV version 1.13