/* -*- mode: C -*- */ /* vim:set ts=4 sw=4 sts=4 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 02110-1301 USA */ #include "igraph_arpack.h" #include "igraph_arpack_internal.h" #include "igraph_memory.h" #include #include #include /* The ARPACK example file dssimp.f is used as a template */ int igraph_i_arpack_err_dsaupd(int error) { switch (error) { case 1: return IGRAPH_ARPACK_MAXIT; case 3: return IGRAPH_ARPACK_NOSHIFT; case -1: return IGRAPH_ARPACK_NPOS; case -2: return IGRAPH_ARPACK_NEVNPOS; case -3: return IGRAPH_ARPACK_NCVSMALL; case -4: return IGRAPH_ARPACK_NONPOSI; case -5: return IGRAPH_ARPACK_WHICHINV; case -6: return IGRAPH_ARPACK_BMATINV; case -7: return IGRAPH_ARPACK_WORKLSMALL; case -8: return IGRAPH_ARPACK_TRIDERR; case -9: return IGRAPH_ARPACK_ZEROSTART; case -10: return IGRAPH_ARPACK_MODEINV; case -11: return IGRAPH_ARPACK_MODEBMAT; case -12: return IGRAPH_ARPACK_ISHIFT; case -13: return IGRAPH_ARPACK_NEVBE; case -9999: return IGRAPH_ARPACK_NOFACT; default: return IGRAPH_ARPACK_UNKNOWN; } } int igraph_i_arpack_err_dseupd(int error) { switch (error) { case -1: return IGRAPH_ARPACK_NPOS; case -2: return IGRAPH_ARPACK_NEVNPOS; case -3: return IGRAPH_ARPACK_NCVSMALL; case -5: return IGRAPH_ARPACK_WHICHINV; case -6: return IGRAPH_ARPACK_BMATINV; case -7: return IGRAPH_ARPACK_WORKLSMALL; case -8: return IGRAPH_ARPACK_TRIDERR; case -9: return IGRAPH_ARPACK_ZEROSTART; case -10: return IGRAPH_ARPACK_MODEINV; case -11: return IGRAPH_ARPACK_MODEBMAT; case -12: return IGRAPH_ARPACK_NEVBE; case -14: return IGRAPH_ARPACK_FAILED; case -15: return IGRAPH_ARPACK_HOWMNY; case -16: return IGRAPH_ARPACK_HOWMNYS; case -17: return IGRAPH_ARPACK_EVDIFF; default: return IGRAPH_ARPACK_UNKNOWN; } } int igraph_i_arpack_err_dnaupd(int error) { switch (error) { case 1: return IGRAPH_ARPACK_MAXIT; case 3: return IGRAPH_ARPACK_NOSHIFT; case -1: return IGRAPH_ARPACK_NPOS; case -2: return IGRAPH_ARPACK_NEVNPOS; case -3: return IGRAPH_ARPACK_NCVSMALL; case -4: return IGRAPH_ARPACK_NONPOSI; case -5: return IGRAPH_ARPACK_WHICHINV; case -6: return IGRAPH_ARPACK_BMATINV; case -7: return IGRAPH_ARPACK_WORKLSMALL; case -8: return IGRAPH_ARPACK_TRIDERR; case -9: return IGRAPH_ARPACK_ZEROSTART; case -10: return IGRAPH_ARPACK_MODEINV; case -11: return IGRAPH_ARPACK_MODEBMAT; case -12: return IGRAPH_ARPACK_ISHIFT; case -9999: return IGRAPH_ARPACK_NOFACT; default: return IGRAPH_ARPACK_UNKNOWN; } } int igraph_i_arpack_err_dneupd(int error) { switch (error) { case 1: return IGRAPH_ARPACK_REORDER; case -1: return IGRAPH_ARPACK_NPOS; case -2: return IGRAPH_ARPACK_NEVNPOS; case -3: return IGRAPH_ARPACK_NCVSMALL; case -5: return IGRAPH_ARPACK_WHICHINV; case -6: return IGRAPH_ARPACK_BMATINV; case -7: return IGRAPH_ARPACK_WORKLSMALL; case -8: return IGRAPH_ARPACK_SHUR; case -9: return IGRAPH_ARPACK_LAPACK; case -10: return IGRAPH_ARPACK_MODEINV; case -11: return IGRAPH_ARPACK_MODEBMAT; case -12: return IGRAPH_ARPACK_HOWMNYS; case -13: return IGRAPH_ARPACK_HOWMNY; case -14: return IGRAPH_ARPACK_FAILED; case -15: return IGRAPH_ARPACK_EVDIFF; default: return IGRAPH_ARPACK_UNKNOWN; } } /** * \function igraph_arpack_options_init * Initialize ARPACK options * * Initializes ARPACK options, set them to default values. * You can always pass the initialized \ref igraph_arpack_options_t * object to built-in igraph functions without any modification. The * built-in igraph functions modify the options to perform their * calculation, e.g. \ref igraph_pagerank() always searches for the * eigenvalue with the largest magnitude, regardless of the supplied * value. * * If you want to implement your own function involving eigenvalue * calculation using ARPACK, however, you will likely need to set up * the fields for yourself. * \param o The \ref igraph_arpack_options_t object to initialize. * * Time complexity: O(1). */ void igraph_arpack_options_init(igraph_arpack_options_t *o) { o->bmat[0] = 'I'; o->n = 0; /* needs to be updated! */ o->which[0] = 'X'; o->which[1] = 'X'; o->nev = 1; o->tol = 0; o->ncv = 0; /* 0 means "automatic" */ o->ldv = o->n; /* will be updated to (real) n */ o->ishift = 1; o->mxiter = 3000; o->nb = 1; o->mode = 1; o->start = 0; o->lworkl = 0; 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; } /** * \function igraph_arpack_storage_init * Initialize ARPACK storage * * You only need this function if you want to run multiple eigenvalue * calculations using ARPACK, and want to spare the memory * allocation/deallocation between each two runs. Otherwise it is safe * to supply a null pointer as the \c storage argument of both \ref * igraph_arpack_rssolve() and \ref igraph_arpack_rnsolve() to make * memory allocated and deallocated automatically. * * Don't forget to call the \ref * igraph_arpack_storage_destroy() function on the storage object if * you don't need it any more. * \param s The \ref igraph_arpack_storage_t object to initialize. * \param maxn The maximum order of the matrices. * \param maxncv The maximum NCV parameter intended to use. * \param maxldv The maximum LDV parameter intended to use. * \param symm Whether symmetric or non-symmetric problems will be * solved using this \ref igraph_arpack_storage_t. (You cannot use * the same storage both with symmetric and non-symmetric solvers.) * \return Error code. * * Time complexity: O(maxncv*(maxldv+maxn)). */ 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; s->maxldv = (int) maxldv; #define CHECKMEM(x) \ if (!x) { \ IGRAPH_ERROR("Cannot allocate memory for ARPACK", IGRAPH_ENOMEM); \ } \ IGRAPH_FINALLY(igraph_free, x); s->v = igraph_Calloc(maxldv * maxncv, igraph_real_t); CHECKMEM(s->v); s->workd = igraph_Calloc(3 * maxn, igraph_real_t); CHECKMEM(s->workd); s->d = igraph_Calloc(2 * maxncv, igraph_real_t); CHECKMEM(s->d); s->resid = igraph_Calloc(maxn, igraph_real_t); CHECKMEM(s->resid); s->ax = igraph_Calloc(maxn, igraph_real_t); CHECKMEM(s->ax); s->select = igraph_Calloc(maxncv, int); CHECKMEM(s->select); if (symm) { s->workl = igraph_Calloc(maxncv * (maxncv + 8), igraph_real_t); CHECKMEM(s->workl); s->di = 0; s->workev = 0; } else { s->workl = igraph_Calloc(3 * maxncv * (maxncv + 2), igraph_real_t); CHECKMEM(s->workl); s->di = igraph_Calloc(2 * maxncv, igraph_real_t); CHECKMEM(s->di); s->workev = igraph_Calloc(3 * maxncv, igraph_real_t); CHECKMEM(s->workev); IGRAPH_FINALLY_CLEAN(2); } #undef CHECKMEM IGRAPH_FINALLY_CLEAN(7); return 0; } /** * \function igraph_arpack_storage_destroy * Deallocate ARPACK storage * * \param s The \ref igraph_arpack_storage_t object for which the * memory will be deallocated. * * Time complexity: operating system dependent. */ 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); igraph_Free(s->resid); igraph_Free(s->d); igraph_Free(s->workd); igraph_Free(s->v); } /** * "Solver" for 1x1 eigenvalue problems since ARPACK sometimes blows up with * these. */ int 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) { igraph_real_t a, b; int nev = options->nev; if (nev <= 0) { IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_NEVNPOS); } /* Probe the value in the matrix */ a = 1; if (fun(&b, &a, 1, extra)) { IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product", IGRAPH_ARPACK_PROD); } options->nconv = nev; if (values != 0) { IGRAPH_CHECK(igraph_vector_resize(values, 1)); VECTOR(*values)[0] = b; } if (vectors != 0) { IGRAPH_CHECK(igraph_matrix_resize(vectors, 1, 1)); MATRIX(*vectors, 0, 0) = 1; } return IGRAPH_SUCCESS; } /** * "Solver" for 1x1 eigenvalue problems since ARPACK sometimes blows up with * these. */ int 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) { igraph_real_t a, b; int nev = options->nev; if (nev <= 0) { IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_NEVNPOS); } /* Probe the value in the matrix */ a = 1; if (fun(&b, &a, 1, extra)) { IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product", IGRAPH_ARPACK_PROD); } options->nconv = nev; if (values != 0) { IGRAPH_CHECK(igraph_matrix_resize(values, 1, 2)); MATRIX(*values, 0, 0) = b; MATRIX(*values, 0, 1) = 0; } if (vectors != 0) { IGRAPH_CHECK(igraph_matrix_resize(vectors, 1, 1)); MATRIX(*vectors, 0, 0) = 1; } return IGRAPH_SUCCESS; } /** * "Solver" for 2x2 nonsymmetric eigenvalue problems since ARPACK sometimes * blows up with these. */ int 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) { igraph_real_t vec[2], mat[4]; igraph_real_t a, b, c, d; igraph_real_t trace, det, tsq4_minus_d; igraph_complex_t eval1, eval2; igraph_complex_t evec1[2], evec2[2]; igraph_bool_t swap_evals = 0; igraph_bool_t complex_evals = 0; int nev = options->nev; if (nev <= 0) { IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_NEVNPOS); } if (nev > 2) { nev = 2; } /* Probe the values in the matrix */ vec[0] = 1; vec[1] = 0; if (fun(mat, vec, 2, extra)) { IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product", IGRAPH_ARPACK_PROD); } vec[0] = 0; vec[1] = 1; if (fun(mat + 2, vec, 2, extra)) { IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product", IGRAPH_ARPACK_PROD); } a = mat[0]; b = mat[2]; c = mat[1]; d = mat[3]; /* Get the trace and the determinant */ trace = a + d; det = a * d - b * c; tsq4_minus_d = trace * trace / 4 - det; /* Calculate the eigenvalues */ complex_evals = tsq4_minus_d < 0; eval1 = igraph_complex_sqrt_real(tsq4_minus_d); if (complex_evals) { eval2 = igraph_complex_mul_real(eval1, -1); } else { /* to avoid having -0 in the imaginary part */ eval2 = igraph_complex(-IGRAPH_REAL(eval1), 0); } eval1 = igraph_complex_add_real(eval1, trace / 2); eval2 = igraph_complex_add_real(eval2, trace / 2); if (c != 0) { evec1[0] = igraph_complex_sub_real(eval1, d); evec1[1] = igraph_complex(c, 0); evec2[0] = igraph_complex_sub_real(eval2, d); evec2[1] = igraph_complex(c, 0); } else if (b != 0) { evec1[0] = igraph_complex(b, 0); evec1[1] = igraph_complex_sub_real(eval1, a); evec2[0] = igraph_complex(b, 0); evec2[1] = igraph_complex_sub_real(eval2, a); } else { evec1[0] = igraph_complex(1, 0); evec1[1] = igraph_complex(0, 0); evec2[0] = igraph_complex(0, 0); evec2[1] = igraph_complex(1, 0); } /* Sometimes we have to swap eval1 with eval2 and evec1 with eval2; * determine whether we have to do it now */ if (options->which[0] == 'S') { if (options->which[1] == 'M') { /* eval1 must be the one with the smallest magnitude */ swap_evals = (igraph_complex_mod(eval1) > igraph_complex_mod(eval2)); } else if (options->which[1] == 'R') { /* eval1 must be the one with the smallest real part */ swap_evals = (IGRAPH_REAL(eval1) > IGRAPH_REAL(eval2)); } else if (options->which[1] == 'I') { /* eval1 must be the one with the smallest imaginary part */ swap_evals = (IGRAPH_IMAG(eval1) > IGRAPH_IMAG(eval2)); } else { IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_WHICHINV); } } else if (options->which[0] == 'L') { if (options->which[1] == 'M') { /* eval1 must be the one with the largest magnitude */ swap_evals = (igraph_complex_mod(eval1) < igraph_complex_mod(eval2)); } else if (options->which[1] == 'R') { /* eval1 must be the one with the largest real part */ swap_evals = (IGRAPH_REAL(eval1) < IGRAPH_REAL(eval2)); } else if (options->which[1] == 'I') { /* eval1 must be the one with the largest imaginary part */ swap_evals = (IGRAPH_IMAG(eval1) < IGRAPH_IMAG(eval2)); } else { IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_WHICHINV); } } else if (options->which[0] == 'X' && options->which[1] == 'X') { /* No preference on the ordering of eigenvectors */ } else { /* fprintf(stderr, "%c%c\n", options->which[0], options->which[1]); */ IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_WHICHINV); } options->nconv = nev; if (swap_evals) { igraph_complex_t dummy; dummy = eval1; eval1 = eval2; eval2 = dummy; dummy = evec1[0]; evec1[0] = evec2[0]; evec2[0] = dummy; dummy = evec1[1]; evec1[1] = evec2[1]; evec2[1] = dummy; } if (complex_evals) { /* The eigenvalues are conjugate pairs, so we store only the * one with positive imaginary part */ if (IGRAPH_IMAG(eval1) < 0) { eval1 = eval2; evec1[0] = evec2[0]; evec1[1] = evec2[1]; } } if (values != 0) { IGRAPH_CHECK(igraph_matrix_resize(values, nev, 2)); MATRIX(*values, 0, 0) = IGRAPH_REAL(eval1); MATRIX(*values, 0, 1) = IGRAPH_IMAG(eval1); if (nev > 1) { MATRIX(*values, 1, 0) = IGRAPH_REAL(eval2); MATRIX(*values, 1, 1) = IGRAPH_IMAG(eval2); } } if (vectors != 0) { if (complex_evals) { IGRAPH_CHECK(igraph_matrix_resize(vectors, 2, 2)); MATRIX(*vectors, 0, 0) = IGRAPH_REAL(evec1[0]); MATRIX(*vectors, 1, 0) = IGRAPH_REAL(evec1[1]); MATRIX(*vectors, 0, 1) = IGRAPH_IMAG(evec1[0]); MATRIX(*vectors, 1, 1) = IGRAPH_IMAG(evec1[1]); } else { IGRAPH_CHECK(igraph_matrix_resize(vectors, 2, nev)); MATRIX(*vectors, 0, 0) = IGRAPH_REAL(evec1[0]); MATRIX(*vectors, 1, 0) = IGRAPH_REAL(evec1[1]); if (nev > 1) { MATRIX(*vectors, 0, 1) = IGRAPH_REAL(evec2[0]); MATRIX(*vectors, 1, 1) = IGRAPH_REAL(evec2[1]); } } } return IGRAPH_SUCCESS; } /** * "Solver" for symmetric 2x2 eigenvalue problems since ARPACK sometimes blows * up with these. */ int 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) { igraph_real_t vec[2], mat[4]; igraph_real_t a, b, c, d; igraph_real_t trace, det, tsq4_minus_d; igraph_real_t eval1, eval2; int nev = options->nev; if (nev <= 0) { IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_NEVNPOS); } if (nev > 2) { nev = 2; } /* Probe the values in the matrix */ vec[0] = 1; vec[1] = 0; if (fun(mat, vec, 2, extra)) { IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product", IGRAPH_ARPACK_PROD); } vec[0] = 0; vec[1] = 1; if (fun(mat + 2, vec, 2, extra)) { IGRAPH_ERROR("ARPACK error while evaluating matrix-vector product", IGRAPH_ARPACK_PROD); } a = mat[0]; b = mat[2]; c = mat[1]; d = mat[3]; /* Get the trace and the determinant */ trace = a + d; det = a * d - b * c; tsq4_minus_d = trace * trace / 4 - det; if (tsq4_minus_d >= 0) { /* Both eigenvalues are real */ eval1 = trace / 2 + sqrt(tsq4_minus_d); eval2 = trace / 2 - sqrt(tsq4_minus_d); if (c != 0) { mat[0] = eval1 - d; mat[2] = eval2 - d; mat[1] = c; mat[3] = c; } else if (b != 0) { mat[0] = b; mat[2] = b; mat[1] = eval1 - a; mat[3] = eval2 - a; } else { mat[0] = 1; mat[2] = 0; mat[1] = 0; mat[3] = 1; } } else { /* Both eigenvalues are complex. Should not happen with symmetric * matrices. */ IGRAPH_ERROR("ARPACK error, 2x2 matrix is not symmetric", IGRAPH_EINVAL); } /* eval1 is always the larger eigenvalue. If we want the smaller * one, we have to swap eval1 with eval2 and also the columns of mat */ if (options->which[0] == 'S') { trace = eval1; eval1 = eval2; eval2 = trace; trace = mat[0]; mat[0] = mat[2]; mat[2] = trace; trace = mat[1]; mat[1] = mat[3]; mat[3] = trace; } else if (options->which[0] == 'L' || options->which[0] == 'B') { /* Nothing to do here */ } else if (options->which[0] == 'X' && options->which[1] == 'X') { /* No preference on the ordering of eigenvectors */ } else { IGRAPH_ERROR("ARPACK error", IGRAPH_ARPACK_WHICHINV); } options->nconv = nev; if (values != 0) { IGRAPH_CHECK(igraph_vector_resize(values, nev)); VECTOR(*values)[0] = eval1; if (nev > 1) { VECTOR(*values)[1] = eval2; } } if (vectors != 0) { IGRAPH_CHECK(igraph_matrix_resize(vectors, 2, nev)); MATRIX(*vectors, 0, 0) = mat[0]; MATRIX(*vectors, 1, 0) = mat[1]; if (nev > 1) { MATRIX(*vectors, 0, 1) = mat[2]; MATRIX(*vectors, 1, 1) = mat[3]; } } return IGRAPH_SUCCESS; } int 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) { igraph_vector_t order; char sort[2]; int apply = 1; unsigned int n = (unsigned int) options->n; int nconv = options->nconv; int nev = options->nev; unsigned int nans = (unsigned int) (nconv < nev ? nconv : nev); #define which(a,b) (options->which[0]==a && options->which[1]==b) if (which('L', 'A')) { sort[0] = 'S'; sort[1] = 'A'; } else if (which('S', 'A')) { sort[0] = 'L'; sort[1] = 'A'; } else if (which('L', 'M')) { sort[0] = 'S'; sort[1] = 'M'; } else if (which('S', 'M')) { sort[0] = 'L'; sort[1] = 'M'; } else if (which('B', 'E')) { sort[0] = 'L'; sort[1] = 'A'; } IGRAPH_CHECK(igraph_vector_init_seq(&order, 0, nconv - 1)); IGRAPH_FINALLY(igraph_vector_destroy, &order); #ifdef HAVE_GFORTRAN igraphdsortr_(sort, &apply, &nconv, d, VECTOR(order), /*which_len=*/ 2); #else igraphdsortr_(sort, &apply, &nconv, d, VECTOR(order)); #endif /* BE is special */ if (which('B', 'E')) { int w = 0, l1 = 0, l2 = nev - 1; igraph_vector_t order2, d2; 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]; w++; l1++; if (l1 <= l2) { VECTOR(order2)[w] = VECTOR(order)[l2]; VECTOR(d2)[w] = d[l2]; w++; l2--; } } igraph_vector_update(&order, &order2); igraph_vector_copy_to(&d2, d); igraph_vector_destroy(&order2); igraph_vector_destroy(&d2); IGRAPH_FINALLY_CLEAN(2); } #undef which /* Copy values */ if (values) { IGRAPH_CHECK(igraph_vector_resize(values, nans)); memcpy(VECTOR(*values), d, sizeof(igraph_real_t) * nans); } /* Reorder vectors */ if (vectors) { int i; IGRAPH_CHECK(igraph_matrix_resize(vectors, n, nans)); 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); return 0; } 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, igraph_real_t *v) { igraph_vector_t order; char sort[2]; int apply = 1, i; unsigned int n = (unsigned int) options->n; int nconv = options->nconv; int nev = options->nev; unsigned int nans = (unsigned int) (nconv < nev ? nconv : nev); #define which(a,b) (options->which[0]==a && options->which[1]==b) if (which('L', 'M')) { sort[0] = 'S'; sort[1] = 'M'; } else if (which('S', 'M')) { sort[0] = 'L'; sort[1] = 'M'; } else if (which('L', 'R')) { sort[0] = 'S'; sort[1] = 'R'; } else if (which('S', 'R')) { sort[0] = 'L'; sort[1] = 'R'; } else if (which('L', 'I')) { sort[0] = 'S'; sort[1] = 'I'; } else if (which('S', 'I')) { sort[0] = 'L'; sort[1] = 'I'; } #undef which IGRAPH_CHECK(igraph_vector_init_seq(&order, 0, nconv - 1)); IGRAPH_FINALLY(igraph_vector_destroy, &order); #ifdef HAVE_GFORTRAN igraphdsortc_(sort, &apply, &nconv, dr, di, VECTOR(order), /*which_len=*/ 2); #else igraphdsortc_(sort, &apply, &nconv, dr, di, VECTOR(order)); #endif if (values) { IGRAPH_CHECK(igraph_matrix_resize(values, nans, 2)); memcpy(&MATRIX(*values, 0, 0), dr, sizeof(igraph_real_t) * nans); memcpy(&MATRIX(*values, 0, 1), di, sizeof(igraph_real_t) * nans); } if (vectors) { 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; 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 (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 { /* 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++; } } } 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; } /** * \function igraph_i_arpack_auto_ncv * \brief Tries to set up the value of \c ncv in an \c igraph_arpack_options_t * automagically. */ void igraph_i_arpack_auto_ncv(igraph_arpack_options_t* options) { /* This is similar to how Octave determines the value of ncv, with some * modifications. */ int min_ncv = options->nev * 2 + 1; /* Use twice the number of desired eigenvectors plus one by default */ options->ncv = min_ncv; /* ...but use at least 20 Lanczos vectors... */ if (options->ncv < 20) { options->ncv = 20; } /* ...but having ncv close to n leads to some problems with small graphs * (example: PageRank of "A <--> C, D <--> E, B"), so we don't let it * to be larger than n / 2... */ if (options->ncv > options->n / 2) { options->ncv = options->n / 2; } /* ...but we need at least min_ncv. */ if (options->ncv < min_ncv) { options->ncv = min_ncv; } /* ...but at most n */ if (options->ncv > options->n) { options->ncv = options->n; } } /** * \function igraph_i_arpack_report_no_convergence * \brief Prints a warning that informs the user that the ARPACK solver * did not converge. */ void igraph_i_arpack_report_no_convergence(const igraph_arpack_options_t* options) { char buf[1024]; snprintf(buf, sizeof(buf), "ARPACK solver failed to converge (%d iterations, " "%d/%d eigenvectors converged)", options->iparam[2], options->iparam[4], options->nev); IGRAPH_WARNING(buf); } /** * \function igraph_arpack_rssolve * \brief ARPACK solver for symmetric matrices * * This is the ARPACK solver for symmetric matrices. Please use * \ref igraph_arpack_rnsolve() for non-symmetric matrices. * \param fun Pointer to an \ref igraph_arpack_function_t object, * the function that performs the matrix-vector multiplication. * \param extra An extra argument to be passed to \c fun. * \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 * 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 * vector will be resized as needed. * \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. * 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. * \return Error code. * * Time complexity: depends on the matrix-vector * multiplication. Usually a small number of iterations is enough, so * if the matrix is sparse and the matrix-vector multiplication can be * done in O(n) time (the number of vertices), then the eigenvalues * are found in O(n) time as well. */ 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; int ido = 0; int rvec = vectors || storage ? 1 : 0; /* calculate eigenvectors? */ 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; /* 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->mode == 1 && options->n == 2) { return igraph_i_arpack_rssolve_2x2(fun, extra, options, values, vectors); } /* Brush up options if needed */ if (options->ldv == 0) { options->ldv = options->n; } if (options->ncv == 0) { igraph_i_arpack_auto_ncv(options); } 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) { IGRAPH_ERROR("Not enough storage for ARPACK (`n')", IGRAPH_EINVAL); } if (storage->maxncv < options->ncv) { IGRAPH_ERROR("Not enough storage for ARPACK (`ncv')", IGRAPH_EINVAL); } if (storage->maxldv < options->ldv) { IGRAPH_ERROR("Not enough storage for ARPACK (`ldv')", IGRAPH_EINVAL); } v = storage->v; workl = storage->workl; workd = storage->workd; d = storage->d; 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); \ } \ IGRAPH_FINALLY(igraph_free, x); v = igraph_Calloc(options->ldv * options->ncv, igraph_real_t); CHECKMEM(v); workl = igraph_Calloc(options->lworkl, igraph_real_t); CHECKMEM(workl); workd = igraph_Calloc(3 * options->n, igraph_real_t); CHECKMEM(workd); d = igraph_Calloc(2 * options->ncv, igraph_real_t); CHECKMEM(d); resid = igraph_Calloc(options->n, igraph_real_t); CHECKMEM(resid); ax = igraph_Calloc(options->n, igraph_real_t); CHECKMEM(ax); select = igraph_Calloc(options->ncv, int); CHECKMEM(select); #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] = 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 || igraph_matrix_ncol(vectors) != 1)) { IGRAPH_ERROR("Invalid starting vector size", IGRAPH_EINVAL); } if (vectors) { for (i = 0; i < options->n; i++) { resid[i] = MATRIX(*vectors, i, 0); } } } /* Ok, we have everything */ while (1) { #ifdef HAVE_GFORTRAN igraphdsaupd_(&ido, options->bmat, &options->n, options->which, &options->nev, &options->tol, resid, &options->ncv, v, &options->ldv, options->iparam, options->ipntr, workd, workl, &options->lworkl, &options->info, /*bmat_len=*/ 1, /*which_len=*/ 2); #else igraphdsaupd_(&ido, options->bmat, &options->n, options->which, &options->nev, &options->tol, resid, &options->ncv, v, &options->ldv, 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; if (fun(to, from, options->n, extra) != 0) { 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, &options->sigma, options->bmat, &options->n, options->which, &options->nev, &options->tol, resid, &options->ncv, v, &options->ldv, options->iparam, options->ipntr, workd, workl, &options->lworkl, &options->ierr, /*howmny_len=*/ 1, /*bmat_len=*/ 1, /*which_len=*/ 2); #else igraphdseupd_(&rvec, all, select, d, v, &options->ldv, &options->sigma, options->bmat, &options->n, options->which, &options->nev, &options->tol, resid, &options->ncv, v, &options->ldv, options->iparam, options->ipntr, workd, workl, &options->lworkl, &options->ierr); #endif 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]; options->numopb = options->iparam[9]; options->numreo = options->iparam[10]; if (options->nconv < options->nev) { IGRAPH_WARNING("Not enough eigenvalues/vectors in symmetric ARPACK " "solver"); } 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); igraph_Free(ax); igraph_Free(resid); igraph_Free(d); igraph_Free(workd); igraph_Free(workl); igraph_Free(v); IGRAPH_FINALLY_CLEAN(7); } return 0; } /** * \function igraph_arpack_rnsolve * \brief ARPACK solver for non-symmetric matrices * * Please always consider calling \ref igraph_arpack_rssolve() if your * matrix is symmetric, it is much faster. * \ref igraph_arpack_rnsolve() for non-symmetric matrices. * * Note that ARPACK is not called for 2x2 matrices as an exact algebraic * solution exists in these cases. * * \param fun Pointer to an \ref igraph_arpack_function_t object, * the function that performs the matrix-vector multiplication. * \param extra An extra argument to be passed to \c fun. * \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. * \param values If not a null pointer, then it should be a pointer to an * initialized matrix. The (possibly complex) eigenvalues will be * stored here. The matrix will have two columns, the first column * contains the real, the second the imaginary parts of the * eigenvalues. * The matrix will be resized as needed. * \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 * multiplication. Usually a small number of iterations is enough, so * if the matrix is sparse and the matrix-vector multiplication can be * done in O(n) time (the number of vertices), then the eigenvalues * are found in O(n) time as well. */ int 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) { 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 in mode 1 */ if (options->mode == 1 && options->n == 1) { return igraph_i_arpack_rnsolve_1x1(fun, extra, options, values, vectors); } else if (options->mode == 1 && options->n == 2) { return igraph_i_arpack_rnsolve_2x2(fun, extra, options, values, vectors); } /* Brush up options if needed */ if (options->ldv == 0) { options->ldv = options->n; } if (options->ncv == 0) { igraph_i_arpack_auto_ncv(options); } 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) { IGRAPH_ERROR("Not enough storage for ARPACK (`n')", IGRAPH_EINVAL); } if (storage->maxncv < options->ncv) { IGRAPH_ERROR("Not enough storage for ARPACK (`ncv')", IGRAPH_EINVAL); } if (storage->maxldv < options->ldv) { IGRAPH_ERROR("Not enough storage for ARPACK (`ldv')", IGRAPH_EINVAL); } v = storage->v; workl = storage->workl; workd = storage->workd; workev = storage->workev; dr = storage->d; di = storage->di; 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); \ } \ IGRAPH_FINALLY(igraph_free, x); v = igraph_Calloc(options->n * options->ncv, igraph_real_t); CHECKMEM(v); workl = igraph_Calloc(options->lworkl, igraph_real_t); CHECKMEM(workl); workd = igraph_Calloc(3 * options->n, igraph_real_t); CHECKMEM(workd); d_size = 2 * options->nev + 1 > options->ncv ? 2 * options->nev + 1 : options->ncv; dr = igraph_Calloc(d_size, igraph_real_t); CHECKMEM(dr); di = igraph_Calloc(d_size, igraph_real_t); CHECKMEM(di); 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] = 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) { IGRAPH_ERROR("Invalid starting vector size", IGRAPH_EINVAL); } for (i = 0; i < options->n; i++) { resid[i] = MATRIX(*vectors, i, 0); } } /* Ok, we have everything */ while (1) { #ifdef HAVE_GFORTRAN igraphdnaupd_(&ido, options->bmat, &options->n, options->which, &options->nev, &options->tol, resid, &options->ncv, v, &options->ldv, 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, resid, &options->ncv, v, &options->ldv, 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; if (fun(to, from, options->n, extra) != 0) { 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 && options->info != -9999) { IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dnaupd(options->info)); } options->ierr = 0; #ifdef HAVE_GFORTRAN igraphdneupd_(&rvec, all, select, dr, di, v, &options->ldv, &options->sigma, &options->sigmai, workev, options->bmat, &options->n, options->which, &options->nev, &options->tol, resid, &options->ncv, v, &options->ldv, options->iparam, options->ipntr, workd, workl, &options->lworkl, &options->ierr, /*howmny_len=*/ 1, /*bmat_len=*/ 1, /*which_len=*/ 2); #else igraphdneupd_(&rvec, all, select, dr, di, v, &options->ldv, &options->sigma, &options->sigmai, workev, options->bmat, &options->n, options->which, &options->nev, &options->tol, resid, &options->ncv, v, &options->ldv, options->iparam, options->ipntr, workd, workl, &options->lworkl, &options->ierr); #endif if (options->ierr != 0) { IGRAPH_ERROR("ARPACK error", igraph_i_arpack_err_dneupd(options->info)); } /* Save the result */ options->noiter = options->iparam[2]; options->nconv = options->iparam[4]; options->numop = options->iparam[8]; options->numopb = options->iparam[9]; options->numreo = options->iparam[10]; if (options->nconv < options->nev) { IGRAPH_WARNING("Not enough eigenvalues/vectors in ARPACK " "solver"); } /* 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); igraph_Free(select); igraph_Free(resid); igraph_Free(di); igraph_Free(dr); igraph_Free(workd); igraph_Free(workl); igraph_Free(v); IGRAPH_FINALLY_CLEAN(8); } return 0; } /** * \function igraph_arpack_unpack_complex * \brief Make the result of the non-symmetric ARPACK solver more readable * * This function works on the output of \ref igraph_arpack_rnsolve and * brushes it up a bit: it only keeps \p nev eigenvalues/vectors and * every eigenvector is stored in two columns of the \p vectors * matrix. * * * The output of the non-symmetric ARPACK solver is somewhat hard to * parse, as real eigenvectors occupy only one column in the matrix, * and the complex conjugate eigenvectors are not stored at all * (usually). The other problem is that the solver might return more * eigenvalues than requested. The common use of this function is to * call it directly after \ref igraph_arpack_rnsolve with its \p * vectors and \p values argument and \c options->nev as \p nev. * \param vectors The eigenvector matrix, as returned by \ref * igraph_arpack_rnsolve. It will be resized, typically it will be * larger. * \param values The eigenvalue matrix, as returned by \ref * igraph_arpack_rnsolve. It will be resized, typically extra, * unneeded rows (=eigenvalues) will be removed. * \param nev The number of eigenvalues/vectors to keep. Can be less * or equal than the number originally requested from ARPACK. * \return Error code. * * Time complexity: linear in the number of elements in the \p vectors * matrix. */ int igraph_arpack_unpack_complex(igraph_matrix_t *vectors, igraph_matrix_t *values, long int nev) { long int nodes = igraph_matrix_nrow(vectors); long int no_evs = igraph_matrix_nrow(values); long int i, j, k, wh; size_t colsize = (unsigned) nodes * sizeof(igraph_real_t); /* Error checks */ if (nev < 0) { IGRAPH_ERROR("`nev' cannot be negative", IGRAPH_EINVAL); } if (nev > no_evs) { IGRAPH_ERROR("`nev' too large, we don't have that many in `values'", IGRAPH_EINVAL); } IGRAPH_CHECK(igraph_matrix_resize(vectors, nodes, nev * 2)); 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 ???? */ /* real */ j++; } else { /* complex */ if (wh == 0) { j += 2; wh = 1 - wh; } } } j--; /* if (j>=origcol) { */ /* IGRAPH_WARNING("Too few columns in `vectors', ARPACK results are likely wrong"); */ /* } */ /* We copy the j-th eigenvector to the (k-1)-th and k-th column */ k = nev * 2 - 1; for (i = nev - 1; i >= 0; i--) { if (MATRIX(*values, i, 1) == 0) { /* real */ memset( &MATRIX(*vectors, 0, k), 0, colsize); if (k - 1 != j) { memcpy( &MATRIX(*vectors, 0, k - 1), &MATRIX(*vectors, 0, j), colsize); } k -= 2; j -= 1; } else { /* complex */ if (k != j) { /* Separate copy required, otherwise 'from' and 'to' might overlap */ memcpy( &MATRIX(*vectors, 0, k), &MATRIX(*vectors, 0, j), colsize); memcpy( &MATRIX(*vectors, 0, k - 1), &MATRIX(*vectors, 0, j - 1), colsize); } if (i > 1 && MATRIX(*values, i, 1) != -MATRIX(*values, i - 1, 1)) { /* The next one is not a conjugate of this one */ j -= 2; } else { /* Conjugate */ int l; for (l = 0; l < nodes; l++) { MATRIX(*vectors, l, k) = - MATRIX(*vectors, l, k); } } k -= 2; } } return 0; }