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_real_t triplets[] = { 1, 0, 1 / 4.0,       0, 1, 1 / 3.0,
32                                  2, 0, 1 / 4.0,       0, 2, 1 / 3.0,
33                                  3, 0, 1.0,         0, 3, 1 / 3.0,
34                                  4, 1, 1.0,         1, 4, 1 / 4.0,
35                                  5, 1, 1.0,         1, 5, 1 / 4.0,
36                                  6, 1, 1.0,         1, 6, 1 / 4.0,
37                                  7, 2, 1.0,         2, 7, 1 / 4.0,
38                                  8, 2, 1.0,         2, 8, 1 / 4.0,
39                                  9, 2, 1.0,         2, 9, 1 / 4.0
40                                };
41 
42     igraph_sparsemat_t mat;
43     int i, n = sizeof(triplets) / sizeof(igraph_real_t);
44     igraph_eigen_which_t which;
45     igraph_vector_complex_t values, values2;
46     igraph_matrix_complex_t vectors, vectors2;
47     igraph_matrix_t mat2;
48 
49     igraph_sparsemat_init(&mat, nodes, nodes, n / 3);
50     for (i = 0; i < n; i += 3) {
51         igraph_sparsemat_entry(&mat, triplets[i], triplets[i + 1], triplets[i + 2]);
52     }
53 
54     which.pos = IGRAPH_EIGEN_LM;
55     which.howmany = 1;
56 
57     igraph_vector_complex_init(&values, 0);
58     igraph_matrix_complex_init(&vectors, 0, 0);
59 
60     igraph_eigen_matrix(/*matrix=*/ 0, /*sparsemat=*/ &mat, /*fun=*/ 0,
61                                     nodes, /*extra=*/ 0, IGRAPH_EIGEN_LAPACK, &which,
62                                     /*options=*/ 0, /*storage=*/ 0, &values, &vectors);
63 
64     if (IGRAPH_REAL(MATRIX(vectors, 0, 0)) < 0) {
65         igraph_matrix_complex_scale(&vectors, igraph_complex(-1.0, -0.0 ));
66     }
67 
68     igraph_vector_complex_print(&values);
69     igraph_matrix_complex_print(&vectors);
70 
71     igraph_sparsemat_destroy(&mat);
72 
73     /* Calcualate all eigenvalues, using SM and LM and then check that they
74        are the same, in opposite order. We use a random matrix this time. */
75 
76     igraph_rng_seed(igraph_rng_default(), 42);
77     igraph_matrix_init(&mat2, nodes, nodes);
78     for (i = 0; i < nodes; i++) {
79         int j;
80         for (j = 0; j < nodes; j++) {
81             MATRIX(mat2, i, j) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
82         }
83     }
84 
85     which.pos = IGRAPH_EIGEN_LM;
86     which.howmany = nodes;
87     igraph_eigen_matrix(&mat2, /*sparsemat=*/ 0, /*fun=*/ 0, nodes,
88                         /*extra=*/ 0, IGRAPH_EIGEN_LAPACK, &which,
89                         /*options=*/ 0, /*storage=*/ 0, &values, &vectors);
90 
91     which.pos = IGRAPH_EIGEN_SM;
92     which.howmany = nodes;
93     igraph_vector_complex_init(&values2, 0);
94     igraph_matrix_complex_init(&vectors2, 0, 0);
95     igraph_eigen_matrix(&mat2, /*sparsemat=*/ 0, /*fun=*/ 0, nodes,
96                         /*extra=*/ 0, IGRAPH_EIGEN_LAPACK, &which,
97                         /*options=*/ 0, /*storage=*/ 0, &values2, &vectors2);
98 
99 #define DUMP() do {             \
100         igraph_vector_complex_print(&values);   \
101         igraph_vector_complex_print(&values2);  \
102     } while(0)
103 
104     for (i = 0; i < nodes; i++) {
105         int j;
106         igraph_real_t d =
107             igraph_complex_abs(igraph_complex_sub(VECTOR(values)[i],
108                                VECTOR(values2)[nodes - i - 1]));
109         if (d > 1e-15) {
110             DUMP();
111             return 2;
112         }
113         for (j = 0; j < nodes; j++) {
114             igraph_real_t d =
115                 igraph_complex_abs(igraph_complex_sub(MATRIX(vectors, j, i),
116                                    MATRIX(vectors2, j,
117                                           nodes - i - 1)));
118             if (d > 1e-15) {
119                 DUMP();
120                 return 3;
121             }
122         }
123     }
124 
125     igraph_vector_complex_destroy(&values);
126     igraph_matrix_complex_destroy(&vectors);
127     igraph_vector_complex_destroy(&values2);
128     igraph_matrix_complex_destroy(&vectors2);
129 
130     igraph_matrix_destroy(&mat2);
131 
132     VERIFY_FINALLY_STACK();
133 
134     return 0;
135 }
136