LCOV - code coverage report
Current view: top level - tests - matrix.cpp (source / functions) Hit Total Coverage
Test: coverage.info Lines: 3251 3260 99.7 %
Date: 2021-08-21 10:14:39 Functions: 8 8 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13