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