Justine Devin / gammalib

Commit 453e69321e92054568041cee52cf15f40c276321

Authored by knodlseder
1 parent b3bbe291

start code cleaning

Inline Side-by-side
Showing 1 changed file with 175 additions and 150 deletions
test/test_GSparseMatrix.cpp
  View file @453e693
1 1 /***************************************************************************
2 2 * test_GSparseMatrix.cpp - test GSparseMatrix class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright : (C) 2006 by Jurgen Knodlseder *
  4 + * copyright (C) 2006-2010 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software; you can redistribute it and/or modify *
... ... @@ -10,8 +10,14 @@
10 10 * (at your option) any later version. *
11 11 * *
12 12 ***************************************************************************/
  13 +/**
  14 + * @file test_GSparseMatrix.cpp
  15 + * @brief Testing of GSparseMatrix class
  16 + * @author J. Knodlseder
  17 + */
13 18  
14 19 /* __ Includes ___________________________________________________________ */
  20 +#include <stdlib.h>
15 21 #include "test_GSparseMatrix.hpp"
16 22  
17 23 /* __ Namespaces _________________________________________________________ */
... ... @@ -997,36 +1003,41 @@ void test_comparison(const GSparseMatrix&amp; m_test)
997 1003 ***************************************************************************/
998 1004 void test_conversion(void)
999 1005 {
1000   - cout << "Test GSparseMatrix: Matrix conversions: ";
1001   - try {
1002   - //
1003   - // Setup a symmetric sparse matrix
1004   - int num = 10;
1005   - GSparseMatrix sparse(num,num);
1006   - for (int i = 0; i < num; ++i)
1007   - sparse(i,i) = (i+1);
1008   - for (int j = 1; j < num; ++j) {
1009   - sparse(0,j) = (j+1);
1010   - sparse(j,0) = (j+1);
1011   - }
1012   - cout << ".";
1013   - //
1014   - // Convert GSparseMatrix matrix into GSymMatrix object
1015   - GSymMatrix converted = sym_matrix(sparse);
1016   - cout << ".";
1017   - //
1018   - // Convert GSymMatrix back to GSparseMatrix matrix
1019   - GSparseMatrix back_convert = sparse_matrix(converted);
1020   - cout << ".";
1021   - //
1022   - // Compare back converted matrix to original one. They should be identical
1023   - if (sparse != back_convert) {
1024   - cout << endl << "TEST ERROR: Unable to convert matrixes (symmetric)." << endl;
1025   - cout << "Original matrix " << sparse << endl;
1026   - cout << "GSymMatrix matrix " << converted << endl;
1027   - cout << "Back converted matrix " << back_convert << endl;
1028   - throw;
1029   - }
  1006 + // Dump header
  1007 + std::cout << "Test GSparseMatrix: Matrix conversions: ";
  1008 +
  1009 + // Test 1
  1010 + try {
  1011 + //
  1012 + // Setup a symmetric sparse matrix
  1013 + int num = 10;
  1014 + GSparseMatrix sparse(num,num);
  1015 + for (int i = 0; i < num; ++i)
  1016 + sparse(i,i) = (i+1);
  1017 + for (int j = 1; j < num; ++j) {
  1018 + sparse(0,j) = (j+1);
  1019 + sparse(j,0) = (j+1);
  1020 + }
  1021 + std::cout << ".";
  1022 + //
  1023 + // Convert GSparseMatrix matrix into GSymMatrix object
  1024 + GSymMatrix converted = sym_matrix(sparse);
  1025 + std::cout << ".";
  1026 + //
  1027 + // Convert GSymMatrix back to GSparseMatrix matrix
  1028 + GSparseMatrix back_convert = sparse_matrix(converted);
  1029 + std::cout << ".";
  1030 + //
  1031 + // Compare back converted matrix to original one. They should be identical
  1032 + if (sparse != back_convert) {
  1033 + std::cout << std::endl
  1034 + << "TEST ERROR: Unable to convert matrixes (symmetric)."
  1035 + << std::endl;
  1036 + std::cout << "Original matrix " << sparse << std::endl;
  1037 + std::cout << "GSymMatrix matrix " << converted << std::endl;
  1038 + std::cout << "Back converted matrix " << back_convert << std::endl;
  1039 + throw;
  1040 + }
1030 1041 //
1031 1042 // Determine the fill of the matrix. It should be 0.28
1032 1043 double fill = back_convert.fill();
... ... @@ -1136,70 +1147,76 @@ void test_conversion(void)
1136 1147 ***************************************************************************/
1137 1148 void test_cholesky(void)
1138 1149 {
1139   - cout << "Test GSparseMatrix: Cholesky decomposition, solver and inverter: ";
1140   - try {
1141   - //
1142   - // Setup matrix for Cholesky decomposition
1143   - GSparseMatrix m_chol_test(5,5);
1144   - m_chol_test(0,0) = 1.0;
1145   - m_chol_test(0,1) = 0.2;
1146   - m_chol_test(0,2) = 0.2;
1147   - m_chol_test(0,3) = 0.2;
1148   - m_chol_test(0,4) = 0.2;
1149   - m_chol_test(1,0) = 0.2;
1150   - m_chol_test(2,0) = 0.2;
1151   - m_chol_test(3,0) = 0.2;
1152   - m_chol_test(4,0) = 0.2;
1153   - m_chol_test(1,1) = 1.0;
1154   - m_chol_test(2,2) = 1.0;
1155   - m_chol_test(3,3) = 1.0;
1156   - m_chol_test(4,4) = 1.0;
1157   - //
1158   - // Try to solve now (should not work)
1159   - try {
1160   - GVector v_test12(5);
1161   - v_test12 = m_chol_test.cholesky_solver(v_test12);
1162   - }
1163   - catch (GException::matrix_not_factorised) {
1164   - cout << ".";
1165   - }
1166   - catch (exception &e) {
1167   - cout << e.what() << endl;
1168   - throw;
1169   - }
1170   - //
1171   - // Test Cholesky decomposition
1172   - GSparseMatrix cd = cholesky_decompose(m_chol_test);
1173   - cout << ".";
1174   - //
1175   - // Test inplace Cholesky decomposition
1176   - cd = m_chol_test;
1177   - cd.cholesky_decompose();
1178   - cout << ".";
1179   - //
1180   - // Test Cholesky solver
1181   - GVector e0(5);
1182   - GVector a0(5);
1183   - e0(0) = 1.0;
1184   - a0(0) = 1.0;
1185   - a0(1) = 0.2;
1186   - a0(2) = 0.2;
1187   - a0(3) = 0.2;
1188   - a0(4) = 0.2;
1189   - GVector s0 = cd.cholesky_solver(a0);
1190   - double res = max(fabs(s0-e0));
1191   - if (res < 1.0e-15)
1192   - #if defined(DUMP_RESIDUALS)
1193   - cout << " Res(S0)=" << res << " ";
1194   - #else
1195   - cout << ".";
1196   - #endif
1197   - else {
1198   - cout << endl << "TEST ERROR: Corrupt cholesky_solver(GVector) function." << endl;
1199   - cout << "Residual vector (all elements should be zero):" << endl;
1200   - cout << s0 << endl;
1201   - throw;
1202   - }
  1150 + // Dump header
  1151 + std::cout << "Test GSparseMatrix: Cholesky decomposition, solver and inverter: ";
  1152 +
  1153 + // Test 1
  1154 + try {
  1155 + //
  1156 + // Setup matrix for Cholesky decomposition
  1157 + GSparseMatrix m_chol_test(5,5);
  1158 + m_chol_test(0,0) = 1.0;
  1159 + m_chol_test(0,1) = 0.2;
  1160 + m_chol_test(0,2) = 0.2;
  1161 + m_chol_test(0,3) = 0.2;
  1162 + m_chol_test(0,4) = 0.2;
  1163 + m_chol_test(1,0) = 0.2;
  1164 + m_chol_test(2,0) = 0.2;
  1165 + m_chol_test(3,0) = 0.2;
  1166 + m_chol_test(4,0) = 0.2;
  1167 + m_chol_test(1,1) = 1.0;
  1168 + m_chol_test(2,2) = 1.0;
  1169 + m_chol_test(3,3) = 1.0;
  1170 + m_chol_test(4,4) = 1.0;
  1171 + //
  1172 + // Try to solve now (should not work)
  1173 + try {
  1174 + GVector v_test12(5);
  1175 + v_test12 = m_chol_test.cholesky_solver(v_test12);
  1176 + }
  1177 + catch (GException::matrix_not_factorised) {
  1178 + std::cout << ".";
  1179 + }
  1180 + catch (exception &e) {
  1181 + std::cout << e.what() << std::endl;
  1182 + throw;
  1183 + }
  1184 + //
  1185 + // Test Cholesky decomposition
  1186 + GSparseMatrix cd = cholesky_decompose(m_chol_test);
  1187 + std::cout << ".";
  1188 + //
  1189 + // Test inplace Cholesky decomposition
  1190 + cd = m_chol_test;
  1191 + cd.cholesky_decompose();
  1192 + std::cout << ".";
  1193 + //
  1194 + // Test Cholesky solver
  1195 + GVector e0(5);
  1196 + GVector a0(5);
  1197 + e0(0) = 1.0;
  1198 + a0(0) = 1.0;
  1199 + a0(1) = 0.2;
  1200 + a0(2) = 0.2;
  1201 + a0(3) = 0.2;
  1202 + a0(4) = 0.2;
  1203 + GVector s0 = cd.cholesky_solver(a0);
  1204 + double res = max(fabs(s0-e0));
  1205 + if (res < 1.0e-15) {
  1206 + #if defined(DUMP_RESIDUALS)
  1207 + std::cout << " Res(S0)=" << res << " ";
  1208 + #else
  1209 + std::cout << ".";
  1210 + #endif
  1211 + }
  1212 + else {
  1213 + std::cout << std::endl
  1214 + << "TEST ERROR: Corrupt cholesky_solver(GVector) function."
  1215 + << std::endl;
  1216 + std::cout << "Residual vector (all elements should be zero):" << std::endl;
  1217 + std::cout << s0 << std::endl;
  1218 + throw;
  1219 + }
1203 1220 //
1204 1221 e0(0) = 0.0;
1205 1222 e0(1) = 1.0;
... ... @@ -1569,36 +1586,44 @@ void test_cholesky(void)
1569 1586 throw;
1570 1587 }
1571 1588 //
1572   - // Test Cholesky inverter for compressed matrix
1573   - unit = GSparseMatrix(6,6);
1574   - unit(0,0) = 1.0;
1575   - unit(1,1) = 1.0;
1576   - unit(2,2) = 1.0;
1577   - unit(4,4) = 1.0;
1578   - unit(5,5) = 1.0;
1579   - GSparseMatrix m_chol_test_zero_inv = m_chol_test_zero;
1580   - m_chol_test_zero_inv.cholesky_invert();
1581   - GSparseMatrix ciz_product = m_chol_test_zero * m_chol_test_zero_inv;
1582   - GSparseMatrix ciz_residuals = ciz_product - unit;
1583   - res = (fabs(ciz_residuals)).max();
1584   - if (res < 1.0e-15)
1585   - #if defined(DUMP_RESIDUALS)
1586   - cout << "Res(CIZ)=" << res << " ";
1587   - #else
1588   - cout << ".";
1589   - #endif
1590   - else {
1591   - cout << endl << "TEST ERROR: Corrupt compressed GSparseMatrix.cholesky_invert() function." << endl;
1592   - cout << "Residual matrix (all elements should be zero):" << endl;
1593   - cout << ciz_residuals << endl;
1594   - throw;
1595   - }
1596   - }
1597   - catch (exception &e) {
1598   - cout << e.what() << endl;
1599   - throw;
1600   - }
1601   - cout << "ok." << endl;
  1589 + // Test Cholesky inverter for compressed matrix
  1590 + unit = GSparseMatrix(6,6);
  1591 + unit(0,0) = 1.0;
  1592 + unit(1,1) = 1.0;
  1593 + unit(2,2) = 1.0;
  1594 + unit(4,4) = 1.0;
  1595 + unit(5,5) = 1.0;
  1596 + GSparseMatrix m_chol_test_zero_inv = m_chol_test_zero;
  1597 + m_chol_test_zero_inv.cholesky_invert();
  1598 + GSparseMatrix ciz_product = m_chol_test_zero * m_chol_test_zero_inv;
  1599 + GSparseMatrix ciz_residuals = ciz_product - unit;
  1600 + res = (fabs(ciz_residuals)).max();
  1601 + if (res < 1.0e-15) {
  1602 + #if defined(DUMP_RESIDUALS)
  1603 + std::cout << "Res(CIZ)=" << res << " ";
  1604 + #else
  1605 + std::cout << ".";
  1606 + #endif
  1607 + }
  1608 + else {
  1609 + std::cout << std::endl
  1610 + << "TEST ERROR: Corrupt compressed GSparseMatrix.cholesky_invert()"
  1611 + " method." << std::endl;
  1612 + std::cout << "Residual matrix (all elements should be zero):" << std::endl;
  1613 + std::cout << ciz_residuals << std::endl;
  1614 + throw;
  1615 + }
  1616 + }
  1617 + catch (exception &e) {
  1618 + std::cout << e.what() << std::endl;
  1619 + throw;
  1620 + }
  1621 +
  1622 + // Dump final test ok
  1623 + std::cout << "ok." << std::endl;
  1624 +
  1625 + // Return
  1626 + return;
1602 1627 }
1603 1628  
1604 1629  
... ... @@ -1747,32 +1772,32 @@ void test_heavy(int block = 512)
1747 1772 ***************************************************************************/
1748 1773 int main(void)
1749 1774 {
1750   - // Dump header
1751   - cout << endl;
1752   - cout << "*******************************" << endl;
1753   - cout << "* GSparseMatrix class testing *" << endl;
1754   - cout << "*******************************" << endl;
  1775 + // Dump header
  1776 + std::cout << std::endl;
  1777 + std::cout << "*******************************" << std::endl;
  1778 + std::cout << "* GSparseMatrix class testing *" << std::endl;
  1779 + std::cout << "*******************************" << std::endl;
1755 1780  
1756   - // Set test matrix and vector
1757   - GSparseMatrix m_test = set_matrix();
1758   - GVector v_test = set_vector();
  1781 + // Set test matrix and vector
  1782 + GSparseMatrix m_test = set_matrix();
  1783 + GVector v_test = set_vector();
1759 1784  
1760   - // Execute the tests
1761   - test_output(m_test);
1762   - test_allocate();
1763   - test_assign_values(m_test);
1764   - test_copy_constructor(m_test);
1765   - test_assign_matrix(m_test);
1766   - test_transpose(m_test);
1767   - test_matrix_vector(m_test, v_test);
1768   - test_matrix_matrix(m_test);
1769   - test_arithmetics(m_test);
1770   - test_functions(m_test);
1771   - test_comparison(m_test);
1772   - test_conversion();
1773   - test_cholesky();
1774   - test_heavy();
  1785 + // Execute the tests
  1786 + test_output(m_test);
  1787 + test_allocate();
  1788 + test_assign_values(m_test);
  1789 + test_copy_constructor(m_test);
  1790 + test_assign_matrix(m_test);
  1791 + test_transpose(m_test);
  1792 + test_matrix_vector(m_test, v_test);
  1793 + test_matrix_matrix(m_test);
  1794 + test_arithmetics(m_test);
  1795 + test_functions(m_test);
  1796 + test_comparison(m_test);
  1797 + test_conversion();
  1798 + test_cholesky();
  1799 + test_heavy();
1775 1800  
1776   - // Return
1777   - return 0;
  1801 + // Return
  1802 + return 0;
1778 1803 }
... ...