1 /* -*- mode: C -*- */
2 /*
3 IGraph library.
4 Copyright (C) 2011-2012 Gabor Csardi <csardi.gabor@gmail.com>
5 334 Harvard street, 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 #include <pthread.h>
26 #include <string.h>
27
28 #include "linalg/arpack_internal.h"
29
30 #include "test_utilities.inc"
31
32 /* Test whether ARPACK is thread-safe. We will create two threads,
33 each calling a different ARPACK eigensolver. We will make sure that
34 the ARPACK calls from the two threads overlap */
35
36 typedef struct thread_data_t {
37 igraph_matrix_t *m;
38 igraph_vector_t *result;
39 pthread_cond_t *cond;
40 pthread_mutex_t *mutex;
41 int *steps, *othersteps;
42 } thread_data_t;
43
arpack_mult(igraph_real_t * to,igraph_real_t * from,int n,igraph_matrix_t * matrix)44 int arpack_mult(igraph_real_t *to, igraph_real_t *from, int n,
45 igraph_matrix_t *matrix) {
46 /* TODO */
47 igraph_blas_dgemv_array(/*transpose=*/ 0, /*alpha=*/ 1.0, matrix,
48 from, /*beta=*/ 0.0, to);
49
50 return 0;
51 }
52
53 /* This is the function performed by each thread. It calls the
54 low-level ARPACK symmetric eigensolver, step by step. After each
55 step, it synchronizes with the other thread.
56
57 The synchronization ensures that the two threads are using the
58 thread-local variables at the same time. If they are really
59 thread-local, then ARPACK still delivers the correct solution for
60 the two matrices. Otherwise the result is undefined: maybe results
61 will be incorrect, or the program will crash.
62
63 This function is basically a simplified copy of igraph_arpack_rssolve.
64 */
65
thread_function(void * arg)66 void *thread_function(void *arg) {
67 thread_data_t *data = (thread_data_t*) arg;
68 igraph_matrix_t *M = data->m;
69 igraph_vector_t *result = data->result;
70 pthread_cond_t *cond = data->cond;
71 pthread_mutex_t *mutex = data->mutex;
72 igraph_arpack_options_t options;
73 igraph_real_t *v, *workl, *workd, *d, *resid, *ax;
74 int *select;
75 int ido = 0;
76 #if IGRAPH_THREAD_SAFE
77 int rvec = 1;
78 char *all = "All";
79 #endif
80 int i;
81
82 igraph_arpack_options_init(&options);
83 options.n = igraph_matrix_nrow(M);
84 options.ldv = options.n;
85 options.nev = 1;
86 options.ncv = 3;
87 options.lworkl = options.ncv * (options.ncv + 8);
88 options.which[0] = 'L';
89 options.which[1] = 'M';
90
91 options.iparam[0] = options.ishift;
92 options.iparam[2] = options.mxiter;
93 options.iparam[3] = options.nb;
94 options.iparam[4] = 0;
95 options.iparam[6] = options.mode;
96 options.info = options.start;
97
98 v = IGRAPH_CALLOC(options.ldv * options.ncv, igraph_real_t);
99 workl = IGRAPH_CALLOC(options.lworkl, igraph_real_t);
100 workd = IGRAPH_CALLOC(3 * options.n, igraph_real_t);
101 d = IGRAPH_CALLOC(2 * options.ncv, igraph_real_t);
102 resid = IGRAPH_CALLOC(options.n, igraph_real_t);
103 ax = IGRAPH_CALLOC(options.n, igraph_real_t);
104 select = IGRAPH_CALLOC(options.ncv, int);
105
106 if (!v || !workl || !workd || !d || !resid || !ax || !select) {
107 printf("Out of memory\n");
108 return 0;
109 }
110
111 while (1) {
112 #if IGRAPH_THREAD_SAFE
113 igraphdsaupd_(&ido, options.bmat, &options.n, options.which,
114 &options.nev, &options.tol, resid, &options.ncv, v,
115 &options.ldv, options.iparam, options.ipntr, workd,
116 workl, &options.lworkl, &options.info);
117 #endif
118
119 if (ido == -1 || ido == 1) {
120
121 igraph_real_t *from = workd + options.ipntr[0] - 1;
122 igraph_real_t *to = workd + options.ipntr[1] - 1;
123 arpack_mult(to, from, options.n, M);
124
125 } else {
126 break;
127 }
128
129 pthread_mutex_lock(mutex);
130 *(data->steps) += 1;
131 if ( *(data->othersteps) == *(data->steps) ) {
132 pthread_cond_signal(cond);
133 }
134
135 while ( *(data->othersteps) < * (data->steps) && *(data->othersteps) != -1 ) {
136 pthread_cond_wait(cond, mutex);
137 }
138 pthread_mutex_unlock(mutex);
139 }
140
141 pthread_mutex_lock(mutex);
142 *data->steps = -1;
143 pthread_cond_signal(cond);
144 pthread_mutex_unlock(mutex);
145
146 if (options.info != 0) {
147 printf("ARPACK error\n");
148 return 0;
149 }
150
151 #if IGRAPH_THREAD_SAFE
152 igraphdseupd_(&rvec, all, select, d, v, &options.ldv,
153 &options.sigma, options.bmat, &options.n,
154 options.which, &options.nev, &options.tol,
155 resid, &options.ncv, v, &options.ldv, options.iparam,
156 options.ipntr, workd, workl, &options.lworkl,
157 &options.ierr);
158 #endif
159
160 if (options.ierr != 0) {
161 printf("ARPACK error\n");
162 return 0;
163 }
164
165 igraph_vector_resize(result, options.n);
166 for (i = 0; i < options.n; i++) {
167 VECTOR(*result)[i] = v[i];
168 }
169
170 free(v);
171 free(workl);
172 free(workd);
173 free(d);
174 free(resid);
175 free(ax);
176 free(select);
177
178 return 0;
179 }
180
main()181 int main() {
182 pthread_t thread_id1, thread_id2;
183 void *exit_status1, *exit_status2;
184 igraph_matrix_t m1, m2;
185 igraph_vector_t result1, result2;
186 pthread_cond_t steps_cond = PTHREAD_COND_INITIALIZER;
187 pthread_mutex_t steps_mutex = PTHREAD_MUTEX_INITIALIZER;
188 int steps1 = 0, steps2 = 0;
189 thread_data_t
190 data1 = { &m1, &result1, &steps_cond, &steps_mutex, &steps1, &steps2 },
191 data2 = { &m2, &result2, &steps_cond, &steps_mutex, &steps2, &steps1 };
192 int i, j;
193
194 /* Skip if igraph is not thread safe */
195 if (!IGRAPH_THREAD_SAFE) {
196 return 77;
197 }
198
199 igraph_matrix_init(&m1, 10, 10);
200 igraph_matrix_init(&m2, 10, 10);
201 igraph_vector_init(&result1, igraph_matrix_nrow(&m1));
202 igraph_vector_init(&result2, igraph_matrix_nrow(&m2));
203
204 igraph_rng_seed(igraph_rng_default(), 42);
205
206 for (i = 0; i < igraph_matrix_nrow(&m1); i++) {
207 for (j = 0; j <= i; j++) {
208 MATRIX(m1, i, j) = MATRIX(m1, j, i) =
209 igraph_rng_get_integer(igraph_rng_default(), 0, 10);
210 }
211 }
212
213 for (i = 0; i < igraph_matrix_nrow(&m2); i++) {
214 for (j = 0; j <= i; j++) {
215 MATRIX(m2, i, j) = MATRIX(m2, j, i) =
216 igraph_rng_get_integer(igraph_rng_default(), 0, 10);
217 }
218 }
219
220 pthread_create(&thread_id1, NULL, thread_function, (void *) &data1);
221 pthread_create(&thread_id2, NULL, thread_function, (void *) &data2);
222
223 pthread_join(thread_id1, &exit_status1);
224 pthread_join(thread_id2, &exit_status2);
225
226 igraph_matrix_print(&m1);
227 igraph_vector_print(&result1);
228 printf("---\n");
229 igraph_matrix_print(&m2);
230 igraph_vector_print(&result2);
231
232 igraph_vector_destroy(&result1);
233 igraph_vector_destroy(&result2);
234 igraph_matrix_destroy(&m1);
235 igraph_matrix_destroy(&m2);
236
237 VERIFY_FINALLY_STACK();
238
239 return 0;
240 }
241