--- /dev/null
+From: Tamas Nepusz <ntamas@gmail.com>
+Date: Wed, 1 Aug 2018 19:18:20 +0200
+Origin: https://github.com/szhorvat/igraph/commit/6391d0c.patch
+Bug-Debian: https://bugs.debian.org/902760
+Subject: [PATCH] fixing issues with ARPACK 3.6, related to #1107
+
+Remark: This does not solve the build issue
+
+---
+ examples/simple/igraph_arpack_rnsolve.c | 161 +++++++++++----
+ examples/simple/igraph_arpack_rnsolve.out | 72 ++-----
+ src/arpack.c | 230 ++++++++++++++--------
+ 3 files changed, 277 insertions(+), 186 deletions(-)
+
+--- 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 <csardi.gabor@gmail.com>
+ 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");
+
+ /* -------------- */
+
+--- 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
++===
+--- 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 <csardi.gabor@gmail.com>
+ 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_a
+ 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_a
+ 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_ar
+ }
+
+ #undef CHECKMEM
+-
++
+ IGRAPH_FINALLY_CLEAN(7);
+ return 0;
+ }
+@@ -237,14 +237,14 @@ int igraph_arpack_storage_init(igraph_ar
+ */
+
+ 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_a
+
+ 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_a
+ }
+
+ 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
+ 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
+ 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
+ if (vectors) {
+ int i;
+ IGRAPH_CHECK(igraph_matrix_resize(vectors, n, nans));
+- for (i=0; i<nans; i++) {
++ for (i=0; i<nans; i++) {
+ unsigned int idx=(unsigned int) VECTOR(order)[i];
+ const igraph_real_t *ptr=v + n * idx;
+ memcpy(&MATRIX(*vectors, 0, i), ptr, sizeof(igraph_real_t) * n);
+ }
+ }
+-
++
+ igraph_vector_destroy(&order);
+ IGRAPH_FINALLY_CLEAN(1);
+
+@@ -656,13 +656,13 @@ int igraph_arpack_rssort(igraph_vector_t
+ }
+
+ int 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,
++ const igraph_arpack_options_t *options,
++ igraph_real_t *dr, igraph_real_t *di,
+ igraph_real_t *v) {
+
+ igraph_vector_t order;
+ char sort[2];
+- int apply=1;
++ int apply=1, i;
+ unsigned int n=(unsigned int) options->n;
+ int nconv=options->nconv;
+ int nev=options->nev;
+@@ -685,7 +685,7 @@ int igraph_arpack_rnsort(igraph_matrix_t
+ }
+
+ #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
+ }
+
+ 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<nans; i++) {
+ if (di[i] == 0) { nr++; } else { nc++; }
+ }
+ ncol=(nc/2)*2 + (nc%2)*2 + nr;
+ IGRAPH_CHECK(igraph_matrix_resize(vectors, n, ncol));
++
+ for (i=0; i<nans; i++) {
+- unsigned int idx=(unsigned int) VECTOR(order)[i];
+- igraph_real_t *ptr=v + n * idx;
+- if (di[i]==0) {
+- memcpy(&MATRIX(*vectors, 0, vx), ptr, sizeof(igraph_real_t) * n);
++ unsigned int idx;
++
++ idx = (unsigned int) VECTOR(order)[i];
++
++ if (di[i] == 0) {
++ /* real eigenvalue, single eigenvector */
++ memcpy(&MATRIX(*vectors, 0, vx), v + n * idx, sizeof(igraph_real_t) * n);
+ vx++;
+- } else if (wh==0) {
+- if (di[i] < 0) { ptr -= n; }
+- memcpy(&MATRIX(*vectors, 0, vx), ptr, sizeof(igraph_real_t) * n * 2);
+- wh=1-wh;
+- vx+=2;
++ } else if (di[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
+ 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_arp
+ 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_convergen
+ * \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_
+ 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_
+ 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_
+ }
+ 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_
+ 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_
+ }
+
+ /* 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_
+ 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_
+ 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_
+ "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_
+ * \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_
+ 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_
+ }
+ 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_
+ 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_
+ 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_
+ 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_
+ 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_
+ 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_
+ IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product",
+ IGRAPH_ARPACK_PROD);
+ }
+-
+ } else {
+ break;
+ }
+@@ -1170,7 +1230,7 @@ int igraph_arpack_rnsolve(igraph_arpack_
+ 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_
+ }
+
+ /* 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_
+ "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_
+ for (i=nev; i<igraph_matrix_nrow(values); i++) {
+ IGRAPH_CHECK(igraph_matrix_remove_row(values, i));
+ }
+-
++
+ /* Calculate where to start copying */
+ for (i=0, j=0, wh=0; i<nev; i++) {
+ if (MATRIX(*values,i,1) == 0) { /* TODO: == 0.0 ???? */
+@@ -1333,7 +1395,7 @@ int igraph_arpack_unpack_complex(igraph_
+ /* Conjugate */
+ int l;
+ for (l=0; l<nodes; l++) {
+- MATRIX(*vectors,l,k) = - MATRIX(*vectors,l,k);
++ MATRIX(*vectors,l,k) = - MATRIX(*vectors,l,k);
+ }
+ }
+ k-=2;
projects / igraph.git / commitdiff