1 /* -*- mode: C -*- */
2 /*
3 IGraph library.
4 Copyright (C) 2011-2012 Gabor Csardi <csardi.gabor@gmail.com>
5 334 Harvard st, Cambridge MA, 02139 USA
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
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 02110-1301 USA
21
22 */
23
24 #include <igraph.h>
25
26 #include "test_utilities.inc"
27
main()28 int main() {
29
30 int nodes = 10;
31 igraph_t tree;
32 igraph_matrix_t sto;
33 igraph_matrix_t hess;
34 igraph_matrix_complex_t evec1, evec2;
35 igraph_vector_complex_t eval1, eval2;
36 igraph_eigen_which_t which;
37 int i;
38
39 igraph_tree(&tree, nodes, /* children= */ 3, IGRAPH_TREE_UNDIRECTED);
40
41 igraph_matrix_init(&sto, nodes, nodes);
42 igraph_get_stochastic(&tree, &sto, /*column_wise=*/ 0);
43 igraph_matrix_transpose(&sto);
44
45 igraph_matrix_init(&hess, nodes, nodes);
46 igraph_lapack_dgehrd(&sto, 1, nodes, &hess);
47
48 igraph_matrix_complex_init(&evec1, 0, 0);
49 igraph_vector_complex_init(&eval1, 0);
50 which.pos = IGRAPH_EIGEN_ALL;
51 igraph_eigen_matrix(&sto, 0, 0, nodes, 0, IGRAPH_EIGEN_LAPACK, &which, 0, 0,
52 &eval1, &evec1);
53
54 igraph_matrix_complex_init(&evec2, 0, 0);
55 igraph_vector_complex_init(&eval2, 0);
56 igraph_eigen_matrix(&hess, 0, 0, nodes, 0, IGRAPH_EIGEN_LAPACK, &which, 0,
57 0, &eval2, &evec2);
58
59 for (i = 0; i < nodes; i++) {
60 igraph_real_t d = igraph_complex_abs(igraph_complex_sub(VECTOR(eval1)[i],
61 VECTOR(eval2)[i]));
62 if (d > 1e-14) {
63 printf("Difference: %g\n", d);
64 return 1;
65 }
66 }
67
68 igraph_matrix_complex_destroy(&evec2);
69 igraph_vector_complex_destroy(&eval2);
70
71 igraph_matrix_complex_destroy(&evec1);
72 igraph_vector_complex_destroy(&eval1);
73
74 igraph_matrix_destroy(&hess);
75 igraph_matrix_destroy(&sto);
76 igraph_destroy(&tree);
77
78 VERIFY_FINALLY_STACK();
79
80 return 0;
81 }
82