1 /* -*- mode: C -*- */
2 /* vim:set ts=4 sw=4 sts=4 noet: */
3 /*
4 IGraph library.
5 Copyright (C) 2007-2012 Gabor Csardi <csardi.gabor@gmail.com>
6 334 Harvard street, Cambridge, MA 02139 USA
7
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2 of the License, or
11 (at your option) any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 02110-1301 USA
22
23 */
24
25 #include "igraph_arpack.h"
26 #include "igraph_memory.h"
27
28 #include "linalg/arpack_internal.h"
29
30 #include <math.h>
31 #include <stdio.h>
32 #include <string.h>
33
34 /* The ARPACK example file dssimp.f is used as a template */
35
igraph_i_arpack_err_dsaupd(int error)36 static int igraph_i_arpack_err_dsaupd(int error) {
37 switch (error) {
38 case 1: return IGRAPH_ARPACK_MAXIT;
39 case 3: return IGRAPH_ARPACK_NOSHIFT;
40 case -1: return IGRAPH_ARPACK_NPOS;
41 case -2: return IGRAPH_ARPACK_NEVNPOS;
42 case -3: return IGRAPH_ARPACK_NCVSMALL;
43 case -4: return IGRAPH_ARPACK_NONPOSI;
44 case -5: return IGRAPH_ARPACK_WHICHINV;
45 case -6: return IGRAPH_ARPACK_BMATINV;
46 case -7: return IGRAPH_ARPACK_WORKLSMALL;
47 case -8: return IGRAPH_ARPACK_TRIDERR;
48 case -9: return IGRAPH_ARPACK_ZEROSTART;
49 case -10: return IGRAPH_ARPACK_MODEINV;
50 case -11: return IGRAPH_ARPACK_MODEBMAT;
51 case -12: return IGRAPH_ARPACK_ISHIFT;
52 case -13: return IGRAPH_ARPACK_NEVBE;
53 case -9999: return IGRAPH_ARPACK_NOFACT;
54 default: return IGRAPH_ARPACK_UNKNOWN;
55 }
56 }
57
igraph_i_arpack_err_dseupd(int error)58 static int igraph_i_arpack_err_dseupd(int error) {
59 switch (error) {
60 case -1: return IGRAPH_ARPACK_NPOS;
61 case -2: return IGRAPH_ARPACK_NEVNPOS;
62 case -3: return IGRAPH_ARPACK_NCVSMALL;
63 case -5: return IGRAPH_ARPACK_WHICHINV;
64 case -6: return IGRAPH_ARPACK_BMATINV;
65 case -7: return IGRAPH_ARPACK_WORKLSMALL;
66 case -8: return IGRAPH_ARPACK_TRIDERR;
67 case -9: return IGRAPH_ARPACK_ZEROSTART;
68 case -10: return IGRAPH_ARPACK_MODEINV;
69 case -11: return IGRAPH_ARPACK_MODEBMAT;
70 case -12: return IGRAPH_ARPACK_NEVBE;
71 case -14: return IGRAPH_ARPACK_FAILED;
72 case -15: return IGRAPH_ARPACK_HOWMNY;
73 case -16: return IGRAPH_ARPACK_HOWMNYS;
74 case -17: return IGRAPH_ARPACK_EVDIFF;
75 default: return IGRAPH_ARPACK_UNKNOWN;
76 }
77
78 }
79
igraph_i_arpack_err_dnaupd(int error)80 static int igraph_i_arpack_err_dnaupd(int error) {
81 switch (error) {
82 case 1: return IGRAPH_ARPACK_MAXIT;
83 case 3: return IGRAPH_ARPACK_NOSHIFT;
84 case -1: return IGRAPH_ARPACK_NPOS;
85 case -2: return IGRAPH_ARPACK_NEVNPOS;
86 case -3: return IGRAPH_ARPACK_NCVSMALL;
87 case -4: return IGRAPH_ARPACK_NONPOSI;
88 case -5: return IGRAPH_ARPACK_WHICHINV;
89 case -6: return IGRAPH_ARPACK_BMATINV;
90 case -7: return IGRAPH_ARPACK_WORKLSMALL;
91 case -8: return IGRAPH_ARPACK_TRIDERR;
92 case -9: return IGRAPH_ARPACK_ZEROSTART;
93 case -10: return IGRAPH_ARPACK_MODEINV;
94 case -11: return IGRAPH_ARPACK_MODEBMAT;
95 case -12: return IGRAPH_ARPACK_ISHIFT;
96 case -9999: return IGRAPH_ARPACK_NOFACT;
97 default: return IGRAPH_ARPACK_UNKNOWN;
98 }
99 }
100
igraph_i_arpack_err_dneupd(int error)101 static int igraph_i_arpack_err_dneupd(int error) {
102 switch (error) {
103 case 1: return IGRAPH_ARPACK_REORDER;
104 case -1: return IGRAPH_ARPACK_NPOS;
105 case -2: return IGRAPH_ARPACK_NEVNPOS;
106 case -3: return IGRAPH_ARPACK_NCVSMALL;
107 case -5: return IGRAPH_ARPACK_WHICHINV;
108 case -6: return IGRAPH_ARPACK_BMATINV;
109 case -7: return IGRAPH_ARPACK_WORKLSMALL;
110 case -8: return IGRAPH_ARPACK_SHUR;
111 case -9: return IGRAPH_ARPACK_LAPACK;
112 case -10: return IGRAPH_ARPACK_MODEINV;
113 case -11: return IGRAPH_ARPACK_MODEBMAT;
114 case -12: return IGRAPH_ARPACK_HOWMNYS;
115 case -13: return IGRAPH_ARPACK_HOWMNY;
116 case -14: return IGRAPH_ARPACK_FAILED;
117 case -15: return IGRAPH_ARPACK_EVDIFF;
118 default: return IGRAPH_ARPACK_UNKNOWN;
119 }
120 }
121
122 /**
123 * \function igraph_arpack_options_init
124 * Initialize ARPACK options
125 *
126 * Initializes ARPACK options, set them to default values.
127 * You can always pass the initialized \ref igraph_arpack_options_t
128 * object to built-in igraph functions without any modification. The
129 * built-in igraph functions modify the options to perform their
130 * calculation, e.g. \ref igraph_pagerank() always searches for the
131 * eigenvalue with the largest magnitude, regardless of the supplied
132 * value.
133 * </para><para>
134 * If you want to implement your own function involving eigenvalue
135 * calculation using ARPACK, however, you will likely need to set up
136 * the fields for yourself.
137 * \param o The \ref igraph_arpack_options_t object to initialize.
138 *
139 * Time complexity: O(1).
140 */
141
igraph_arpack_options_init(igraph_arpack_options_t * o)142 void igraph_arpack_options_init(igraph_arpack_options_t *o) {
143 o->bmat[0] = 'I';
144 o->n = 0; /* needs to be updated! */
145 o->which[0] = 'X'; o->which[1] = 'X';
146 o->nev = 1;
147 o->tol = 0;
148 o->ncv = 0; /* 0 means "automatic" */
149 o->ldv = o->n; /* will be updated to (real) n */
150 o->ishift = 1;
151 o->mxiter = 3000;
152 o->nb = 1;
153 o->mode = 1;
154 o->start = 0;
155 o->lworkl = 0;
156 o->sigma = 0;
157 o->sigmai = 0;
158 o->info = o->start;
159
160 o->iparam[0] = o->ishift; o->iparam[1] = 0; o->iparam[2] = o->mxiter; o->iparam[3] = o->nb;
161 o->iparam[4] = 0; o->iparam[5] = 0; o->iparam[6] = o->mode; o->iparam[7] = 0;
162 o->iparam[8] = 0; o->iparam[9] = 0; o->iparam[10] = 0;
163 }
164
165 /**
166 * \function igraph_arpack_storage_init
167 * Initialize ARPACK storage
168 *
169 * You only need this function if you want to run multiple eigenvalue
170 * calculations using ARPACK, and want to spare the memory
171 * allocation/deallocation between each two runs. Otherwise it is safe
172 * to supply a null pointer as the \c storage argument of both \ref
173 * igraph_arpack_rssolve() and \ref igraph_arpack_rnsolve() to make
174 * memory allocated and deallocated automatically.
175 *
176 * </para><para>Don't forget to call the \ref
177 * igraph_arpack_storage_destroy() function on the storage object if
178 * you don't need it any more.
179 * \param s The \ref igraph_arpack_storage_t object to initialize.
180 * \param maxn The maximum order of the matrices.
181 * \param maxncv The maximum NCV parameter intended to use.
182 * \param maxldv The maximum LDV parameter intended to use.
183 * \param symm Whether symmetric or non-symmetric problems will be
184 * solved using this \ref igraph_arpack_storage_t. (You cannot use
185 * the same storage both with symmetric and non-symmetric solvers.)
186 * \return Error code.
187 *
188 * Time complexity: O(maxncv*(maxldv+maxn)).
189 */
190
igraph_arpack_storage_init(igraph_arpack_storage_t * s,long int maxn,long int maxncv,long int maxldv,igraph_bool_t symm)191 int igraph_arpack_storage_init(igraph_arpack_storage_t *s, long int maxn,
192 long int maxncv, long int maxldv,
193 igraph_bool_t symm) {
194
195 /* TODO: check arguments */
196 s->maxn = (int) maxn;
197 s->maxncv = (int) maxncv;
198 s->maxldv = (int) maxldv;
199
200 #define CHECKMEM(x) \
201 if (!x) { \
202 IGRAPH_ERROR("Cannot allocate memory for ARPACK", IGRAPH_ENOMEM); \
203 } \
204 IGRAPH_FINALLY(igraph_free, x);
205
206 s->v = IGRAPH_CALLOC(maxldv * maxncv, igraph_real_t); CHECKMEM(s->v);
207 s->workd = IGRAPH_CALLOC(3 * maxn, igraph_real_t); CHECKMEM(s->workd);
208 s->d = IGRAPH_CALLOC(2 * maxncv, igraph_real_t); CHECKMEM(s->d);
209 s->resid = IGRAPH_CALLOC(maxn, igraph_real_t); CHECKMEM(s->resid);
210 s->ax = IGRAPH_CALLOC(maxn, igraph_real_t); CHECKMEM(s->ax);
211 s->select = IGRAPH_CALLOC(maxncv, int); CHECKMEM(s->select);
212
213 if (symm) {
214 s->workl = IGRAPH_CALLOC(maxncv * (maxncv + 8), igraph_real_t); CHECKMEM(s->workl);
215 s->di = 0;
216 s->workev = 0;
217 } else {
218 s->workl = IGRAPH_CALLOC(3 * maxncv * (maxncv + 2), igraph_real_t); CHECKMEM(s->workl);
219 s->di = IGRAPH_CALLOC(2 * maxncv, igraph_real_t); CHECKMEM(s->di);
220 s->workev = IGRAPH_CALLOC(3 * maxncv, igraph_real_t); CHECKMEM(s->workev);
221 IGRAPH_FINALLY_CLEAN(2);
222 }
223
224 #undef CHECKMEM
225
226 IGRAPH_FINALLY_CLEAN(7);
227 return 0;
228 }
229
230 /**
231 * \function igraph_arpack_storage_destroy
232 * Deallocate ARPACK storage
233 *
234 * \param s The \ref igraph_arpack_storage_t object for which the
235 * memory will be deallocated.
236 *
237 * Time complexity: operating system dependent.
238 */
239
igraph_arpack_storage_destroy(igraph_arpack_storage_t * s)240 void igraph_arpack_storage_destroy(igraph_arpack_storage_t *s) {
241
242 if (s->di) {
243 IGRAPH_FREE(s->di);
244 }
245 if (s->workev) {
246 IGRAPH_FREE(s->workev);
247 }
248
249 IGRAPH_FREE(s->workl);
250 IGRAPH_FREE(s->select);
251 IGRAPH_FREE(s->ax);
252 IGRAPH_FREE(s->resid);
253 IGRAPH_FREE(s->d);
254 IGRAPH_FREE(s->workd);
255 IGRAPH_FREE(s->v);
256 }
257
258 /**
259 * "Solver" for 1x1 eigenvalue problems since ARPACK sometimes blows up with
260 * these.
261 */
igraph_i_arpack_rssolve_1x1(igraph_arpack_function_t * fun,void * extra,igraph_arpack_options_t * options,igraph_vector_t * values,igraph_matrix_t * vectors)262 static int igraph_i_arpack_rssolve_1x1(igraph_arpack_function_t *fun, void *extra,
263 igraph_arpack_options_t* options,
264 igraph_vector_t* values, igraph_matrix_t* vectors) {
265 igraph_real_t a, b;
266 int nev = options->nev;
267
268 if (nev <= 0) {
269 IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_NEVNPOS);
270 }
271
272 /* Probe the value in the matrix */
273 a = 1;
274 if (fun(&b, &a, 1, extra)) {
275 IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
276 IGRAPH_ARPACK_PROD);
277 }
278
279 options->nconv = nev;
280
281 if (values != 0) {
282 IGRAPH_CHECK(igraph_vector_resize(values, 1));
283 VECTOR(*values)[0] = b;
284 }
285
286 if (vectors != 0) {
287 IGRAPH_CHECK(igraph_matrix_resize(vectors, 1, 1));
288 MATRIX(*vectors, 0, 0) = 1;
289 }
290
291 return IGRAPH_SUCCESS;
292 }
293
294 /**
295 * "Solver" for 1x1 eigenvalue problems since ARPACK sometimes blows up with
296 * these.
297 */
igraph_i_arpack_rnsolve_1x1(igraph_arpack_function_t * fun,void * extra,igraph_arpack_options_t * options,igraph_matrix_t * values,igraph_matrix_t * vectors)298 static int igraph_i_arpack_rnsolve_1x1(igraph_arpack_function_t *fun, void *extra,
299 igraph_arpack_options_t* options,
300 igraph_matrix_t* values, igraph_matrix_t* vectors) {
301 igraph_real_t a, b;
302 int nev = options->nev;
303
304 if (nev <= 0) {
305 IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_NEVNPOS);
306 }
307
308 /* Probe the value in the matrix */
309 a = 1;
310 if (fun(&b, &a, 1, extra)) {
311 IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
312 IGRAPH_ARPACK_PROD);
313 }
314
315 options->nconv = nev;
316
317 if (values != 0) {
318 IGRAPH_CHECK(igraph_matrix_resize(values, 1, 2));
319 MATRIX(*values, 0, 0) = b; MATRIX(*values, 0, 1) = 0;
320 }
321
322 if (vectors != 0) {
323 IGRAPH_CHECK(igraph_matrix_resize(vectors, 1, 1));
324 MATRIX(*vectors, 0, 0) = 1;
325 }
326
327 return IGRAPH_SUCCESS;
328 }
329
330 /**
331 * "Solver" for 2x2 nonsymmetric eigenvalue problems since ARPACK sometimes
332 * blows up with these.
333 */
igraph_i_arpack_rnsolve_2x2(igraph_arpack_function_t * fun,void * extra,igraph_arpack_options_t * options,igraph_matrix_t * values,igraph_matrix_t * vectors)334 static int igraph_i_arpack_rnsolve_2x2(igraph_arpack_function_t *fun, void *extra,
335 igraph_arpack_options_t* options, igraph_matrix_t* values,
336 igraph_matrix_t* vectors) {
337 igraph_real_t vec[2], mat[4];
338 igraph_real_t a, b, c, d;
339 igraph_real_t trace, det, tsq4_minus_d;
340 igraph_complex_t eval1, eval2;
341 igraph_complex_t evec1[2], evec2[2];
342 igraph_bool_t swap_evals = 0;
343 igraph_bool_t complex_evals = 0;
344 int nev = options->nev;
345
346 if (nev <= 0) {
347 IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_NEVNPOS);
348 }
349 if (nev > 2) {
350 nev = 2;
351 }
352
353 /* Probe the values in the matrix */
354 vec[0] = 1; vec[1] = 0;
355 if (fun(mat, vec, 2, extra)) {
356 IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
357 IGRAPH_ARPACK_PROD);
358 }
359 vec[0] = 0; vec[1] = 1;
360 if (fun(mat + 2, vec, 2, extra)) {
361 IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
362 IGRAPH_ARPACK_PROD);
363 }
364 a = mat[0]; b = mat[2]; c = mat[1]; d = mat[3];
365
366 /* Get the trace and the determinant */
367 trace = a + d;
368 det = a * d - b * c;
369 tsq4_minus_d = trace * trace / 4 - det;
370
371 /* Calculate the eigenvalues */
372 complex_evals = tsq4_minus_d < 0;
373 eval1 = igraph_complex_sqrt_real(tsq4_minus_d);
374 if (complex_evals) {
375 eval2 = igraph_complex_mul_real(eval1, -1);
376 } else {
377 /* to avoid having -0 in the imaginary part */
378 eval2 = igraph_complex(-IGRAPH_REAL(eval1), 0);
379 }
380 eval1 = igraph_complex_add_real(eval1, trace / 2);
381 eval2 = igraph_complex_add_real(eval2, trace / 2);
382
383 if (c != 0) {
384 evec1[0] = igraph_complex_sub_real(eval1, d);
385 evec1[1] = igraph_complex(c, 0);
386 evec2[0] = igraph_complex_sub_real(eval2, d);
387 evec2[1] = igraph_complex(c, 0);
388 } else if (b != 0) {
389 evec1[0] = igraph_complex(b, 0);
390 evec1[1] = igraph_complex_sub_real(eval1, a);
391 evec2[0] = igraph_complex(b, 0);
392 evec2[1] = igraph_complex_sub_real(eval2, a);
393 } else {
394 evec1[0] = igraph_complex(1, 0);
395 evec1[1] = igraph_complex(0, 0);
396 evec2[0] = igraph_complex(0, 0);
397 evec2[1] = igraph_complex(1, 0);
398 }
399
400 /* Sometimes we have to swap eval1 with eval2 and evec1 with eval2;
401 * determine whether we have to do it now */
402 if (options->which[0] == 'S') {
403 if (options->which[1] == 'M') {
404 /* eval1 must be the one with the smallest magnitude */
405 swap_evals = (igraph_complex_mod(eval1) > igraph_complex_mod(eval2));
406 } else if (options->which[1] == 'R') {
407 /* eval1 must be the one with the smallest real part */
408 swap_evals = (IGRAPH_REAL(eval1) > IGRAPH_REAL(eval2));
409 } else if (options->which[1] == 'I') {
410 /* eval1 must be the one with the smallest imaginary part */
411 swap_evals = (IGRAPH_IMAG(eval1) > IGRAPH_IMAG(eval2));
412 } else {
413 IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_WHICHINV);
414 }
415 } else if (options->which[0] == 'L') {
416 if (options->which[1] == 'M') {
417 /* eval1 must be the one with the largest magnitude */
418 swap_evals = (igraph_complex_mod(eval1) < igraph_complex_mod(eval2));
419 } else if (options->which[1] == 'R') {
420 /* eval1 must be the one with the largest real part */
421 swap_evals = (IGRAPH_REAL(eval1) < IGRAPH_REAL(eval2));
422 } else if (options->which[1] == 'I') {
423 /* eval1 must be the one with the largest imaginary part */
424 swap_evals = (IGRAPH_IMAG(eval1) < IGRAPH_IMAG(eval2));
425 } else {
426 IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_WHICHINV);
427 }
428 } else if (options->which[0] == 'X' && options->which[1] == 'X') {
429 /* No preference on the ordering of eigenvectors */
430 } else {
431 /* fprintf(stderr, "%c%c\n", options->which[0], options->which[1]); */
432 IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_WHICHINV);
433 }
434
435 options->nconv = nev;
436
437 if (swap_evals) {
438 igraph_complex_t dummy;
439 dummy = eval1; eval1 = eval2; eval2 = dummy;
440 dummy = evec1[0]; evec1[0] = evec2[0]; evec2[0] = dummy;
441 dummy = evec1[1]; evec1[1] = evec2[1]; evec2[1] = dummy;
442 }
443
444 if (complex_evals) {
445 /* The eigenvalues are conjugate pairs, so we store only the
446 * one with positive imaginary part */
447 if (IGRAPH_IMAG(eval1) < 0) {
448 eval1 = eval2;
449 evec1[0] = evec2[0]; evec1[1] = evec2[1];
450 }
451 }
452
453 if (values != 0) {
454 IGRAPH_CHECK(igraph_matrix_resize(values, nev, 2));
455 MATRIX(*values, 0, 0) = IGRAPH_REAL(eval1);
456 MATRIX(*values, 0, 1) = IGRAPH_IMAG(eval1);
457 if (nev > 1) {
458 MATRIX(*values, 1, 0) = IGRAPH_REAL(eval2);
459 MATRIX(*values, 1, 1) = IGRAPH_IMAG(eval2);
460 }
461 }
462
463 if (vectors != 0) {
464 if (complex_evals) {
465 IGRAPH_CHECK(igraph_matrix_resize(vectors, 2, 2));
466 MATRIX(*vectors, 0, 0) = IGRAPH_REAL(evec1[0]);
467 MATRIX(*vectors, 1, 0) = IGRAPH_REAL(evec1[1]);
468 MATRIX(*vectors, 0, 1) = IGRAPH_IMAG(evec1[0]);
469 MATRIX(*vectors, 1, 1) = IGRAPH_IMAG(evec1[1]);
470 } else {
471 IGRAPH_CHECK(igraph_matrix_resize(vectors, 2, nev));
472 MATRIX(*vectors, 0, 0) = IGRAPH_REAL(evec1[0]);
473 MATRIX(*vectors, 1, 0) = IGRAPH_REAL(evec1[1]);
474 if (nev > 1) {
475 MATRIX(*vectors, 0, 1) = IGRAPH_REAL(evec2[0]);
476 MATRIX(*vectors, 1, 1) = IGRAPH_REAL(evec2[1]);
477 }
478 }
479 }
480
481 return IGRAPH_SUCCESS;
482 }
483
484 /**
485 * "Solver" for symmetric 2x2 eigenvalue problems since ARPACK sometimes blows
486 * up with these.
487 */
igraph_i_arpack_rssolve_2x2(igraph_arpack_function_t * fun,void * extra,igraph_arpack_options_t * options,igraph_vector_t * values,igraph_matrix_t * vectors)488 static int igraph_i_arpack_rssolve_2x2(igraph_arpack_function_t *fun, void *extra,
489 igraph_arpack_options_t* options, igraph_vector_t* values,
490 igraph_matrix_t* vectors) {
491 igraph_real_t vec[2], mat[4];
492 igraph_real_t a, b, c, d;
493 igraph_real_t trace, det, tsq4_minus_d;
494 igraph_real_t eval1, eval2;
495 int nev = options->nev;
496
497 if (nev <= 0) {
498 IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_NEVNPOS);
499 }
500 if (nev > 2) {
501 nev = 2;
502 }
503
504 /* Probe the values in the matrix */
505 vec[0] = 1; vec[1] = 0;
506 if (fun(mat, vec, 2, extra)) {
507 IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
508 IGRAPH_ARPACK_PROD);
509 }
510 vec[0] = 0; vec[1] = 1;
511 if (fun(mat + 2, vec, 2, extra)) {
512 IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
513 IGRAPH_ARPACK_PROD);
514 }
515 a = mat[0]; b = mat[2]; c = mat[1]; d = mat[3];
516
517 /* Get the trace and the determinant */
518 trace = a + d;
519 det = a * d - b * c;
520 tsq4_minus_d = trace * trace / 4 - det;
521
522 if (tsq4_minus_d >= 0) {
523 /* Both eigenvalues are real */
524 eval1 = trace / 2 + sqrt(tsq4_minus_d);
525 eval2 = trace / 2 - sqrt(tsq4_minus_d);
526 if (c != 0) {
527 mat[0] = eval1 - d; mat[2] = eval2 - d;
528 mat[1] = c; mat[3] = c;
529 } else if (b != 0) {
530 mat[0] = b; mat[2] = b;
531 mat[1] = eval1 - a; mat[3] = eval2 - a;
532 } else {
533 mat[0] = 1; mat[2] = 0;
534 mat[1] = 0; mat[3] = 1;
535 }
536 } else {
537 /* Both eigenvalues are complex. Should not happen with symmetric
538 * matrices. */
539 IGRAPH_ERROR("ARPACK error, 2x2 matrix is not symmetric", IGRAPH_EINVAL);
540 }
541
542 /* eval1 is always the larger eigenvalue. If we want the smaller
543 * one, we have to swap eval1 with eval2 and also the columns of mat */
544 if (options->which[0] == 'S') {
545 trace = eval1; eval1 = eval2; eval2 = trace;
546 trace = mat[0]; mat[0] = mat[2]; mat[2] = trace;
547 trace = mat[1]; mat[1] = mat[3]; mat[3] = trace;
548 } else if (options->which[0] == 'L' || options->which[0] == 'B') {
549 /* Nothing to do here */
550 } else if (options->which[0] == 'X' && options->which[1] == 'X') {
551 /* No preference on the ordering of eigenvectors */
552 } else {
553 IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_WHICHINV);
554 }
555
556 options->nconv = nev;
557
558 if (values != 0) {
559 IGRAPH_CHECK(igraph_vector_resize(values, nev));
560 VECTOR(*values)[0] = eval1;
561 if (nev > 1) {
562 VECTOR(*values)[1] = eval2;
563 }
564 }
565
566 if (vectors != 0) {
567 IGRAPH_CHECK(igraph_matrix_resize(vectors, 2, nev));
568 MATRIX(*vectors, 0, 0) = mat[0];
569 MATRIX(*vectors, 1, 0) = mat[1];
570 if (nev > 1) {
571 MATRIX(*vectors, 0, 1) = mat[2];
572 MATRIX(*vectors, 1, 1) = mat[3];
573 }
574 }
575
576 return IGRAPH_SUCCESS;
577 }
578
igraph_arpack_rssort(igraph_vector_t * values,igraph_matrix_t * vectors,const igraph_arpack_options_t * options,igraph_real_t * d,const igraph_real_t * v)579 int igraph_arpack_rssort(igraph_vector_t *values, igraph_matrix_t *vectors,
580 const igraph_arpack_options_t *options,
581 igraph_real_t *d, const igraph_real_t *v) {
582
583 igraph_vector_t order;
584 char sort[2];
585 int apply = 1;
586 unsigned int n = (unsigned int) options->n;
587 int nconv = options->nconv;
588 int nev = options->nev;
589 unsigned int nans = (unsigned int) (nconv < nev ? nconv : nev);
590 unsigned int i;
591
592 #define which(a,b) (options->which[0]==a && options->which[1]==b)
593
594 if (which('L', 'A')) {
595 sort[0] = 'S'; sort[1] = 'A';
596 } else if (which('S', 'A')) {
597 sort[0] = 'L'; sort[1] = 'A';
598 } else if (which('L', 'M')) {
599 sort[0] = 'S'; sort[1] = 'M';
600 } else if (which('S', 'M')) {
601 sort[0] = 'L'; sort[1] = 'M';
602 } else if (which('B', 'E')) {
603 sort[0] = 'L'; sort[1] = 'A';
604 }
605
606 IGRAPH_CHECK(igraph_vector_init_seq(&order, 0, nconv - 1));
607 IGRAPH_FINALLY(igraph_vector_destroy, &order);
608 #ifdef HAVE_GFORTRAN
609 igraphdsortr_(sort, &apply, &nconv, d, VECTOR(order), /*which_len=*/ 2);
610 #else
611 igraphdsortr_(sort, &apply, &nconv, d, VECTOR(order));
612 #endif
613
614 /* BE is special */
615 if (which('B', 'E')) {
616 int w = 0, l1 = 0, l2 = nev - 1;
617 igraph_vector_t order2, d2;
618 IGRAPH_VECTOR_INIT_FINALLY(&order2, nev);
619 IGRAPH_VECTOR_INIT_FINALLY(&d2, nev);
620 while (l1 <= l2) {
621 VECTOR(order2)[w] = VECTOR(order)[l1];
622 VECTOR(d2)[w] = d[l1];
623 w++; l1++;
624 if (l1 <= l2) {
625 VECTOR(order2)[w] = VECTOR(order)[l2];
626 VECTOR(d2)[w] = d[l2];
627 w++; l2--;
628 }
629 }
630 igraph_vector_update(&order, &order2);
631 igraph_vector_copy_to(&d2, d);
632 igraph_vector_destroy(&order2);
633 igraph_vector_destroy(&d2);
634 IGRAPH_FINALLY_CLEAN(2);
635 }
636
637 #undef which
638
639 /* Copy values */
640 if (values) {
641 IGRAPH_CHECK(igraph_vector_resize(values, nans));
642 memcpy(VECTOR(*values), d, sizeof(igraph_real_t) * nans);
643 }
644
645 /* Reorder vectors */
646 if (vectors) {
647 IGRAPH_CHECK(igraph_matrix_resize(vectors, n, nans));
648 for (i = 0; i < nans; i++) {
649 unsigned int idx = (unsigned int) VECTOR(order)[i];
650 const igraph_real_t *ptr = v + n * idx;
651 memcpy(&MATRIX(*vectors, 0, i), ptr, sizeof(igraph_real_t) * n);
652 }
653 }
654
655 igraph_vector_destroy(&order);
656 IGRAPH_FINALLY_CLEAN(1);
657
658 return 0;
659 }
660
igraph_arpack_rnsort(igraph_matrix_t * values,igraph_matrix_t * vectors,const igraph_arpack_options_t * options,igraph_real_t * dr,igraph_real_t * di,igraph_real_t * v)661 int igraph_arpack_rnsort(igraph_matrix_t *values, igraph_matrix_t *vectors,
662 const igraph_arpack_options_t *options,
663 igraph_real_t *dr, igraph_real_t *di,
664 igraph_real_t *v) {
665
666 igraph_vector_t order;
667 char sort[2];
668 int apply = 1;
669 unsigned int n = (unsigned int) options->n;
670 int nconv = options->nconv;
671 int nev = options->nev;
672 unsigned int nans = (unsigned int) (nconv < nev ? nconv : nev);
673 unsigned int i;
674
675 #define which(a,b) (options->which[0]==a && options->which[1]==b)
676
677 if (which('L', 'M')) {
678 sort[0] = 'S'; sort[1] = 'M';
679 } else if (which('S', 'M')) {
680 sort[0] = 'L'; sort[1] = 'M';
681 } else if (which('L', 'R')) {
682 sort[0] = 'S'; sort[1] = 'R';
683 } else if (which('S', 'R')) {
684 sort[0] = 'L'; sort[1] = 'R';
685 } else if (which('L', 'I')) {
686 sort[0] = 'S'; sort[1] = 'I';
687 } else if (which('S', 'I')) {
688 sort[0] = 'L'; sort[1] = 'I';
689 }
690
691 #undef which
692
693 IGRAPH_CHECK(igraph_vector_init_seq(&order, 0, nconv - 1));
694 IGRAPH_FINALLY(igraph_vector_destroy, &order);
695 #ifdef HAVE_GFORTRAN
696 igraphdsortc_(sort, &apply, &nconv, dr, di, VECTOR(order), /*which_len=*/ 2);
697 #else
698 igraphdsortc_(sort, &apply, &nconv, dr, di, VECTOR(order));
699 #endif
700
701 if (values) {
702 IGRAPH_CHECK(igraph_matrix_resize(values, nans, 2));
703 memcpy(&MATRIX(*values, 0, 0), dr, sizeof(igraph_real_t) * nans);
704 memcpy(&MATRIX(*values, 0, 1), di, sizeof(igraph_real_t) * nans);
705 }
706
707 if (vectors) {
708 int nc = 0, nr = 0, ncol, vx = 0;
709 for (i = 0; i < nans; i++) {
710 if (di[i] == 0) {
711 nr++;
712 } else {
713 nc++;
714 }
715 }
716 ncol = (nc / 2) * 2 + (nc % 2) * 2 + nr;
717 IGRAPH_CHECK(igraph_matrix_resize(vectors, n, ncol));
718
719 for (i = 0; i < nans; i++) {
720 unsigned int idx;
721
722 idx = (unsigned int) VECTOR(order)[i];
723
724 if (di[i] == 0) {
725 /* real eigenvalue, single eigenvector */
726 memcpy(&MATRIX(*vectors, 0, vx), v + n * idx, sizeof(igraph_real_t) * n);
727 vx++;
728 } else if (di[i] > 0) {
729 /* complex eigenvalue, positive imaginary part encountered first.
730 * ARPACK stores its eigenvector directly in two consecutive columns.
731 * The complex conjugate pair of the eigenvalue (if any) will be in
732 * the next column and we will skip it because we advance 'i' below */
733 memcpy(&MATRIX(*vectors, 0, vx), v + n * idx, sizeof(igraph_real_t) * 2 * n);
734 vx += 2;
735 i++;
736 } else {
737 /* complex eigenvalue, negative imaginary part encountered first.
738 * The positive one will be the next one, but we need to copy the
739 * eigenvector corresponding to the eigenvalue with the positive
740 * imaginary part. */
741 idx = (unsigned int) VECTOR(order)[i + 1];
742 memcpy(&MATRIX(*vectors, 0, vx), v + n * idx, sizeof(igraph_real_t) * 2 * n);
743 vx += 2;
744 i++;
745 }
746 }
747 }
748
749 igraph_vector_destroy(&order);
750 IGRAPH_FINALLY_CLEAN(1);
751
752 if (values) {
753 /* Strive to include complex conjugate eigenvalue pairs in a way that the
754 * positive imaginary part comes first */
755 for (i = 0; i < nans; i++) {
756 if (MATRIX(*values, i, 1) == 0) {
757 /* Real eigenvalue, nothing to do */
758 } else if (MATRIX(*values, i, 1) < 0) {
759 /* Negative imaginary part came first; negate the imaginary part for
760 * this eigenvalue and the next one (which is the complex conjugate
761 * pair), and skip it */
762 MATRIX(*values, i, 1) *= -1;
763 i++;
764 if (i < nans) {
765 MATRIX(*values, i, 1) *= -1;
766 }
767 } else {
768 /* Positive imaginary part; skip the next eigenvalue, which is the
769 * complex conjugate pair */
770 i++;
771 }
772 }
773 }
774
775 return 0;
776 }
777
778 /**
779 * \function igraph_i_arpack_auto_ncv
780 * \brief Tries to set up the value of \c ncv in an \c igraph_arpack_options_t
781 * automagically.
782 */
igraph_i_arpack_auto_ncv(igraph_arpack_options_t * options)783 static void igraph_i_arpack_auto_ncv(igraph_arpack_options_t* options) {
784 /* This is similar to how Octave determines the value of ncv, with some
785 * modifications. */
786 int min_ncv = options->nev * 2 + 1;
787
788 /* Use twice the number of desired eigenvectors plus one by default */
789 options->ncv = min_ncv;
790 /* ...but use at least 20 Lanczos vectors... */
791 if (options->ncv < 20) {
792 options->ncv = 20;
793 }
794 /* ...but having ncv close to n leads to some problems with small graphs
795 * (example: PageRank of "A <--> C, D <--> E, B"), so we don't let it
796 * to be larger than n / 2...
797 */
798 if (options->ncv > options->n / 2) {
799 options->ncv = options->n / 2;
800 }
801 /* ...but we need at least min_ncv. */
802 if (options->ncv < min_ncv) {
803 options->ncv = min_ncv;
804 }
805 /* ...but at most n */
806 if (options->ncv > options->n) {
807 options->ncv = options->n;
808 }
809 }
810
811 /**
812 * \function igraph_i_arpack_report_no_convergence
813 * \brief Prints a warning that informs the user that the ARPACK solver
814 * did not converge.
815 */
igraph_i_arpack_report_no_convergence(const igraph_arpack_options_t * options)816 static void igraph_i_arpack_report_no_convergence(const igraph_arpack_options_t* options) {
817 char buf[1024];
818 snprintf(buf, sizeof(buf), "ARPACK solver failed to converge (%d iterations, "
819 "%d/%d eigenvectors converged)", options->iparam[2],
820 options->iparam[4], options->nev);
821 IGRAPH_WARNING(buf);
822 }
823
824 /**
825 * \function igraph_arpack_rssolve
826 * \brief ARPACK solver for symmetric matrices
827 *
828 * This is the ARPACK solver for symmetric matrices. Please use
829 * \ref igraph_arpack_rnsolve() for non-symmetric matrices.
830 * \param fun Pointer to an \ref igraph_arpack_function_t object,
831 * the function that performs the matrix-vector multiplication.
832 * \param extra An extra argument to be passed to \c fun.
833 * \param options An \ref igraph_arpack_options_t object.
834 * \param storage An \ref igraph_arpack_storage_t object, or a null
835 * pointer. In the latter case memory allocation and deallocation
836 * is performed automatically. Either this or the \p vectors argument
837 * must be non-null if the ARPACK iteration is started from a
838 * given starting vector. If both are given \p vectors take
839 * precedence.
840 * \param values If not a null pointer, then it should be a pointer to an
841 * initialized vector. The eigenvalues will be stored here. The
842 * vector will be resized as needed.
843 * \param vectors If not a null pointer, then it must be a pointer to
844 * an initialized matrix. The eigenvectors will be stored in the
845 * columns of the matrix. The matrix will be resized as needed.
846 * Either this or the \p vectors argument must be non-null if the
847 * ARPACK iteration is started from a given starting vector. If
848 * both are given \p vectors take precedence.
849 * \return Error code.
850 *
851 * Time complexity: depends on the matrix-vector
852 * multiplication. Usually a small number of iterations is enough, so
853 * if the matrix is sparse and the matrix-vector multiplication can be
854 * done in O(n) time (the number of vertices), then the eigenvalues
855 * are found in O(n) time as well.
856 */
857
igraph_arpack_rssolve(igraph_arpack_function_t * fun,void * extra,igraph_arpack_options_t * options,igraph_arpack_storage_t * storage,igraph_vector_t * values,igraph_matrix_t * vectors)858 int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra,
859 igraph_arpack_options_t *options,
860 igraph_arpack_storage_t *storage,
861 igraph_vector_t *values, igraph_matrix_t *vectors) {
862
863 igraph_real_t *v, *workl, *workd, *d, *resid, *ax;
864 igraph_bool_t free_them = 0;
865 int *select, i;
866
867 int ido = 0;
868 int rvec = vectors || storage ? 1 : 0; /* calculate eigenvectors? */
869 char *all = "All";
870
871 int origldv = options->ldv, origlworkl = options->lworkl,
872 orignev = options->nev, origncv = options->ncv;
873 igraph_real_t origtol = options->tol;
874 char origwhich[2];
875
876 origwhich[0] = options->which[0];
877 origwhich[1] = options->which[1];
878
879 /* Special case for 1x1 and 2x2 matrices in mode 1 */
880 if (options->mode == 1 && options->n == 1) {
881 return igraph_i_arpack_rssolve_1x1(fun, extra, options, values, vectors);
882 } else if (options->mode == 1 && options->n == 2) {
883 return igraph_i_arpack_rssolve_2x2(fun, extra, options, values, vectors);
884 }
885
886 /* Brush up options if needed */
887 if (options->ldv == 0) {
888 options->ldv = options->n;
889 }
890 if (options->ncv == 0) {
891 igraph_i_arpack_auto_ncv(options);
892 }
893 if (options->lworkl == 0) {
894 options->lworkl = options->ncv * (options->ncv + 8);
895 }
896 if (options->which[0] == 'X') {
897 options->which[0] = 'L';
898 options->which[1] = 'M';
899 }
900
901 if (storage) {
902 /* Storage provided */
903 if (storage->maxn < options->n) {
904 IGRAPH_ERROR("Not enough storage for ARPACK (`n')", IGRAPH_EINVAL);
905 }
906 if (storage->maxncv < options->ncv) {
907 IGRAPH_ERROR("Not enough storage for ARPACK (`ncv')", IGRAPH_EINVAL);
908 }
909 if (storage->maxldv < options->ldv) {
910 IGRAPH_ERROR("Not enough storage for ARPACK (`ldv')", IGRAPH_EINVAL);
911 }
912
913 v = storage->v;
914 workl = storage->workl;
915 workd = storage->workd;
916 d = storage->d;
917 resid = storage->resid;
918 ax = storage->ax;
919 select = storage->select;
920
921 } else {
922 /* Storage not provided */
923 free_them = 1;
924
925 #define CHECKMEM(x) \
926 if (!x) { \
927 IGRAPH_ERROR("Cannot allocate memory for ARPACK", IGRAPH_ENOMEM); \
928 } \
929 IGRAPH_FINALLY(igraph_free, x);
930
931 v = IGRAPH_CALLOC(options->ldv * options->ncv, igraph_real_t); CHECKMEM(v);
932 workl = IGRAPH_CALLOC(options->lworkl, igraph_real_t); CHECKMEM(workl);
933 workd = IGRAPH_CALLOC(3 * options->n, igraph_real_t); CHECKMEM(workd);
934 d = IGRAPH_CALLOC(2 * options->ncv, igraph_real_t); CHECKMEM(d);
935 resid = IGRAPH_CALLOC(options->n, igraph_real_t); CHECKMEM(resid);
936 ax = IGRAPH_CALLOC(options->n, igraph_real_t); CHECKMEM(ax);
937 select = IGRAPH_CALLOC(options->ncv, int); CHECKMEM(select);
938
939 #undef CHECKMEM
940
941 }
942
943 /* Set final bits */
944 options->bmat[0] = 'I';
945 options->iparam[0] = options->ishift;
946 options->iparam[1] = 0; // not referenced
947 options->iparam[2] = options->mxiter;
948 options->iparam[3] = 1; // currently dsaupd() works only for nb=1
949 options->iparam[4] = 0;
950 options->iparam[5] = 0; // not referenced
951 options->iparam[6] = options->mode;
952 options->iparam[7] = 0; // return value
953 options->iparam[8] = 0; // return value
954 options->iparam[9] = 0; // return value
955 options->iparam[10] = 0; // return value
956 options->info = options->start;
957 if (options->start) {
958 if (!storage && !vectors) {
959 IGRAPH_ERROR("Starting vector not given", IGRAPH_EINVAL);
960 }
961 if (vectors && (igraph_matrix_nrow(vectors) != options->n ||
962 igraph_matrix_ncol(vectors) != 1)) {
963 IGRAPH_ERROR("Invalid starting vector size", IGRAPH_EINVAL);
964 }
965 if (vectors) {
966 for (i = 0; i < options->n; i++) {
967 resid[i] = MATRIX(*vectors, i, 0);
968 }
969 }
970 }
971
972 /* Ok, we have everything */
973 while (1) {
974 #ifdef HAVE_GFORTRAN
975 igraphdsaupd_(&ido, options->bmat, &options->n, options->which,
976 &options->nev, &options->tol,
977 resid, &options->ncv, v, &options->ldv,
978 options->iparam, options->ipntr,
979 workd, workl, &options->lworkl, &options->info,
980 /*bmat_len=*/ 1, /*which_len=*/ 2);
981 #else
982 igraphdsaupd_(&ido, options->bmat, &options->n, options->which,
983 &options->nev, &options->tol,
984 resid, &options->ncv, v, &options->ldv,
985 options->iparam, options->ipntr,
986 workd, workl, &options->lworkl, &options->info);
987 #endif
988
989 if (ido == -1 || ido == 1) {
990 igraph_real_t *from = workd + options->ipntr[0] - 1;
991 igraph_real_t *to = workd + options->ipntr[1] - 1;
992 if (fun(to, from, options->n, extra) != 0) {
993 IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
994 IGRAPH_ARPACK_PROD);
995 }
996
997 } else {
998 break;
999 }
1000 }
1001
1002 if (options->info == 1) {
1003 igraph_i_arpack_report_no_convergence(options);
1004 }
1005 if (options->info != 0) {
1006 IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dsaupd(options->info));
1007 }
1008
1009 options->ierr = 0;
1010 #ifdef HAVE_GFORTRAN
1011 igraphdseupd_(&rvec, all, select, d, v, &options->ldv,
1012 &options->sigma, options->bmat, &options->n,
1013 options->which, &options->nev, &options->tol,
1014 resid, &options->ncv, v, &options->ldv, options->iparam,
1015 options->ipntr, workd, workl, &options->lworkl,
1016 &options->ierr, /*howmny_len=*/ 1, /*bmat_len=*/ 1,
1017 /*which_len=*/ 2);
1018 #else
1019 igraphdseupd_(&rvec, all, select, d, v, &options->ldv,
1020 &options->sigma, options->bmat, &options->n,
1021 options->which, &options->nev, &options->tol,
1022 resid, &options->ncv, v, &options->ldv, options->iparam,
1023 options->ipntr, workd, workl, &options->lworkl,
1024 &options->ierr);
1025 #endif
1026
1027 if (options->ierr != 0) {
1028 IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dseupd(options->ierr));
1029 }
1030
1031 /* Save the result */
1032
1033 options->noiter = options->iparam[2];
1034 options->nconv = options->iparam[4];
1035 options->numop = options->iparam[8];
1036 options->numopb = options->iparam[9];
1037 options->numreo = options->iparam[10];
1038
1039 if (options->nconv < options->nev) {
1040 IGRAPH_WARNING("Not enough eigenvalues/vectors in symmetric ARPACK "
1041 "solver");
1042 }
1043
1044 if (values || vectors) {
1045 IGRAPH_CHECK(igraph_arpack_rssort(values, vectors, options, d, v));
1046 }
1047
1048 options->ldv = origldv;
1049 options->ncv = origncv;
1050 options->lworkl = origlworkl;
1051 options->which[0] = origwhich[0]; options->which[1] = origwhich[1];
1052 options->tol = origtol;
1053 options->nev = orignev;
1054
1055 /* Clean up if needed */
1056 if (free_them) {
1057 IGRAPH_FREE(select);
1058 IGRAPH_FREE(ax);
1059 IGRAPH_FREE(resid);
1060 IGRAPH_FREE(d);
1061 IGRAPH_FREE(workd);
1062 IGRAPH_FREE(workl);
1063 IGRAPH_FREE(v);
1064 IGRAPH_FINALLY_CLEAN(7);
1065 }
1066 return 0;
1067 }
1068
1069 /**
1070 * \function igraph_arpack_rnsolve
1071 * \brief ARPACK solver for non-symmetric matrices
1072 *
1073 * Please always consider calling \ref igraph_arpack_rssolve() if your
1074 * matrix is symmetric, it is much faster.
1075 * \ref igraph_arpack_rnsolve() for non-symmetric matrices.
1076 * </para><para>
1077 * Note that ARPACK is not called for 2x2 matrices as an exact algebraic
1078 * solution exists in these cases.
1079 *
1080 * \param fun Pointer to an \ref igraph_arpack_function_t object,
1081 * the function that performs the matrix-vector multiplication.
1082 * \param extra An extra argument to be passed to \c fun.
1083 * \param options An \ref igraph_arpack_options_t object.
1084 * \param storage An \ref igraph_arpack_storage_t object, or a null
1085 * pointer. In the latter case memory allocation and deallocation
1086 * is performed automatically.
1087 * \param values If not a null pointer, then it should be a pointer to an
1088 * initialized matrix. The (possibly complex) eigenvalues will be
1089 * stored here. The matrix will have two columns, the first column
1090 * contains the real, the second the imaginary parts of the
1091 * eigenvalues.
1092 * The matrix will be resized as needed.
1093 * \param vectors If not a null pointer, then it must be a pointer to
1094 * an initialized matrix. The eigenvectors will be stored in the
1095 * columns of the matrix. The matrix will be resized as needed.
1096 * Note that real eigenvalues will have real eigenvectors in a single
1097 * column in this matrix; however, complex eigenvalues come in conjugate
1098 * pairs and the result matrix will store the eigenvector corresponding to
1099 * the eigenvalue with \em positive imaginary part only. Since in this case
1100 * the eigenvector is also complex, it will occupy \em two columns in the
1101 * eigenvector matrix (the real and the imaginary parts, in this order).
1102 * Caveat: if the eigenvalue vector returns only the eigenvalue with the
1103 * \em negative imaginary part for a complex conjugate eigenvalue pair, the
1104 * result vector will \em still store the eigenvector corresponding to the
1105 * eigenvalue with the positive imaginary part (since this is how ARPACK
1106 * works).
1107 * \return Error code.
1108 *
1109 * Time complexity: depends on the matrix-vector
1110 * multiplication. Usually a small number of iterations is enough, so
1111 * if the matrix is sparse and the matrix-vector multiplication can be
1112 * done in O(n) time (the number of vertices), then the eigenvalues
1113 * are found in O(n) time as well.
1114 */
1115
igraph_arpack_rnsolve(igraph_arpack_function_t * fun,void * extra,igraph_arpack_options_t * options,igraph_arpack_storage_t * storage,igraph_matrix_t * values,igraph_matrix_t * vectors)1116 int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra,
1117 igraph_arpack_options_t *options,
1118 igraph_arpack_storage_t *storage,
1119 igraph_matrix_t *values, igraph_matrix_t *vectors) {
1120
1121 igraph_real_t *v, *workl, *workd, *dr, *di, *resid, *workev;
1122 igraph_bool_t free_them = 0;
1123 int *select, i;
1124
1125 int ido = 0;
1126 int rvec = vectors || storage ? 1 : 0;
1127 char *all = "All";
1128
1129 int origldv = options->ldv, origlworkl = options->lworkl,
1130 orignev = options->nev, origncv = options->ncv;
1131 igraph_real_t origtol = options->tol;
1132 int d_size;
1133 char origwhich[2];
1134
1135 origwhich[0] = options->which[0];
1136 origwhich[1] = options->which[1];
1137
1138 /* Special case for 1x1 and 2x2 matrices in mode 1 */
1139 if (options->mode == 1 && options->n == 1) {
1140 return igraph_i_arpack_rnsolve_1x1(fun, extra, options, values, vectors);
1141 } else if (options->mode == 1 && options->n == 2) {
1142 return igraph_i_arpack_rnsolve_2x2(fun, extra, options, values, vectors);
1143 }
1144
1145 /* Brush up options if needed */
1146 if (options->ldv == 0) {
1147 options->ldv = options->n;
1148 }
1149 if (options->ncv == 0) {
1150 igraph_i_arpack_auto_ncv(options);
1151 }
1152 if (options->lworkl == 0) {
1153 options->lworkl = 3 * options->ncv * (options->ncv + 2);
1154 }
1155 if (options->which[0] == 'X') {
1156 options->which[0] = 'L';
1157 options->which[1] = 'M';
1158 }
1159
1160 if (storage) {
1161 /* Storage provided */
1162 if (storage->maxn < options->n) {
1163 IGRAPH_ERROR("Not enough storage for ARPACK (`n')", IGRAPH_EINVAL);
1164 }
1165 if (storage->maxncv < options->ncv) {
1166 IGRAPH_ERROR("Not enough storage for ARPACK (`ncv')", IGRAPH_EINVAL);
1167 }
1168 if (storage->maxldv < options->ldv) {
1169 IGRAPH_ERROR("Not enough storage for ARPACK (`ldv')", IGRAPH_EINVAL);
1170 }
1171
1172 v = storage->v;
1173 workl = storage->workl;
1174 workd = storage->workd;
1175 workev = storage->workev;
1176 dr = storage->d;
1177 di = storage->di;
1178 d_size = options->n;
1179 resid = storage->resid;
1180 select = storage->select;
1181
1182 } else {
1183 /* Storage not provided */
1184 free_them = 1;
1185
1186 #define CHECKMEM(x) \
1187 if (!x) { \
1188 IGRAPH_ERROR("Cannot allocate memory for ARPACK", IGRAPH_ENOMEM); \
1189 } \
1190 IGRAPH_FINALLY(igraph_free, x);
1191
1192 v = IGRAPH_CALLOC(options->n * options->ncv, igraph_real_t); CHECKMEM(v);
1193 workl = IGRAPH_CALLOC(options->lworkl, igraph_real_t); CHECKMEM(workl);
1194 workd = IGRAPH_CALLOC(3 * options->n, igraph_real_t); CHECKMEM(workd);
1195 d_size = 2 * options->nev + 1 > options->ncv ? 2 * options->nev + 1 : options->ncv;
1196 dr = IGRAPH_CALLOC(d_size, igraph_real_t); CHECKMEM(dr);
1197 di = IGRAPH_CALLOC(d_size, igraph_real_t); CHECKMEM(di);
1198 resid = IGRAPH_CALLOC(options->n, igraph_real_t); CHECKMEM(resid);
1199 select = IGRAPH_CALLOC(options->ncv, int); CHECKMEM(select);
1200 workev = IGRAPH_CALLOC(3 * options->ncv, igraph_real_t); CHECKMEM(workev);
1201
1202 #undef CHECKMEM
1203
1204 }
1205
1206 /* Set final bits */
1207 options->bmat[0] = 'I';
1208 options->iparam[0] = options->ishift;
1209 options->iparam[1] = 0; // not referenced
1210 options->iparam[2] = options->mxiter;
1211 options->iparam[3] = 1; // currently dnaupd() works only for nb=1
1212 options->iparam[4] = 0;
1213 options->iparam[5] = 0; // not referenced
1214 options->iparam[6] = options->mode;
1215 options->iparam[7] = 0; // return value
1216 options->iparam[8] = 0; // return value
1217 options->iparam[9] = 0; // return value
1218 options->iparam[10] = 0; // return value
1219 options->info = options->start;
1220 if (options->start) {
1221 if (!storage && !vectors) {
1222 IGRAPH_ERROR("Starting vector not given", IGRAPH_EINVAL);
1223 }
1224 if (vectors && (igraph_matrix_nrow(vectors) != options->n ||
1225 igraph_matrix_ncol(vectors) != 1)) {
1226 IGRAPH_ERROR("Invalid starting vector size", IGRAPH_EINVAL);
1227 }
1228 if (vectors) {
1229 for (i = 0; i < options->n; i++) {
1230 resid[i] = MATRIX(*vectors, i, 0);
1231 }
1232 }
1233 }
1234
1235 /* Ok, we have everything */
1236 while (1) {
1237 #ifdef HAVE_GFORTRAN
1238 igraphdnaupd_(&ido, options->bmat, &options->n, options->which,
1239 &options->nev, &options->tol,
1240 resid, &options->ncv, v, &options->ldv,
1241 options->iparam, options->ipntr,
1242 workd, workl, &options->lworkl, &options->info,
1243 /*bmat_len=*/ 1, /*which_len=*/ 2);
1244 #else
1245 igraphdnaupd_(&ido, options->bmat, &options->n, options->which,
1246 &options->nev, &options->tol,
1247 resid, &options->ncv, v, &options->ldv,
1248 options->iparam, options->ipntr,
1249 workd, workl, &options->lworkl, &options->info);
1250 #endif
1251
1252 if (ido == -1 || ido == 1) {
1253 igraph_real_t *from = workd + options->ipntr[0] - 1;
1254 igraph_real_t *to = workd + options->ipntr[1] - 1;
1255 if (fun(to, from, options->n, extra) != 0) {
1256 IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
1257 IGRAPH_ARPACK_PROD);
1258 }
1259 } else {
1260 break;
1261 }
1262 }
1263
1264 if (options->info == 1) {
1265 igraph_i_arpack_report_no_convergence(options);
1266 }
1267 if (options->info != 0 && options->info != -9999) {
1268 IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dnaupd(options->info));
1269 }
1270
1271 options->ierr = 0;
1272 #ifdef HAVE_GFORTRAN
1273 igraphdneupd_(&rvec, all, select, dr, di, v, &options->ldv,
1274 &options->sigma, &options->sigmai, workev, options->bmat,
1275 &options->n, options->which, &options->nev, &options->tol,
1276 resid, &options->ncv, v, &options->ldv, options->iparam,
1277 options->ipntr, workd, workl, &options->lworkl,
1278 &options->ierr, /*howmny_len=*/ 1, /*bmat_len=*/ 1,
1279 /*which_len=*/ 2);
1280 #else
1281 igraphdneupd_(&rvec, all, select, dr, di, v, &options->ldv,
1282 &options->sigma, &options->sigmai, workev, options->bmat,
1283 &options->n, options->which, &options->nev, &options->tol,
1284 resid, &options->ncv, v, &options->ldv, options->iparam,
1285 options->ipntr, workd, workl, &options->lworkl,
1286 &options->ierr);
1287 #endif
1288
1289 if (options->ierr != 0) {
1290 IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dneupd(options->info));
1291 }
1292
1293 /* Save the result */
1294
1295 options->noiter = options->iparam[2];
1296 options->nconv = options->iparam[4];
1297 options->numop = options->iparam[8];
1298 options->numopb = options->iparam[9];
1299 options->numreo = options->iparam[10];
1300
1301 if (options->nconv < options->nev) {
1302 IGRAPH_WARNING("Not enough eigenvalues/vectors in ARPACK "
1303 "solver");
1304 }
1305
1306 /* ARPACK might modify stuff in 'options' so reset everything that could
1307 * potentially get modified */
1308 options->ldv = origldv;
1309 options->ncv = origncv;
1310 options->lworkl = origlworkl;
1311 options->which[0] = origwhich[0]; options->which[1] = origwhich[1];
1312 options->tol = origtol;
1313 options->nev = orignev;
1314
1315 if (values || vectors) {
1316 IGRAPH_CHECK(igraph_arpack_rnsort(values, vectors, options,
1317 dr, di, v));
1318 }
1319
1320 /* Clean up if needed */
1321 if (free_them) {
1322 IGRAPH_FREE(workev);
1323 IGRAPH_FREE(select);
1324 IGRAPH_FREE(resid);
1325 IGRAPH_FREE(di);
1326 IGRAPH_FREE(dr);
1327 IGRAPH_FREE(workd);
1328 IGRAPH_FREE(workl);
1329 IGRAPH_FREE(v);
1330 IGRAPH_FINALLY_CLEAN(8);
1331 }
1332 return 0;
1333 }
1334
1335 /**
1336 * \function igraph_arpack_unpack_complex
1337 * \brief Make the result of the non-symmetric ARPACK solver more readable
1338 *
1339 * This function works on the output of \ref igraph_arpack_rnsolve and
1340 * brushes it up a bit: it only keeps \p nev eigenvalues/vectors and
1341 * every eigenvector is stored in two columns of the \p vectors
1342 * matrix.
1343 *
1344 * </para><para>
1345 * The output of the non-symmetric ARPACK solver is somewhat hard to
1346 * parse, as real eigenvectors occupy only one column in the matrix,
1347 * and the complex conjugate eigenvectors are not stored at all
1348 * (usually). The other problem is that the solver might return more
1349 * eigenvalues than requested. The common use of this function is to
1350 * call it directly after \ref igraph_arpack_rnsolve with its \p
1351 * vectors and \p values argument and \c options->nev as \p nev.
1352 * This will add the vectors for eigenvalues with a negative imaginary
1353 * part and return all vectors as 2 columns, a real and imaginary part.
1354 * \param vectors The eigenvector matrix, as returned by \ref
1355 * igraph_arpack_rnsolve. It will be resized, typically it will be
1356 * larger.
1357 * \param values The eigenvalue matrix, as returned by \ref
1358 * igraph_arpack_rnsolve. It will be resized, typically extra,
1359 * unneeded rows (=eigenvalues) will be removed.
1360 * \param nev The number of eigenvalues/vectors to keep. Can be less
1361 * or equal than the number originally requested from ARPACK.
1362 * \return Error code.
1363 *
1364 * Time complexity: linear in the number of elements in the \p vectors
1365 * matrix.
1366 */
1367
igraph_arpack_unpack_complex(igraph_matrix_t * vectors,igraph_matrix_t * values,long int nev)1368 int igraph_arpack_unpack_complex(igraph_matrix_t *vectors, igraph_matrix_t *values,
1369 long int nev) {
1370
1371 long int nodes = igraph_matrix_nrow(vectors);
1372 long int no_evs = igraph_matrix_nrow(values);
1373 long int i, j;
1374 long int new_vector_pos;
1375 long int vector_pos;
1376 igraph_matrix_t new_vectors;
1377
1378 /* Error checks */
1379 if (nev < 0) {
1380 IGRAPH_ERROR("`nev' cannot be negative", IGRAPH_EINVAL);
1381 }
1382 if (nev > no_evs) {
1383 IGRAPH_ERROR("`nev' too large, we don't have that many in `values'",
1384 IGRAPH_EINVAL);
1385 }
1386
1387 for (i = no_evs -1; i >= nev; i--) {
1388 IGRAPH_CHECK(igraph_matrix_remove_row(values, i));
1389 }
1390
1391 IGRAPH_CHECK(igraph_matrix_init(&new_vectors, nodes, nev * 2));
1392 IGRAPH_FINALLY(igraph_matrix_destroy, &new_vectors);
1393
1394 new_vector_pos = 0;
1395 vector_pos = 0;
1396 for (i = 0; i < nev && vector_pos < igraph_matrix_ncol(vectors); i++) {
1397 if (MATRIX(*values, i, 1) == 0) {
1398 /* Real eigenvalue */
1399 for (j = 0; j < nodes; j++) {
1400 MATRIX(new_vectors, j, new_vector_pos) = MATRIX(*vectors, j, vector_pos);
1401 }
1402 new_vector_pos += 2;
1403 vector_pos += 1;
1404 } else {
1405 /* complex eigenvalue */
1406 for (j = 0; j < nodes; j++) {
1407 MATRIX(new_vectors, j, new_vector_pos) = MATRIX(*vectors, j, vector_pos);
1408 MATRIX(new_vectors, j, new_vector_pos + 1) = MATRIX(*vectors, j, vector_pos + 1);
1409 }
1410
1411 /* handle the conjugate */
1412
1413 /* first check if the conjugate eigenvalue is there */
1414 i++;
1415 if (i >= nev) {
1416 break;
1417 }
1418
1419 if (MATRIX(*values, i, 1) != -MATRIX(*values, i-1, 1)) {
1420 IGRAPH_ERROR("Complex eigenvalue not followed by its conjugate.", IGRAPH_EINVAL);
1421 }
1422
1423 /* then copy and negate */
1424 for (j = 0; j < nodes; j++) {
1425 MATRIX(new_vectors, j, new_vector_pos + 2) = MATRIX(*vectors, j, vector_pos);
1426 MATRIX(new_vectors, j, new_vector_pos + 3) = -MATRIX(*vectors, j, vector_pos + 1);
1427 }
1428 new_vector_pos += 4;
1429 vector_pos += 2;
1430 }
1431 }
1432 igraph_matrix_destroy(vectors);
1433 IGRAPH_CHECK(igraph_matrix_copy(vectors, &new_vectors));
1434 igraph_matrix_destroy(&new_vectors);
1435 IGRAPH_FINALLY_CLEAN(1);
1436
1437 return IGRAPH_SUCCESS;
1438 }
1439