From b600e85a82a5d60cf1c9d3f03d4ae129ae39073a Mon Sep 17 00:00:00 2001 From: Tamas Nepusz Date: Wed, 1 Aug 2018 19:18:20 +0200 Subject: [PATCH] [PATCH] fixing issues with ARPACK 3.6, related to #1107 Remark: This does not solve the build issue Gbp-Pq: Name fix_test_arpack-3.6.patch --- examples/simple/igraph_arpack_rnsolve.c | 153 ++++++++++---- examples/simple/igraph_arpack_rnsolve.out | 72 ++----- src/arpack.c | 230 ++++++++++++++-------- 3 files changed, 270 insertions(+), 185 deletions(-) diff --git a/examples/simple/igraph_arpack_rnsolve.c b/examples/simple/igraph_arpack_rnsolve.c index 5ca2c69..3141426 100644 --- a/examples/simple/igraph_arpack_rnsolve.c +++ b/examples/simple/igraph_arpack_rnsolve.c @@ -1,22 +1,22 @@ /* -*- mode: C -*- */ -/* +/* IGraph library. Copyright (C) 2012 Gabor Csardi 334 Harvard st, Cambridge MA, 02139 USA - + This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. - + This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - + You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software - Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ @@ -29,18 +29,112 @@ typedef struct cb2_data_t { int cb2(igraph_real_t *to, const igraph_real_t *from, int n, void *extra) { cb2_data_t *data=(cb2_data_t*) extra; - igraph_blas_dgemv_array(/*transpose=*/ 0, /*alpha=*/ 1.0, + igraph_blas_dgemv_array(/*transpose=*/ 0, /*alpha=*/ 1.0, data->A, from, /*beta=*/ 0.0, to); return 0; } +int check_eigenvector( + const char* test_name, + igraph_matrix_t* A, igraph_matrix_t* values, igraph_matrix_t* vectors, + int eval_idx, int evec_col_idx +) { + igraph_complex_t eval, prod; + igraph_complex_t *evec; + int i, j, n = igraph_matrix_nrow(A); + + eval = igraph_complex(MATRIX(*values, eval_idx, 0), MATRIX(*values, eval_idx, 1)); + evec = (igraph_complex_t*) calloc(n, sizeof(igraph_complex_t)); + if (IGRAPH_IMAG(eval) == 0) { + /* Real eigenvalue, so we have a real eigenvector */ + for (i = 0; i < n; i++) { + evec[i] = igraph_complex(MATRIX(*vectors, i, evec_col_idx), 0); + } + } else { + /* Complex eigenvalue pair, so we have a complex eigenvector pair */ + /* ARPACK always stores the eigenvector corresponding to the eigenvalue + * with a positive imaginary part. If the imaginary part is negative, we + * need to multiply the imaginary part of the eigenvector by -1 */ + for (i = 0; i < n; i++) { + evec[i] = igraph_complex( + MATRIX(*vectors, i, evec_col_idx), + MATRIX(*vectors, i, evec_col_idx+1) * ( + IGRAPH_IMAG(eval) < 0 ? -1 : 1 + ) + ); + } + } + + /* Multiply matrix with eigenvector */ + for (i = 0; i < n; i++) { + prod = igraph_complex(0, 0); + for (j = 0; j < n; j++) { + prod = igraph_complex_add( + igraph_complex_mul_real(evec[j], MATRIX(*A, i, j)), + prod + ); + } + prod = igraph_complex_div(prod, eval); + if (!igraph_complex_eq_tol(prod, evec[i], 1e-6)) { + prod = igraph_complex_sub(prod, evec[i]); + printf("%s: vector corresponding to eigenvalue (%.4f + %.4f*i) is not an " + "eigenvector, coordinate %d differs by %.4f + %.4f*i\n", + test_name, IGRAPH_REAL(eval), IGRAPH_IMAG(eval), + i, IGRAPH_REAL(prod), IGRAPH_IMAG(prod)); + return 1; + } + } + + /* Free stuff */ + free(evec); + + return 0; +} + +int check_eigenvectors( + const char* test_name, + igraph_matrix_t* A, igraph_matrix_t* values, igraph_matrix_t* vectors +) { + int i, j; + int nev = igraph_matrix_nrow(values); + int errors = 0; + igraph_bool_t conjugate_pair_will_come = 0; + + for (i = 0, j = 0; i < nev; i++) { + errors += check_eigenvector(test_name, A, values, vectors, i, j); + if (MATRIX(*values, i, 1) != 0) { + /* Complex eigenvalue */ + if (conjugate_pair_will_come) { + j += 2; + conjugate_pair_will_come = 0; + } else { + conjugate_pair_will_come = 1; + } + } else { + /* Real eigenvalue */ + j++; + } + } + return (errors > 0) ? 1 : 0; +} + +void print_debug_output( + igraph_matrix_t* values, igraph_matrix_t* vectors +) { + printf("---\n"); + igraph_matrix_print(values); + printf("---\n"); + igraph_matrix_print(vectors); + printf("---\n"); +} + #define DIM 10 int main() { igraph_matrix_t A; igraph_matrix_t values, vectors; igraph_arpack_options_t options; - cb2_data_t data = { &A }; + cb2_data_t data = { &A }; int i, j; igraph_rng_seed(igraph_rng_default(), 42 * 42); @@ -52,7 +146,10 @@ int main() { MATRIX(A, i, j) = igraph_rng_get_integer(igraph_rng_default(), -10, 10); } } - + + igraph_matrix_print(&A); + printf("===\n"); + igraph_arpack_options_init(&options); options.n=DIM; options.start=0; @@ -66,13 +163,8 @@ int main() { igraph_arpack_rnsolve(cb2, /*extra=*/ &data, &options, /*storage=*/ 0, &values, &vectors); - while (options.nev < igraph_matrix_nrow(&values)) { - igraph_matrix_remove_row(&values, igraph_matrix_nrow(&values)-1); - } - - if (MATRIX(values, 2, 1) > 0) { - MATRIX(values, 2, 1) = -MATRIX(values, 2, 1); - MATRIX(values, 3, 1) = -MATRIX(values, 3, 1); + if (check_eigenvectors("LM #1", &A, &values, &vectors)) { + print_debug_output(&values, &vectors); } igraph_matrix_print(&values); @@ -88,19 +180,10 @@ int main() { igraph_arpack_rnsolve(cb2, /*extra=*/ &data, &options, /*storage=*/ 0, &values, &vectors); - while (options.nev < igraph_matrix_nrow(&values)) { - igraph_matrix_remove_row(&values, igraph_matrix_nrow(&values)-1); - } - - if (MATRIX(values, 2, 1) > 0) { - MATRIX(values, 2, 1) = -MATRIX(values, 2, 1); + if (check_eigenvectors("LM #1", &A, &values, &vectors)) { + print_debug_output(&values, &vectors); } - igraph_matrix_print(&values); - printf("---\n"); - igraph_matrix_print(&vectors); - printf("---\n"); - /* -------------- */ options.nev=3; @@ -109,14 +192,9 @@ int main() { igraph_arpack_rnsolve(cb2, /*extra=*/ &data, &options, /*storage=*/ 0, &values, &vectors); - while (options.nev < igraph_matrix_nrow(&values)) { - igraph_matrix_remove_row(&values, igraph_matrix_nrow(&values)-1); + if (check_eigenvectors("SR", &A, &values, &vectors)) { + print_debug_output(&values, &vectors); } - - igraph_matrix_print(&values); - printf("---\n"); - igraph_matrix_print(&vectors); - printf("---\n"); /* -------------- */ @@ -126,14 +204,9 @@ int main() { igraph_arpack_rnsolve(cb2, /*extra=*/ &data, &options, /*storage=*/ 0, &values, &vectors); - while (options.nev < igraph_matrix_nrow(&values)) { - igraph_matrix_remove_row(&values, igraph_matrix_nrow(&values)-1); + if (check_eigenvectors("LI", &A, &values, &vectors)) { + print_debug_output(&values, &vectors); } - - igraph_matrix_print(&values); - printf("---\n"); - igraph_matrix_print(&vectors); - printf("---\n"); /* -------------- */ diff --git a/examples/simple/igraph_arpack_rnsolve.out b/examples/simple/igraph_arpack_rnsolve.out index ec38497..93e0a44 100644 --- a/examples/simple/igraph_arpack_rnsolve.out +++ b/examples/simple/igraph_arpack_rnsolve.out @@ -1,61 +1,11 @@ -22.0483 0 --21.3281 0 --3.00735 -19.2957 --3.00735 19.2957 ---- --0.373224 0.226696 0.0473383 -0.204213 -0.289145 -0.296079 0.156365 0.0479785 --0.539167 -0.046114 0.189081 0.152399 --0.155332 0.358022 -0.364566 0.0248164 --0.217492 -0.618375 -0.30221 -0.084056 -0.395358 0.276259 -0.0368916 -0.429673 --0.323702 0.196714 0.132731 -0.336104 --0.177548 0.177994 0.382879 -0.0409697 -0.348174 0.0512606 -0.305212 0.141094 -0.0335622 0.446014 0.239725 0.0551197 ---- -22.0483 0 --21.3281 0 --3.00735 -19.2957 ---- -0.373224 0.226696 0.204213 0.0473383 --0.289145 -0.296079 -0.0479785 0.156365 -0.539167 -0.046114 -0.152399 0.189081 -0.155332 0.358022 -0.0248164 -0.364566 -0.217492 -0.618375 0.084056 -0.30221 --0.395358 0.276259 0.429673 -0.0368916 -0.323702 0.196714 0.336104 0.132731 -0.177548 0.177994 0.0409697 0.382879 --0.348174 0.0512606 -0.141094 -0.305212 --0.0335622 0.446014 -0.0551197 0.239725 ---- --21.3281 0 --12.4527 0 --3.00735 -19.2957 ---- --0.226696 0.695866 -0.204213 -0.0473383 -0.296079 0.120213 0.0479785 -0.156365 -0.046114 0.0628274 0.152399 -0.189081 --0.358022 -0.0817567 0.0248164 0.364566 -0.618375 -0.354011 -0.084056 0.30221 --0.276259 0.33649 -0.429673 0.0368916 --0.196714 0.284155 -0.336104 -0.132731 --0.177994 0.0164834 -0.0409697 -0.382879 --0.0512606 0.175732 0.141094 0.305212 --0.446014 0.374487 0.0551197 -0.239725 ---- --3.00735 19.2957 --3.00735 -19.2957 -12.1099 6.27293 ---- -0.0768616 -0.195028 -0.152389 0.21912 -0.147607 0.0704569 0.346547 0.125122 -0.164607 0.178554 -0.153007 -0.188931 --0.364251 -0.0290787 0.0368355 0.21496 --0.286559 -0.127595 -0.277981 0.264759 -0.0267116 -0.430426 0.246402 0.129704 -0.180726 -0.312925 -0.262001 -0.40777 -0.384741 0.0157948 -0.296988 0.134903 --0.322646 0.0946649 -0.148076 0.121394 -0.229009 0.089782 -0.271714 0.0980886 ---- +-6 0 10 3 8 1 -4 10 -8 0 +-6 1 0 8 -4 4 -7 1 1 6 +7 -7 8 6 -4 -8 -1 -7 -3 -7 +6 8 -4 -1 10 3 7 7 -3 -8 +1 -7 -4 9 0 5 5 6 -8 10 +-9 10 -5 -9 5 3 -5 7 -7 10 +-3 0 8 -6 -2 -7 1 -3 -8 1 +2 0 9 -3 0 -9 -4 0 10 0 +-9 1 -6 -1 7 10 9 9 8 -2 +-7 1 9 -7 10 -1 -2 -5 7 6 +=== diff --git a/src/arpack.c b/src/arpack.c index 500db65..24de03f 100644 --- a/src/arpack.c +++ b/src/arpack.c @@ -1,23 +1,23 @@ /* -*- mode: C -*- */ -/* vim:set sw=2 ts=2 sts=2 et: */ -/* +/* vim:set sw=2 ts=8 sts=2 noet: */ +/* IGraph library. Copyright (C) 2007-2012 Gabor Csardi 334 Harvard street, Cambridge, MA 02139 USA - + This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. - + This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - + You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software - Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ @@ -73,7 +73,7 @@ int igraph_i_arpack_err_dseupd(int error) { case -17: return IGRAPH_ARPACK_EVDIFF; default: return IGRAPH_ARPACK_UNKNOWN; } - + } int igraph_i_arpack_err_dnaupd(int error) { @@ -155,7 +155,7 @@ void igraph_arpack_options_init(igraph_arpack_options_t *o) { o->sigma=0; o->sigmai=0; o->info=o->start; - + o->iparam[0]=o->ishift; o->iparam[1]=0; o->iparam[2]=o->mxiter; o->iparam[3]=o->nb; o->iparam[4]=0; o->iparam[5]=0; o->iparam[6]=o->mode; o->iparam[7]=0; o->iparam[8]=0; o->iparam[9]=0; o->iparam[10]=0; @@ -190,7 +190,7 @@ void igraph_arpack_options_init(igraph_arpack_options_t *o) { int igraph_arpack_storage_init(igraph_arpack_storage_t *s, long int maxn, long int maxncv, long int maxldv, igraph_bool_t symm) { - + /* TODO: check arguments */ s->maxn=(int) maxn; s->maxncv=(int) maxncv; @@ -221,7 +221,7 @@ int igraph_arpack_storage_init(igraph_arpack_storage_t *s, long int maxn, } #undef CHECKMEM - + IGRAPH_FINALLY_CLEAN(7); return 0; } @@ -237,14 +237,14 @@ int igraph_arpack_storage_init(igraph_arpack_storage_t *s, long int maxn, */ void igraph_arpack_storage_destroy(igraph_arpack_storage_t *s) { - + if (s->di) { igraph_Free(s->di); } if (s->workev) { igraph_Free(s->workev); } - + igraph_Free(s->workl); igraph_Free(s->select); igraph_Free(s->ax); @@ -320,7 +320,7 @@ int igraph_i_arpack_rnsolve_1x1(igraph_arpack_function_t *fun, void *extra, if (vectors != 0) { IGRAPH_CHECK(igraph_matrix_resize(vectors, 1, 1)); - MATRIX(*vectors, 0, 0) = 1; + MATRIX(*vectors, 0, 0) = 1; } return IGRAPH_SUCCESS; @@ -574,7 +574,7 @@ int igraph_i_arpack_rssolve_2x2(igraph_arpack_function_t *fun, void *extra, } int igraph_arpack_rssort(igraph_vector_t *values, igraph_matrix_t *vectors, - const igraph_arpack_options_t *options, + const igraph_arpack_options_t *options, igraph_real_t *d, const igraph_real_t *v) { igraph_vector_t order; @@ -614,12 +614,12 @@ int igraph_arpack_rssort(igraph_vector_t *values, igraph_matrix_t *vectors, IGRAPH_VECTOR_INIT_FINALLY(&order2, nev); IGRAPH_VECTOR_INIT_FINALLY(&d2, nev); while (l1 <= l2) { - VECTOR(order2)[w] = VECTOR(order)[l1]; - VECTOR(d2)[w]=d[l1]; + VECTOR(order2)[w] = VECTOR(order)[l1]; + VECTOR(d2)[w]=d[l1]; w++; l1++; if (l1 <= l2) { - VECTOR(order2)[w] = VECTOR(order)[l2]; - VECTOR(d2)[w]=d[l2]; + VECTOR(order2)[w] = VECTOR(order)[l2]; + VECTOR(d2)[w]=d[l2]; w++; l2--; } } @@ -627,13 +627,13 @@ int igraph_arpack_rssort(igraph_vector_t *values, igraph_matrix_t *vectors, igraph_vector_copy_to(&d2, d); igraph_vector_destroy(&order2); igraph_vector_destroy(&d2); - IGRAPH_FINALLY_CLEAN(2); + IGRAPH_FINALLY_CLEAN(2); } #undef which /* Copy values */ - if (values) { + if (values) { IGRAPH_CHECK(igraph_vector_resize(values, nans)); memcpy(VECTOR(*values), d, sizeof(igraph_real_t) * nans); } @@ -642,13 +642,13 @@ int igraph_arpack_rssort(igraph_vector_t *values, igraph_matrix_t *vectors, if (vectors) { int i; IGRAPH_CHECK(igraph_matrix_resize(vectors, n, nans)); - for (i=0; in; int nconv=options->nconv; int nev=options->nev; @@ -685,7 +685,7 @@ int igraph_arpack_rnsort(igraph_matrix_t *values, igraph_matrix_t *vectors, } #undef which - + IGRAPH_CHECK(igraph_vector_init_seq(&order, 0, nconv-1)); IGRAPH_FINALLY(igraph_vector_destroy, &order); #ifdef HAVE_GFORTRAN @@ -701,25 +701,39 @@ int igraph_arpack_rnsort(igraph_matrix_t *values, igraph_matrix_t *vectors, } if (vectors) { - int i, nc=0, nr=0, ncol, wh=0, vx=0; + int nc=0, nr=0, ncol, vx=0; for (i=0; i 0) { + /* complex eigenvalue, positive imaginary part encountered first. + * ARPACK stores its eigenvector directly in two consecutive columns. + * The complex conjugate pair of the eigenvalue (if any) will be in + * the next column and we will skip it because we advance 'i' below */ + memcpy(&MATRIX(*vectors, 0, vx), v + n * idx, sizeof(igraph_real_t) * 2 * n); + vx += 2; + i++; } else { - wh=1-wh; + /* complex eigenvalue, negative imaginary part encountered first. + * The positive one will be the next one, but we need to copy the + * eigenvector corresponding to the eigenvalue with the positive + * imaginary part. */ + idx = (unsigned int) VECTOR(order)[i+1]; + memcpy(&MATRIX(*vectors, 0, vx), v + n * idx, sizeof(igraph_real_t) * 2 * n); + vx += 2; + i++; } } } @@ -727,6 +741,29 @@ int igraph_arpack_rnsort(igraph_matrix_t *values, igraph_matrix_t *vectors, igraph_vector_destroy(&order); IGRAPH_FINALLY_CLEAN(1); + if (values) { + /* Strive to include complex conjugate eigenvalue pairs in a way that the + * positive imaginary part comes first */ + for (i = 0; i < nans; i++) { + if (MATRIX(*values, i, 1) == 0) { + /* Real eigenvalue, nothing to do */ + } else if (MATRIX(*values, i, 1) < 0) { + /* Negative imaginary part came first; negate the imaginary part for + * this eigenvalue and the next one (which is the complex conjugate + * pair), and skip it */ + MATRIX(*values, i, 1) *= -1; + i++; + if (i < nans) { + MATRIX(*values, i, 1) *= -1; + } + } else { + /* Positive imaginary part; skip the next eigenvalue, which is the + * complex conjugate pair */ + i++; + } + } + } + return 0; } @@ -758,8 +795,8 @@ void igraph_i_arpack_auto_ncv(igraph_arpack_options_t* options) { options->ncv = min_ncv; } /* ...but at most n-1 */ - if (options->ncv > options->n) { - options->ncv = options->n; + if (options->ncv > options->n - 1) { + options->ncv = options->n - 1; } } @@ -788,9 +825,9 @@ void igraph_i_arpack_report_no_convergence(const igraph_arpack_options_t* option * \param options An \ref igraph_arpack_options_t object. * \param storage An \ref igraph_arpack_storage_t object, or a null * pointer. In the latter case memory allocation and deallocation - * is performed automatically. Either this or the \p vectors argument - * must be non-null if the ARPACK iteration is started from a - * given starting vector. If both are given \p vectors take + * is performed automatically. Either this or the \p vectors argument + * must be non-null if the ARPACK iteration is started from a + * given starting vector. If both are given \p vectors take * precedence. * \param values If not a null pointer, then it should be a pointer to an * initialized vector. The eigenvalues will be stored here. The @@ -814,7 +851,7 @@ int 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) { - + igraph_real_t *v, *workl, *workd, *d, *resid, *ax; igraph_bool_t free_them=0; int *select, i; @@ -828,10 +865,10 @@ int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra, char origwhich[2]={ options->which[0], options->which[1] }; igraph_real_t origtol=options->tol; - /* Special case for 1x1 and 2x2 matrices */ - if (options->n == 1) { + /* Special case for 1x1 and 2x2 matrices in mode 1 */ + if (options->mode == 1 && options->n == 1) { return igraph_i_arpack_rssolve_1x1(fun, extra, options, values, vectors); - } else if (options->n == 2) { + } else if (options->mode == 1 && options->n == 2) { return igraph_i_arpack_rssolve_2x2(fun, extra, options, values, vectors); } @@ -842,7 +879,7 @@ int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra, } if (options->lworkl == 0) { options->lworkl=options->ncv*(options->ncv+8); } if (options->which[0] == 'X') { options->which[0]='L'; options->which[1]='M'; } - + if (storage) { /* Storage provided */ if (storage->maxn < options->n) { @@ -862,11 +899,11 @@ int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra, resid = storage->resid; ax = storage->ax; select = storage->select; - + } else { /* Storage not provided */ free_them=1; - + #define CHECKMEM(x) \ if (!x) { \ IGRAPH_ERROR("Cannot allocate memory for ARPACK", IGRAPH_ENOMEM); \ @@ -886,17 +923,24 @@ int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra, } /* Set final bits */ + options->bmat[0]='I'; options->iparam[0]=options->ishift; + options->iparam[1]=0; // not referenced options->iparam[2]=options->mxiter; - options->iparam[3]=options->nb; + options->iparam[3]=1; // currently dsaupd() works only for nb=1 options->iparam[4]=0; + options->iparam[5]=0; // not referenced options->iparam[6]=options->mode; + options->iparam[7]=0; // return value + options->iparam[8]=0; // return value + options->iparam[9]=0; // return value + options->iparam[10]=0; // return value options->info=options->start; if (options->start) { if (!storage && !vectors) { IGRAPH_ERROR("Starting vector not given", IGRAPH_EINVAL); } - if (vectors && (igraph_matrix_nrow(vectors) != options->n || + if (vectors && (igraph_matrix_nrow(vectors) != options->n || igraph_matrix_ncol(vectors) != 1)) { IGRAPH_ERROR("Invalid starting vector size", IGRAPH_EINVAL); } @@ -931,19 +975,19 @@ int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra, IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product", IGRAPH_ARPACK_PROD); } - + } else { break; } } - + if (options->info == 1) { igraph_i_arpack_report_no_convergence(options); } if (options->info != 0) { IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dsaupd(options->info)); } - + options->ierr=0; #ifdef HAVE_GFORTRAN igraphdseupd_(&rvec, all, select, d, v, &options->ldv, @@ -965,9 +1009,9 @@ int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra, if (options->ierr != 0) { IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dseupd(options->ierr)); } - + /* Save the result */ - + options->noiter=options->iparam[2]; options->nconv=options->iparam[4]; options->numop=options->iparam[8]; @@ -979,17 +1023,17 @@ int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra, "solver"); } - if (values || vectors) { + if (values || vectors) { IGRAPH_CHECK(igraph_arpack_rssort(values, vectors, options, d, v)); } - + options->ldv=origldv; options->ncv=origncv; options->lworkl=origlworkl; options->which[0] = origwhich[0]; options->which[1] = origwhich[1]; options->tol=origtol; options->nev=orignev; - + /* Clean up if needed */ if (free_them) { igraph_Free(select); @@ -1031,6 +1075,17 @@ int igraph_arpack_rssolve(igraph_arpack_function_t *fun, void *extra, * \param vectors If not a null pointer, then it must be a pointer to * an initialized matrix. The eigenvectors will be stored in the * columns of the matrix. The matrix will be resized as needed. + * Note that real eigenvalues will have real eigenvectors in a single + * column in this matrix; however, complex eigenvalues come in conjugate + * pairs and the result matrix will store the eigenvector corresponding to + * the eigenvalue with \em positive imaginary part only. Since in this case + * the eigenvector is also complex, it will occupy \em two columns in the + * eigenvector matrix (the real and the imaginary parts, in this order). + * Caveat: if the eigenvalue vector returns only the eigenvalue with the + * \em negative imaginary part for a complex conjugate eigenvalue pair, the + * result vector will \em still store the eigenvector corresponding to the + * eigenvalue with the positive imaginary part (since this is how ARPACK + * works). * \return Error code. * * Time complexity: depends on the matrix-vector @@ -1048,21 +1103,21 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, igraph_real_t *v, *workl, *workd, *dr, *di, *resid, *workev; igraph_bool_t free_them=0; int *select, i; - + int ido=0; int rvec= vectors || storage ? 1 : 0; char *all="All"; - + int origldv=options->ldv, origlworkl=options->lworkl, orignev=options->nev, origncv=options->ncv; char origwhich[2]={ options->which[0], options->which[1] }; igraph_real_t origtol=options->tol; int d_size; - /* Special case for 1x1 and 2x2 matrices */ - if (options->n == 1) { + /* Special case for 1x1 and 2x2 matrices in mode 1 */ + if (options->mode == 1 && options->n == 1) { return igraph_i_arpack_rnsolve_1x1(fun, extra, options, values, vectors); - } else if (options->n == 2) { + } else if (options->mode == 1 && options->n == 2) { return igraph_i_arpack_rnsolve_2x2(fun, extra, options, values, vectors); } @@ -1073,7 +1128,7 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, } if (options->lworkl == 0) { options->lworkl=3*options->ncv*(options->ncv+2); } if (options->which[0] == 'X') { options->which[0]='L'; options->which[1]='M'; } - + if (storage) { /* Storage provided */ if (storage->maxn < options->n) { @@ -1095,11 +1150,11 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, d_size = options->n; resid = storage->resid; select = storage->select; - + } else { /* Storage not provided */ free_them=1; - + #define CHECKMEM(x) \ if (!x) { \ IGRAPH_ERROR("Cannot allocate memory for ARPACK", IGRAPH_ENOMEM); \ @@ -1115,17 +1170,24 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, resid=igraph_Calloc(options->n, igraph_real_t); CHECKMEM(resid); select=igraph_Calloc(options->ncv, int); CHECKMEM(select); workev=igraph_Calloc(3*options->ncv, igraph_real_t); CHECKMEM(workev); - + #undef CHECKMEM } - + /* Set final bits */ + options->bmat[0]='I'; options->iparam[0]=options->ishift; + options->iparam[1]=0; // not referenced options->iparam[2]=options->mxiter; - options->iparam[3]=options->nb; + options->iparam[3]=1; // currently dnaupd() works only for nb=1 options->iparam[4]=0; + options->iparam[5]=0; // not referenced options->iparam[6]=options->mode; + options->iparam[7]=0; // return value + options->iparam[8]=0; // return value + options->iparam[9]=0; // return value + options->iparam[10]=0; // return value options->info=options->start; if (options->start) { if (igraph_matrix_nrow(vectors) != options->n || igraph_matrix_ncol(vectors) != 1) { @@ -1135,7 +1197,7 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, resid[i]=MATRIX(*vectors, i, 0); } } - + /* Ok, we have everything */ while (1) { #ifdef HAVE_GFORTRAN @@ -1145,7 +1207,6 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, options->iparam, options->ipntr, workd, workl, &options->lworkl, &options->info, /*bmat_len=*/ 1, /*which_len=*/ 2); - #else igraphdnaupd_(&ido, options->bmat, &options->n, options->which, &options->nev, &options->tol, @@ -1153,7 +1214,7 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, options->iparam, options->ipntr, workd, workl, &options->lworkl, &options->info); #endif - + if (ido==-1 || ido==1) { igraph_real_t *from=workd+options->ipntr[0]-1; igraph_real_t *to=workd+options->ipntr[1]-1; @@ -1161,7 +1222,6 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product", IGRAPH_ARPACK_PROD); } - } else { break; } @@ -1170,7 +1230,7 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, if (options->info == 1) { igraph_i_arpack_report_no_convergence(options); } - if (options->info != 0) { + if (options->info != 0 && options->info != -9999) { IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dnaupd(options->info)); } @@ -1197,7 +1257,7 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, } /* Save the result */ - + options->noiter=options->iparam[2]; options->nconv=options->iparam[4]; options->numop=options->iparam[8]; @@ -1209,18 +1269,20 @@ int igraph_arpack_rnsolve(igraph_arpack_function_t *fun, void *extra, "solver"); } - if (values || vectors) { - IGRAPH_CHECK(igraph_arpack_rnsort(values, vectors, options, - dr, di, v)); - } - + /* ARPACK might modify stuff in 'options' so reset everything that could + * potentially get modified */ options->ldv=origldv; options->ncv=origncv; options->lworkl=origlworkl; options->which[0] = origwhich[0]; options->which[1] = origwhich[1]; options->tol=origtol; options->nev=orignev; - + + if (values || vectors) { + IGRAPH_CHECK(igraph_arpack_rnsort(values, vectors, options, + dr, di, v)); + } + /* Clean up if needed */ if (free_them) { igraph_Free(workev); @@ -1288,7 +1350,7 @@ int igraph_arpack_unpack_complex(igraph_matrix_t *vectors, igraph_matrix_t *valu for (i=nev; i