Previous Home Next

Chapter 30. Using BLAS, LAPACK and ARPACK for igraph matrices and graphs

1. BLAS interface in igraph
2. LAPACK interface in igraph
3. ARPACK interface in igraph

1.  BLAS interface in igraph

1.1. igraph_blas_ddot — Dot product of two vectors.
1.2. igraph_blas_dnrm2 — Euclidean norm of a vector.
1.3. igraph_blas_dgemv — Matrix-vector multiplication using BLAS, vector version.
1.4. igraph_blas_dgemv_array — Matrix-vector multiplication using BLAS, array version.

BLAS is a highly optimized library for basic linear algebra operations such as vector-vector, matrix-vector and matrix-matrix product. Please see http://www.netlib.org/blas/ for details and a reference implementation in Fortran. igraph contains some wrapper functions that can be used to call BLAS routines in a somewhat more user-friendly way. Not all BLAS routines are included in igraph, and even those which are included might not have wrappers; the extension of the set of wrapped functions will probably be driven by igraph's internal requirements. The wrapper functions usually substitute double-precision floating point arrays used by BLAS with igraph_vector_t and igraph_matrix_t instances and also remove those parameters (such as the number of rows/columns) that can be inferred from the passed arguments directly.

1.1. igraph_blas_ddot — Dot product of two vectors.

int igraph_blas_ddot(const igraph_vector_t *v1, const igraph_vector_t *v2,
                       igraph_real_t *res);

Arguments: 

v1:

The first vector.

v2:

The second vector.

res:

Pointer to a real, the result will be stored here.

Time complexity: O(n) where n is the length of the vectors.

Example 30.1.  File examples/simple/blas.c

#include <igraph.h>

int main() {
    igraph_matrix_t m;
    igraph_vector_t x, y, z;
    igraph_real_t xz, xx;

    igraph_vector_init_real(&x, 3, 1.0, 2.0, 3.0);
    igraph_vector_init_real(&y, 4, 4.0, 5.0, 6.0, 7.0);
    igraph_vector_init_real(&z, 3, -1.0, 0.0, 0.5);

    igraph_matrix_init(&m, 4, 3);
    MATRIX(m, 0, 0) = 1;
    MATRIX(m, 0, 1) = 2;
    MATRIX(m, 0, 2) = 3;
    MATRIX(m, 1, 0) = 2;
    MATRIX(m, 1, 1) = 3;
    MATRIX(m, 1, 2) = 4;
    MATRIX(m, 2, 0) = 3;
    MATRIX(m, 2, 1) = 4;
    MATRIX(m, 2, 2) = 5;
    MATRIX(m, 3, 0) = 4;
    MATRIX(m, 3, 1) = 5;
    MATRIX(m, 3, 2) = 6;

    /* Compute 2 m.x + 3 y and store it in y. */
    igraph_blas_dgemv(/* transpose= */ 0, /* alpha= */ 2, &m, &x, /* beta= */ 3, &y);
    igraph_vector_print(&y);

    /* Compute the squared norm of x, as well as the dor product of x and z. */
    igraph_blas_ddot(&x, &x, &xx);
    igraph_blas_ddot(&x, &z, &xz);
    printf("x.x = %g, x.z = %g\n", xx, xz);

    igraph_matrix_destroy(&m);
    igraph_vector_destroy(&z);
    igraph_vector_destroy(&y);
    igraph_vector_destroy(&x);

    return 0;
}


1.2. igraph_blas_dnrm2 — Euclidean norm of a vector.

igraph_real_t igraph_blas_dnrm2(const igraph_vector_t *v);

Arguments: 

v:

The vector.

Returns: 

Real value, the norm of v.

Time complexity: O(n) where n is the length of the vector.

1.3. igraph_blas_dgemv — Matrix-vector multiplication using BLAS, vector version.

void igraph_blas_dgemv(igraph_bool_t transpose, igraph_real_t alpha,
                       const igraph_matrix_t* a, const igraph_vector_t* x,
                       igraph_real_t beta, igraph_vector_t* y);

This function is a somewhat more user-friendly interface to the dgemv function in BLAS. dgemv performs the operation y = alpha*A*x + beta*y, where x and y are vectors and A is an appropriately sized matrix (symmetric or non-symmetric).

Arguments: 

transpose:

whether to transpose the matrix A

alpha:

the constant alpha

a:

the matrix A

x:

the vector x

beta:

the constant beta

y:

the vector y (which will be modified in-place)

Time complexity: O(nk) if the matrix is of size n x k

See also: 

igraph_blas_dgemv_array if you have arrays instead of vectors.

Example 30.2.  File examples/simple/blas.c

#include <igraph.h>

int main() {
    igraph_matrix_t m;
    igraph_vector_t x, y, z;
    igraph_real_t xz, xx;

    igraph_vector_init_real(&x, 3, 1.0, 2.0, 3.0);
    igraph_vector_init_real(&y, 4, 4.0, 5.0, 6.0, 7.0);
    igraph_vector_init_real(&z, 3, -1.0, 0.0, 0.5);

    igraph_matrix_init(&m, 4, 3);
    MATRIX(m, 0, 0) = 1;
    MATRIX(m, 0, 1) = 2;
    MATRIX(m, 0, 2) = 3;
    MATRIX(m, 1, 0) = 2;
    MATRIX(m, 1, 1) = 3;
    MATRIX(m, 1, 2) = 4;
    MATRIX(m, 2, 0) = 3;
    MATRIX(m, 2, 1) = 4;
    MATRIX(m, 2, 2) = 5;
    MATRIX(m, 3, 0) = 4;
    MATRIX(m, 3, 1) = 5;
    MATRIX(m, 3, 2) = 6;

    /* Compute 2 m.x + 3 y and store it in y. */
    igraph_blas_dgemv(/* transpose= */ 0, /* alpha= */ 2, &m, &x, /* beta= */ 3, &y);
    igraph_vector_print(&y);

    /* Compute the squared norm of x, as well as the dor product of x and z. */
    igraph_blas_ddot(&x, &x, &xx);
    igraph_blas_ddot(&x, &z, &xz);
    printf("x.x = %g, x.z = %g\n", xx, xz);

    igraph_matrix_destroy(&m);
    igraph_vector_destroy(&z);
    igraph_vector_destroy(&y);
    igraph_vector_destroy(&x);

    return 0;
}


1.4. igraph_blas_dgemv_array — Matrix-vector multiplication using BLAS, array version.

void igraph_blas_dgemv_array(igraph_bool_t transpose, igraph_real_t alpha,
                             const igraph_matrix_t* a, const igraph_real_t* x,
                             igraph_real_t beta, igraph_real_t* y);

This function is a somewhat more user-friendly interface to the dgemv function in BLAS. dgemv performs the operation y = alpha*A*x + beta*y, where x and y are vectors and A is an appropriately sized matrix (symmetric or non-symmetric).

Arguments: 

transpose:

whether to transpose the matrix A

alpha:

the constant alpha

a:

the matrix A

x:

the vector x as a regular C array

beta:

the constant beta

y:

the vector y as a regular C array (which will be modified in-place)

Time complexity: O(nk) if the matrix is of size n x k

See also: 

igraph_blas_dgemv if you have vectors instead of arrays.

2.  LAPACK interface in igraph

2.1. Matrix factorization, solving linear systems
2.2. Eigenvalues and eigenvectors of matrices

LAPACK is written in Fortran90 and provides routines for solving systems of simultaneous linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems. The associated matrix factorizations (LU, Cholesky, QR, SVD, Schur, generalized Schur) are also provided, as are related computations such as reordering of the Schur factorizations and estimating condition numbers. Dense and banded matrices are handled, but not general sparse matrices. In all areas, similar functionality is provided for real and complex matrices, in both single and double precision.

igraph provides an interface to a very limited set of LAPACK functions, using the regular igraph data structures.

See more about LAPACK at http://www.netlib.org/lapack/

2.1. Matrix factorization, solving linear systems

2.1.1. igraph_lapack_dgetrf — LU factorization of a general M-by-N matrix.
2.1.2. igraph_lapack_dgetrs — Solve general system of linear equations using LU factorization.
2.1.3. igraph_lapack_dgesv — Solve system of linear equations with LU factorization

2.1.1. igraph_lapack_dgetrf — LU factorization of a general M-by-N matrix.

int igraph_lapack_dgetrf(igraph_matrix_t *a, igraph_vector_int_t *ipiv,
                         int *info);

The factorization has the form A = P * L * U where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n).

Arguments: 

a:

The input/output matrix. On entry, the M-by-N matrix to be factored. On exit, the factors L and U from the factorization A = P * L * U; the unit diagonal elements of L are not stored.

ipiv:

An integer vector, the pivot indices are stored here, unless it is a null pointer. Row i of the matrix was interchanged with row ipiv[i].

info:

LAPACK error code. Zero on successful exit. If its value is a positive number i, it indicates that U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations. If LAPACK returns an error, i.e. a negative info value, then an igraph error is generated as well.

Returns: 

Error code.

Time complexity: TODO.

2.1.2. igraph_lapack_dgetrs — Solve general system of linear equations using LU factorization.

int igraph_lapack_dgetrs(igraph_bool_t transpose, const igraph_matrix_t *a,
                         const igraph_vector_int_t *ipiv, igraph_matrix_t *b);

This function calls LAPACK to solve a system of linear equations A * X = B or A' * X = B with a general N-by-N matrix A using the LU factorization computed by igraph_lapack_dgetrf.

Arguments: 

transpose:

Logical scalar, whether to transpose the input matrix.

a:

A matrix containing the L and U factors from the factorization A = P*L*U. L is expected to be unitriangular, diagonal entries are those of U. If A is singular, no warning or error wil be given and random output will be returned.

ipiv:

An integer vector, the pivot indices from igraph_lapack_dgetrf() must be given here. Row i of A was interchanged with row ipiv[i].

b:

The right hand side matrix must be given here. The solution will also be placed here.

Returns: 

Error code.

Time complexity: TODO.

2.1.3. igraph_lapack_dgesv — Solve system of linear equations with LU factorization

int igraph_lapack_dgesv(igraph_matrix_t *a, igraph_vector_int_t *ipiv,
                        igraph_matrix_t *b, int *info);

This function computes the solution to a real system of linear equations A * X = B, where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

The LU decomposition with partial pivoting and row interchanges is used to factor A as A = P * L * U, where P is a permutation matrix, L is unit lower triangular, and U is upper triangular. The factored form of A is then used to solve the system of equations A * X = B.

Arguments: 

a:

Matrix. On entry the N-by-N coefficient matrix, on exit, the factors L and U from the factorization A=P*L*U; the unit diagonal elements of L are not stored.

ipiv:

An integer vector or a null pointer. If not a null pointer, then the pivot indices that define the permutation matrix P, are stored here. Row i of the matrix was interchanged with row IPIV(i).

b:

Matrix, on entry the right hand side matrix should be stored here. On exit, if there was no error, and the info argument is zero, then it contains the solution matrix X.

info:

The LAPACK info code. If it is positive, then U(info,info) is exactly zero. In this case the factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.

Returns: 

Error code.

Time complexity: TODO.

Example 30.3.  File examples/simple/igraph_lapack_dgesv.c

/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2010-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
   02110-1301 USA

*/

#include <igraph.h>
#include <stdio.h>

#define DIM 10

void igraph_print_warning(const char *reason, const char *file,
                          int line, int igraph_errno) {
    printf("Warning: %s\n", reason);
}

int main() {

    igraph_matrix_t A, B, RHS;
    int info;
    int i, j;

    /* Identity matrix, you have to start somewhere */

    igraph_matrix_init(&A, DIM, DIM);
    igraph_matrix_init(&B, DIM, 1);
    for (i = 0; i < DIM; i++) {
        MATRIX(A, i, i) = 1.0;
        MATRIX(B, i, 0) = i + 1;
    }

    igraph_matrix_copy(&RHS, &B);
    igraph_lapack_dgesv(&A, /*ipiv=*/ 0, &RHS, &info);

    if (info != 0) {
        return 1;
    }
    if (!igraph_matrix_all_e(&B, &RHS)) {
        return 2;
    }

    igraph_matrix_destroy(&A);
    igraph_matrix_destroy(&B);
    igraph_matrix_destroy(&RHS);

    /* Diagonal matrix */

    igraph_matrix_init(&A, DIM, DIM);
    igraph_matrix_init(&RHS, DIM, 1);
    for (i = 0; i < DIM; i++) {
        MATRIX(A, i, i) = i + 1;
        MATRIX(RHS, i, 0) = i + 1;
    }

    igraph_lapack_dgesv(&A, /*ipiv=*/ 0, &RHS, &info);

    if (info != 0) {
        return 3;
    }
    for (i = 0; i < DIM; i++) {
        if (MATRIX(RHS, i, 0) != 1.0) {
            return 4;
        }
    }

    igraph_matrix_destroy(&A);
    igraph_matrix_destroy(&RHS);

    /* A general matrix */

    igraph_rng_seed(igraph_rng_default(), 42);

    igraph_matrix_init(&A, DIM, DIM);
    igraph_matrix_init(&B, DIM, 1);
    igraph_matrix_init(&RHS, DIM, 1);
    for (i = 0; i < DIM; i++) {
        int j;
        MATRIX(B, i, 0) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
        for (j = 0; j < DIM; j++) {
            MATRIX(A, i, j) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
        }
    }
    igraph_blas_dgemv_array(/*transpose=*/ 0, /*alpha=*/ 1.0, /*a=*/ &A,
                                           /*x-*/ &MATRIX(B, 0, 0), /*beta=*/ 0,
                                           /*y=*/ &MATRIX(RHS, 0, 0));

    igraph_lapack_dgesv(&A, /*ipiv=*/ 0, &RHS, &info);
    if (info != 0) {
        return 5;
    }
    for (i = 0; i < DIM; i++) {
        if (fabs(MATRIX(B, i, 0) - MATRIX(RHS, i, 0)) > 1e-13) {
            return 6;
        }
    }

    igraph_matrix_destroy(&A);
    igraph_matrix_destroy(&B);
    igraph_matrix_destroy(&RHS);

    /* A singular matrix */

    igraph_matrix_init(&A, DIM, DIM);
    igraph_matrix_init(&B, DIM, 1);
    igraph_matrix_init(&RHS, DIM, 1);
    for (i = 0; i < DIM; i++) {
        MATRIX(B, i, 0) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
        for (j = 0; j < DIM; j++) {
            MATRIX(A, i, j) = i == j ? 1 : 0;
        }
    }
    for (i = 0; i < DIM; i++) {
        MATRIX(A, DIM - 1, i) = MATRIX(A, 0, i);
    }

    igraph_blas_dgemv_array(/*transpose=*/ 0, /*alpha=*/ 1.0, /*a=*/ &A,
                                           /*x-*/ &MATRIX(B, 0, 0), /*beta=*/ 0,
                                           /*y=*/ &MATRIX(RHS, 0, 0));

    igraph_set_warning_handler(igraph_print_warning);
    igraph_lapack_dgesv(&A, /*ipiv=*/ 0, &RHS, &info);
    if (info != 10) {
        printf("LAPACK returned info = %d, should have been 10", info);
        return 7;
    }

    igraph_matrix_destroy(&A);
    igraph_matrix_destroy(&B);
    igraph_matrix_destroy(&RHS);

    return 0;
}


2.2. Eigenvalues and eigenvectors of matrices

2.2.1. igraph_lapack_dsyevr — Selected eigenvalues and optionally eigenvectors of a symmetric matrix
2.2.2. igraph_lapack_dgeev — Eigenvalues and optionally eigenvectors of a non-symmetric matrix
2.2.3. igraph_lapack_dgeevx — Eigenvalues/vectors of nonsymmetric matrices, expert mode

2.2.1. igraph_lapack_dsyevr — Selected eigenvalues and optionally eigenvectors of a symmetric matrix

int igraph_lapack_dsyevr(const igraph_matrix_t *A,
                         igraph_lapack_dsyev_which_t which,
                         igraph_real_t vl, igraph_real_t vu, int vestimate,
                         int il, int iu, igraph_real_t abstol,
                         igraph_vector_t *values, igraph_matrix_t *vectors,
                         igraph_vector_int_t *support);

Calls the DSYEVR LAPACK function to compute selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix A. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues.

See more in the LAPACK documentation.

Arguments: 

A:

Matrix, on entry it contains the symmetric input matrix. Only the leading N-by-N upper triangular part is used for the computation.

which:

Constant that gives which eigenvalues (and possibly the corresponding eigenvectors) to calculate. Possible values are IGRAPH_LAPACK_DSYEV_ALL, all eigenvalues; IGRAPH_LAPACK_DSYEV_INTERVAL, all eigenvalues in the half-open interval (vl,vu]; IGRAPH_LAPACK_DSYEV_SELECT, the il-th through iu-th eigenvalues.

vl:

If which is IGRAPH_LAPACK_DSYEV_INTERVAL, then this is the lower bound of the interval to be searched for eigenvalues. See also the vestimate argument.

vu:

If which is IGRAPH_LAPACK_DSYEV_INTERVAL, then this is the upper bound of the interval to be searched for eigenvalues. See also the vestimate argument.

vestimate:

An upper bound for the number of eigenvalues in the (vl,vu] interval, if which is IGRAPH_LAPACK_DSYEV_INTERVAL. Memory is allocated only for the given number of eigenvalues (and eigenvectors), so this upper bound must be correct.

il:

The index of the smallest eigenvalue to return, if which is IGRAPH_LAPACK_DSYEV_SELECT.

iu:

The index of the largets eigenvalue to return, if which is IGRAPH_LAPACK_DSYEV_SELECT.

abstol:

The absolute error tolerance for the eigevalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval [a,b] of width less than or equal to abstol + EPS * max(|a|,|b|), where EPS is the machine precision.

values:

An initialized vector, the eigenvalues are stored here, unless it is a null pointer. It will be resized as needed.

vectors:

An initialized matrix, the eigenvectors are stored in its columns, unless it is a null pointer. It will be resized as needed.

support:

An integer vector. If not a null pointer, then it will be resized to (2*max(1,M)) (M is a the total number of eigenvalues found). Then the support of the eigenvectors in vectors is stored here, i.e., the indices indicating the nonzero elements in vectors. The i-th eigenvector is nonzero only in elements support(2*i-1) through support(2*i).

Returns: 

Error code.

Time complexity: TODO.

Example 30.4.  File examples/simple/igraph_lapack_dsyevr.c

/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2010-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
   02110-1301 USA

*/

#include <igraph.h>

#define DIM 10

igraph_bool_t check_ev(const igraph_matrix_t *A,
                       const igraph_vector_t *values,
                       const igraph_matrix_t *vectors, igraph_real_t tol) {
    igraph_vector_t v, y;
    int i, j;
    int m = igraph_matrix_ncol(vectors);
    int n = igraph_matrix_nrow(A);

    if (igraph_matrix_ncol(A) != n)       {
        return 1;
    }
    if (igraph_vector_size(values) != m)  {
        return 1;
    }
    if (igraph_matrix_nrow(vectors) != n) {
        return 1;
    }

    igraph_vector_init(&y, n);

    for (i = 0; i < m; i++) {
        igraph_vector_view(&v, &MATRIX(*vectors, 0, i), n);
        igraph_vector_update(&y, &v);
        igraph_blas_dgemv(/*transpose=*/ 0, /*alpha=*/ 1.0, A, &v,
                                         /*beta=*/ -VECTOR(*values)[i], &y);
        for (j = 0; j < n; j++) {
            if (fabs(VECTOR(y)[i]) > tol) {
                printf("Matrix:\n");
                igraph_matrix_print(A);
                printf("lambda= %g\n", VECTOR(*values)[i]);
                printf("v= ");
                igraph_vector_print(&v);
                printf("residual: ");
                igraph_vector_print(&y);
                return 1;
            }
        }
    }

    igraph_vector_destroy(&y);
    return 0;
}

int main() {

    igraph_matrix_t A;
    igraph_matrix_t vectors, vectors2;
    igraph_vector_t values, values2;
    int i, j;
    int il, iu;
    igraph_real_t vl, vu;

    igraph_rng_seed(igraph_rng_default(), 42);

    igraph_matrix_init(&A, DIM, DIM);
    igraph_matrix_init(&vectors, 0, 0);
    igraph_vector_init(&values, 0);

    /* All eigenvalues and eigenvectors */

    for (i = 0; i < DIM; i++) {
        for (j = i; j < DIM; j++) {
            MATRIX(A, i, j) = MATRIX(A, j, i) =
                                  igraph_rng_get_integer(igraph_rng_default(), 1, 10);
        }
    }

    igraph_lapack_dsyevr(&A, IGRAPH_LAPACK_DSYEV_ALL, /*vl=*/ 0, /*vu=*/ 0,
                         /*vestimate=*/ 0, /*il=*/ 0, /*iu=*/ 0,
                         /*abstol=*/ 1e-10, &values, &vectors, /*support=*/ 0);

    if (igraph_vector_size(&values) != DIM) {
        return 1;
    }
    if (igraph_matrix_nrow(&vectors) != DIM ||
        igraph_matrix_ncol(&vectors) != DIM) {
        return 2;
    }
    if (check_ev(&A, &values, &vectors, /*tol=*/ 1e-8)) {
        return 3;
    }

    /* Only a subset */

    igraph_matrix_init(&vectors2, 0, 0);
    igraph_vector_init(&values2, 0);

    il = 2;
    iu = 5;
    igraph_lapack_dsyevr(&A, IGRAPH_LAPACK_DSYEV_SELECT, /*vl=*/ 0, /*vu=*/ 0,
                         /*vestimate=*/ 0, /*il=*/ il, /*iu=*/ iu,
                         /*abstol=*/ 1e-10, &values2, &vectors2,
                         /*support=*/ 0);

    if (igraph_vector_size(&values2) != iu - il + 1) {
        return 4;
    }
    if (igraph_matrix_nrow(&vectors2) != DIM ||
        igraph_matrix_ncol(&vectors2) != iu - il + 1) {
        return 5;
    }
    for (i = 0; i < iu - il + 1; i++) {
        igraph_real_t m1 = 1.0;

        if (fabs(VECTOR(values)[il + i - 1] - VECTOR(values2)[i]) > 1e-8) {
            printf("Full:   ");
            igraph_vector_print(&values);
            printf("Subset: ");
            igraph_vector_print(&values2);
            return 6;
        }

        if (MATRIX(vectors, 0, il + i - 1) * MATRIX(vectors2, 0, i) < 0) {
            m1 = -1.0;
        } else {
            m1 = 1.0;
        }

        for (j = 0; j < DIM; j++) {
            if (fabs(MATRIX(vectors, j, il + i - 1) -
                     m1 * MATRIX(vectors2, j, i)) > 1e-8) {
                printf("Full:\n");
                igraph_matrix_print(&vectors);
                printf("Subset:\n");
                igraph_matrix_print(&vectors2);
                return 7;
            }
        }
    }

    igraph_vector_destroy(&values2);
    igraph_matrix_destroy(&vectors2);

    /* Subset based on an interval */

    igraph_matrix_init(&vectors2, 0, 0);
    igraph_vector_init(&values2, 0);

    il = 2;
    iu = 5;
    vl = (VECTOR(values)[il - 1] + VECTOR(values)[il - 2]) / 2.0;
    vu = (VECTOR(values)[iu] + VECTOR(values)[iu - 1]) / 2.0;

    igraph_lapack_dsyevr(&A, IGRAPH_LAPACK_DSYEV_INTERVAL, vl, vu,
                         /*vestimate=*/ iu - il + 1, /*il=*/ 0, /*iu=*/ 0,
                         /*abstol=*/ 1e-10, &values2, &vectors2,
                         /*support=*/ 0);

    if (igraph_vector_size(&values2) != iu - il + 1) {
        return 4;
    }
    if (igraph_matrix_nrow(&vectors2) != DIM ||
        igraph_matrix_ncol(&vectors2) != iu - il + 1) {
        return 5;
    }
    for (i = 0; i < iu - il + 1; i++) {
        igraph_real_t m1 = 1.0;

        if (fabs(VECTOR(values)[il + i - 1] - VECTOR(values2)[i]) > 1e-8) {
            printf("Full:   ");
            igraph_vector_print(&values);
            printf("Subset: ");
            igraph_vector_print(&values2);
            return 6;
        }

        if (MATRIX(vectors, 0, il + i - 1) * MATRIX(vectors2, 0, i) < 0) {
            m1 = -1.0;
        } else {
            m1 = 1.0;
        }

        for (j = 0; j < DIM; j++) {
            if (fabs(MATRIX(vectors, j, il + i - 1) -
                     m1 * MATRIX(vectors2, j, i)) > 1e-8) {
                printf("Full:\n");
                igraph_matrix_print(&vectors);
                printf("Subset:\n");
                igraph_matrix_print(&vectors2);
                return 7;
            }
        }
    }

    igraph_vector_destroy(&values2);
    igraph_matrix_destroy(&vectors2);

    igraph_vector_destroy(&values);
    igraph_matrix_destroy(&vectors);
    igraph_matrix_destroy(&A);

    return 0;
}


2.2.2. igraph_lapack_dgeev — Eigenvalues and optionally eigenvectors of a non-symmetric matrix

int igraph_lapack_dgeev(const igraph_matrix_t *A,
                        igraph_vector_t *valuesreal,
                        igraph_vector_t *valuesimag,
                        igraph_matrix_t *vectorsleft,
                        igraph_matrix_t *vectorsright,
                        int *info);

This function calls LAPACK to compute, for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors.

The right eigenvector v(j) of A satisfies A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)**H * A = lambda(j) * u(j)**H where u(j)**H denotes the conjugate transpose of u(j).

The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real.

Arguments: 

A:

matrix. On entry it contains the N-by-N input matrix.

valuesreal:

Pointer to an initialized vector, or a null pointer. If not a null pointer, then the real parts of the eigenvalues are stored here. The vector will be resized as needed.

valuesimag:

Pointer to an initialized vector, or a null pointer. If not a null pointer, then the imaginary parts of the eigenvalues are stored here. The vector will be resized as needed.

vectorsleft:

Pointer to an initialized matrix, or a null pointer. If not a null pointer, then the left eigenvectors are stored in the columns of the matrix. The matrix will be resized as needed.

vectorsright:

Pointer to an initialized matrix, or a null pointer. If not a null pointer, then the right eigenvectors are stored in the columns of the matrix. The matrix will be resized as needed.

info:

This argument is used for two purposes. As an input argument it gives whether an igraph error should be generated if the QR algorithm fails to compute all eigenvalues. If info is non-zero, then an error is generated, otherwise only a warning is given. On exit it contains the LAPACK error code. Zero means successful exit. A negative values means that some of the arguments had an illegal value, this always triggers an igraph error. An i positive value means that the QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed; element i+1:N of valuesreal and valuesimag contain eigenvalues which have converged. This case only generates an igraph error, if info was non-zero on entry.

Returns: 

Error code.

Time complexity: TODO.

Example 30.5.  File examples/simple/igraph_lapack_dgeev.c

/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2010-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
   02110-1301 USA

*/

#include <igraph.h>
#include <stdio.h>

#define DIM 10

int real_cplx_mult(const igraph_matrix_t *A,
                   const igraph_vector_t *v_real,
                   const igraph_vector_t *v_imag,
                   igraph_vector_t *res_real,
                   igraph_vector_t *res_imag) {

    int n = igraph_vector_size(v_real);
    int r, c;

    if (igraph_matrix_nrow(A) != n ||
        igraph_matrix_ncol(A) != n ||
        igraph_vector_size(v_imag) != n) {
        printf("Wrong matrix or vector size");
        return 1;
    }

    igraph_vector_resize(res_real, n);
    igraph_vector_resize(res_imag, n);

    for (r = 0; r < n; r++) {
        igraph_real_t s_real = 0.0;
        igraph_real_t s_imag = 0.0;
        for (c = 0; c < n; c++) {
            s_real += MATRIX(*A, r, c) * VECTOR(*v_real)[c];
            s_imag += MATRIX(*A, r, c) * VECTOR(*v_imag)[c];
        }
        VECTOR(*res_real)[r] = s_real;
        VECTOR(*res_imag)[r] = s_imag;
    }

    return 0;
}

int sc_cplx_cplx_mult(igraph_real_t lambda_real,
                      igraph_real_t lambda_imag,
                      const igraph_vector_t *v_real,
                      const igraph_vector_t *v_imag,
                      igraph_vector_t *res_real,
                      igraph_vector_t *res_imag) {

    int r;
    int n = igraph_vector_size(v_real);

    if (igraph_vector_size(v_imag) != n) {
        printf("Wrong vector sizes");
        return 1;
    }

    igraph_vector_resize(res_real, n);
    igraph_vector_resize(res_imag, n);

    for (r = 0; r < n; r++) {
        VECTOR(*res_real)[r] = (lambda_real * VECTOR(*v_real)[r] -
                                lambda_imag * VECTOR(*v_imag)[r]);
        VECTOR(*res_imag)[r] = (lambda_imag * VECTOR(*v_real)[r] +
                                lambda_real * VECTOR(*v_imag)[r]);
    }

    return 0;
}

igraph_bool_t check_ev(const igraph_matrix_t *A,
                       const igraph_vector_t *values_real,
                       const igraph_vector_t *values_imag,
                       const igraph_matrix_t *vectors_left,
                       const igraph_matrix_t *vectors_right,
                       igraph_real_t tol) {

    int i, n = igraph_matrix_nrow(A);
    igraph_vector_t v_real, v_imag;
    igraph_vector_t AV_real, AV_imag, lv_real, lv_imag;
    igraph_vector_t null;

    if (igraph_matrix_ncol(A)             != n) {
        return 1;
    }
    if (igraph_vector_size(values_real)   != n) {
        return 1;
    }
    if (igraph_vector_size(values_imag)   != n) {
        return 1;
    }
    if (igraph_matrix_nrow(vectors_left)  != n) {
        return 1;
    }
    if (igraph_matrix_ncol(vectors_left)  != n) {
        return 1;
    }
    if (igraph_matrix_nrow(vectors_right) != n) {
        return 1;
    }
    if (igraph_matrix_ncol(vectors_right) != n) {
        return 1;
    }

    igraph_vector_init(&AV_real, n);
    igraph_vector_init(&AV_imag, n);
    igraph_vector_init(&lv_real, n);
    igraph_vector_init(&lv_imag, n);
    igraph_vector_init(&null, n);
    igraph_vector_null(&null);

    for (i = 0; i < n; i++) {
        if (VECTOR(*values_imag)[i] == 0.0) {
            igraph_vector_view(&v_real, &MATRIX(*vectors_right, 0, i), n);
            igraph_vector_view(&v_imag, VECTOR(null), n);
        } else if (VECTOR(*values_imag)[i] > 0.0) {
            igraph_vector_view(&v_real, &MATRIX(*vectors_right, 0, i), n);
            igraph_vector_view(&v_imag, &MATRIX(*vectors_right, 0, i + 1), n);
        } else if (VECTOR(*values_imag)[i] < 0.0) {
            igraph_vector_view(&v_real, &MATRIX(*vectors_right, 0, i - 1), n);
            igraph_vector_view(&v_imag, &MATRIX(*vectors_right, 0, i), n);
            igraph_vector_scale(&v_imag, -1.0);
        }
        real_cplx_mult(A, &v_real, &v_imag, &AV_real, &AV_imag);
        sc_cplx_cplx_mult(VECTOR(*values_real)[i], VECTOR(*values_imag)[i],
                          &v_real, &v_imag, &lv_real, &lv_imag);

        if (igraph_vector_maxdifference(&AV_real, &lv_real) > tol ||
            igraph_vector_maxdifference(&AV_imag, &lv_imag) > tol) {
            printf("ERROR:\n");
            igraph_vector_print(&AV_real);
            igraph_vector_print(&AV_imag);
            igraph_vector_print(&lv_real);
            igraph_vector_print(&lv_imag);
            return 1;
        }
    }

    igraph_vector_destroy(&null);
    igraph_vector_destroy(&AV_imag);
    igraph_vector_destroy(&AV_real);
    igraph_vector_destroy(&lv_imag);
    igraph_vector_destroy(&lv_real);

    return 0;
}

int main() {

    igraph_matrix_t A;
    igraph_matrix_t vectors_left, vectors_right;
    igraph_vector_t values_real, values_imag;
    int i, j;
    int info = 1;

    igraph_rng_seed(igraph_rng_default(), 42);

    igraph_matrix_init(&A, DIM, DIM);
    igraph_matrix_init(&vectors_left, 0, 0);
    igraph_matrix_init(&vectors_right, 0, 0);
    igraph_vector_init(&values_real, 0);
    igraph_vector_init(&values_imag, 0);

    for (i = 0; i < DIM; i++) {
        for (j = 0; j < DIM; j++) {
            MATRIX(A, i, j) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
        }
    }

    igraph_lapack_dgeev(&A, &values_real, &values_imag,
                        &vectors_left, &vectors_right, &info);

    if (check_ev(&A, &values_real, &values_imag,
                 &vectors_left, &vectors_right, /*tol=*/ 1e-8)) {
        return 1;
    }

    /* ------------------------------------------------------- */

    /* igraph_matrix_resize(&A, 10, 10); */
    /* igraph_matrix_null(&A); */
    /* for (i=0; i<10; i++) { MATRIX(A, i, i) = 1.0; }  */
    /* MATRIX(A,0,1) = 1.0; */

    /* igraph_lapack_dgeev(&A, &values_real, &values_imag,  */
    /*              &vectors_left, &vectors_right, &info);   */

    /* if (check_ev(&A, &values_real, &values_imag, */
    /*           &vectors_left, &vectors_right, /\*tol=*\/ 1e-8)) { */
    /*   return 2; */
    /* } */

    /* ------------------------------------------------------- */

    igraph_matrix_resize(&A, 10, 10);
    igraph_matrix_null(&A);
    MATRIX(A, 0, 1) = MATRIX(A, 0, 2) = MATRIX(A, 0, 3) = 1 / 3.0;
    MATRIX(A, 1, 0) = MATRIX(A, 1, 4) = MATRIX(A, 1, 5) = MATRIX(A, 1, 6) = 1 / 4.0;
    MATRIX(A, 2, 0) = MATRIX(A, 2, 7) = MATRIX(A, 2, 8) = MATRIX(A, 2, 9) = 1 / 4.0;
    MATRIX(A, 3, 0) = 1.0;
    MATRIX(A, 4, 1) = 1.0;
    MATRIX(A, 5, 1) = 1.0;
    MATRIX(A, 6, 1) = 1.0;
    MATRIX(A, 7, 2) = 1.0;
    MATRIX(A, 8, 2) = 1.0;
    MATRIX(A, 9, 2) = 1.0;

    info = 0;
    igraph_lapack_dgeev(&A, &values_real, &values_imag,
                        &vectors_left, &vectors_right, &info);

    /* igraph_matrix_print(&A); */
    /* printf("---\n"); */
    /* igraph_vector_print(&values_real); */
    /* igraph_vector_print(&values_imag); */
    /* igraph_matrix_print(&vectors_left); */

    if (check_ev(&A, &values_real, &values_imag,
                 &vectors_left, &vectors_right, /*tol=*/ 1e-8)) {
        return 3;
    }

    igraph_vector_destroy(&values_imag);
    igraph_vector_destroy(&values_real);
    igraph_matrix_destroy(&vectors_right);
    igraph_matrix_destroy(&vectors_left);
    igraph_matrix_destroy(&A);

    return 0;
}


2.2.3. igraph_lapack_dgeevx — Eigenvalues/vectors of nonsymmetric matrices, expert mode

int igraph_lapack_dgeevx(igraph_lapack_dgeevx_balance_t balance,
                         const igraph_matrix_t *A,
                         igraph_vector_t *valuesreal,
                         igraph_vector_t *valuesimag,
                         igraph_matrix_t *vectorsleft,
                         igraph_matrix_t *vectorsright,
                         int *ilo, int *ihi, igraph_vector_t *scale,
                         igraph_real_t *abnrm,
                         igraph_vector_t *rconde,
                         igraph_vector_t *rcondv,
                         int *info);

This function calculates the eigenvalues and optionally the left and/or right eigenvectors of a nonsymmetric N-by-N real matrix.

Optionally also, it computes a balancing transformation to improve the conditioning of the eigenvalues and eigenvectors (ilo, ihi, scale, and abnrm), reciprocal condition numbers for the eigenvalues (rconde), and reciprocal condition numbers for the right eigenvectors (rcondv).

The right eigenvector v(j) of A satisfies A * v(j) = lambda(j) * v(j) where lambda(j) is its eigenvalue. The left eigenvector u(j) of A satisfies u(j)^H * A = lambda(j) * u(j)^H where u(j)^H denotes the conjugate transpose of u(j).

The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real.

Balancing a matrix means permuting the rows and columns to make it more nearly upper triangular, and applying a diagonal similarity transformation D * A * D^(-1), where D is a diagonal matrix, to make its rows and columns closer in norm and the condition numbers of its eigenvalues and eigenvectors smaller. The computed reciprocal condition numbers correspond to the balanced matrix. Permuting rows and columns will not change the condition numbers (in exact arithmetic) but diagonal scaling will. For further explanation of balancing, see section 4.10.2 of the LAPACK Users' Guide.

Arguments: 

balance:

Scalar that indicated, whether the input matrix should be balanced. Possible values:

IGRAPH_LAPACK_DGEEVX_BALANCE_NONE

no not diagonally scale or permute.

IGRAPH_LAPACK_DGEEVX_BALANCE_PERM

perform permutations to make the matrix more nearly upper triangular. Do not diagonally scale.

IGRAPH_LAPACK_DGEEVX_BALANCE_SCALE

diagonally scale the matrix, i.e. replace A by D*A*D^(-1), where D is a diagonal matrix, chosen to make the rows and columns of A more equal in norm. Do not permute.

IGRAPH_LAPACK_DGEEVX_BALANCE_BOTH

both diagonally scale and permute A.

A:

The input matrix, must be square.

valuesreal:

An initialized vector, or a NULL pointer. If not a NULL pointer, then the real parts of the eigenvalues are stored here. The vector will be resized, as needed.

valuesimag:

An initialized vector, or a NULL pointer. If not a NULL pointer, then the imaginary parts of the eigenvalues are stored here. The vector will be resized, as needed.

vectorsleft:

An initialized matrix or a NULL pointer. If not a null pointer, then the left eigenvectors are stored here. The order corresponds to the eigenvalues and the eigenvectors are stored in a compressed form. If the j-th eigenvalue is real then column j contains the corresponding eigenvector. If the j-th and (j+1)-th eigenvalues form a complex conjugate pair, then the j-th and (j+1)-th columns contain their corresponding eigenvectors.

vectorsright:

An initialized matrix or a NULL pointer. If not a null pointer, then the right eigenvectors are stored here. The format is the same, as for the vectorsleft argument.

ilo:

ihi:

ilo and ihi are integer values determined when A was balanced. The balanced A(i,j) = 0 if I>J and J=1,...,ilo-1 or I=ihi+1,...,N.

scale:

Pointer to an initialized vector or a NULL pointer. If not a NULL pointer, then details of the permutations and scaling factors applied when balancing A, are stored here. If P(j) is the index of the row and column interchanged with row and column j, and D(j) is the scaling factor applied to row and column j, then

scale(J) = P(J), for J = 1,...,ilo-1

scale(J) = D(J), for J = ilo,...,ihi

scale(J) = P(J) for J = ihi+1,...,N.

The order in which the interchanges are made is N to ihi+1, then 1 to ilo-1.

abnrm:

Pointer to a real variable, the one-norm of the balanced matrix is stored here. (The one-norm is the maximum of the sum of absolute values of elements in any column.)

rconde:

An initialized vector or a NULL pointer. If not a null pointer, then the reciprocal condition numbers of the eigenvalues are stored here.

rcondv:

An initialized vector or a NULL pointer. If not a null pointer, then the reciprocal condition numbers of the right eigenvectors are stored here.

info:

This argument is used for two purposes. As an input argument it gives whether an igraph error should be generated if the QR algorithm fails to compute all eigenvalues. If info is non-zero, then an error is generated, otherwise only a warning is given. On exit it contains the LAPACK error code. Zero means successful exit. A negative values means that some of the arguments had an illegal value, this always triggers an igraph error. An i positive value means that the QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed; element i+1:N of valuesreal and valuesimag contain eigenvalues which have converged. This case only generated an igraph error, if info was non-zero on entry.

Returns: 

Error code.

Time complexity: TODO

Example 30.6.  File examples/simple/igraph_lapack_dgeevx.c

/* -*- mode: C -*-  */
/*
   IGraph library.
   Copyright (C) 2010-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
   02110-1301 USA

*/

#include <igraph.h>
#include <stdio.h>

#define DIM 10

int real_cplx_mult(const igraph_matrix_t *A,
                   const igraph_vector_t *v_real,
                   const igraph_vector_t *v_imag,
                   igraph_vector_t *res_real,
                   igraph_vector_t *res_imag) {

    int n = igraph_vector_size(v_real);
    int r, c;

    if (igraph_matrix_nrow(A) != n ||
        igraph_matrix_ncol(A) != n ||
        igraph_vector_size(v_imag) != n) {
        printf("Wrong matrix or vector size");
        return 1;
    }

    igraph_vector_resize(res_real, n);
    igraph_vector_resize(res_imag, n);

    for (r = 0; r < n; r++) {
        igraph_real_t s_real = 0.0;
        igraph_real_t s_imag = 0.0;
        for (c = 0; c < n; c++) {
            s_real += MATRIX(*A, r, c) * VECTOR(*v_real)[c];
            s_imag += MATRIX(*A, r, c) * VECTOR(*v_imag)[c];
        }
        VECTOR(*res_real)[r] = s_real;
        VECTOR(*res_imag)[r] = s_imag;
    }

    return 0;
}

int sc_cplx_cplx_mult(igraph_real_t lambda_real,
                      igraph_real_t lambda_imag,
                      const igraph_vector_t *v_real,
                      const igraph_vector_t *v_imag,
                      igraph_vector_t *res_real,
                      igraph_vector_t *res_imag) {

    int r;
    int n = igraph_vector_size(v_real);

    if (igraph_vector_size(v_imag) != n) {
        printf("Wrong vector sizes");
        return 1;
    }

    igraph_vector_resize(res_real, n);
    igraph_vector_resize(res_imag, n);

    for (r = 0; r < n; r++) {
        VECTOR(*res_real)[r] = (lambda_real * VECTOR(*v_real)[r] -
                                lambda_imag * VECTOR(*v_imag)[r]);
        VECTOR(*res_imag)[r] = (lambda_imag * VECTOR(*v_real)[r] +
                                lambda_real * VECTOR(*v_imag)[r]);
    }

    return 0;
}

igraph_bool_t check_ev(const igraph_matrix_t *A,
                       const igraph_vector_t *values_real,
                       const igraph_vector_t *values_imag,
                       const igraph_matrix_t *vectors_left,
                       const igraph_matrix_t *vectors_right,
                       igraph_real_t tol) {

    int n = igraph_matrix_nrow(A);
    igraph_vector_t v_real, v_imag;
    igraph_vector_t AV_real, AV_imag, lv_real, lv_imag;
    igraph_vector_t null;
    int i;

    if (igraph_matrix_ncol(A)             != n) {
        return 1;
    }
    if (igraph_vector_size(values_real)   != n) {
        return 1;
    }
    if (igraph_vector_size(values_imag)   != n) {
        return 1;
    }
    if (igraph_matrix_nrow(vectors_left)  != n) {
        return 1;
    }
    if (igraph_matrix_ncol(vectors_left)  != n) {
        return 1;
    }
    if (igraph_matrix_nrow(vectors_right) != n) {
        return 1;
    }
    if (igraph_matrix_ncol(vectors_right) != n) {
        return 1;
    }

    igraph_vector_init(&AV_real, n);
    igraph_vector_init(&AV_imag, n);
    igraph_vector_init(&lv_real, n);
    igraph_vector_init(&lv_imag, n);
    igraph_vector_init(&null, n);
    igraph_vector_null(&null);

    for (i = 0; i < n; i++) {
        if (VECTOR(*values_imag)[i] == 0.0) {
            igraph_vector_view(&v_real, &MATRIX(*vectors_right, 0, i), n);
            igraph_vector_view(&v_imag, VECTOR(null), n);
        } else if (VECTOR(*values_imag)[i] > 0.0) {
            igraph_vector_view(&v_real, &MATRIX(*vectors_right, 0, i), n);
            igraph_vector_view(&v_imag, &MATRIX(*vectors_right, 0, i + 1), n);
        } else if (VECTOR(*values_imag)[i] < 0.0) {
            igraph_vector_view(&v_real, &MATRIX(*vectors_right, 0, i - 1), n);
            igraph_vector_view(&v_imag, &MATRIX(*vectors_right, 0, i), n);
            igraph_vector_scale(&v_imag, -1.0);
        }
        real_cplx_mult(A, &v_real, &v_imag, &AV_real, &AV_imag);
        sc_cplx_cplx_mult(VECTOR(*values_real)[i], VECTOR(*values_imag)[i],
                          &v_real, &v_imag, &lv_real, &lv_imag);

        if (igraph_vector_maxdifference(&AV_real, &lv_real) > tol ||
            igraph_vector_maxdifference(&AV_imag, &lv_imag) > tol) {
            igraph_vector_print(&AV_real);
            igraph_vector_print(&AV_imag);
            igraph_vector_print(&lv_real);
            igraph_vector_print(&lv_imag);
            return 1;
        }
    }

    igraph_vector_destroy(&null);
    igraph_vector_destroy(&AV_imag);
    igraph_vector_destroy(&AV_real);
    igraph_vector_destroy(&lv_imag);
    igraph_vector_destroy(&lv_real);

    return 0;
}

int main() {

    igraph_matrix_t A;
    igraph_matrix_t vectors_left, vectors_right;
    igraph_vector_t values_real, values_imag;
    int i, j;
    int info = 1;
    int ilo, ihi;
    igraph_real_t abnrm;

    igraph_rng_seed(igraph_rng_default(), 42);

    igraph_matrix_init(&A, DIM, DIM);
    igraph_matrix_init(&vectors_left, 0, 0);
    igraph_matrix_init(&vectors_right, 0, 0);
    igraph_vector_init(&values_real, 0);
    igraph_vector_init(&values_imag, 0);

    for (i = 0; i < DIM; i++) {
        for (j = 0; j < DIM; j++) {
            MATRIX(A, i, j) = igraph_rng_get_integer(igraph_rng_default(), 1, 10);
        }
    }

    igraph_lapack_dgeevx(IGRAPH_LAPACK_DGEEVX_BALANCE_BOTH,
                         &A, &values_real, &values_imag,
                         &vectors_left, &vectors_right, &ilo, &ihi,
                         /*scale=*/ 0, &abnrm, /*rconde=*/ 0,
                         /*rcondv=*/ 0, &info);

    if (check_ev(&A, &values_real, &values_imag,
                 &vectors_left, &vectors_right, /*tol=*/ 1e-8)) {
        return 1;
    }

    /* igraph_matrix_print(&A); */
    /* igraph_vector_print(&values_real); */
    /* igraph_vector_print(&values_imag); */
    /* igraph_matrix_print(&vectors_left); */
    /* igraph_matrix_print(&vectors_right); */

    igraph_vector_destroy(&values_imag);
    igraph_vector_destroy(&values_real);
    igraph_matrix_destroy(&vectors_right);
    igraph_matrix_destroy(&vectors_left);
    igraph_matrix_destroy(&A);

    return 0;
}


3.  ARPACK interface in igraph

3.1. Data structures
3.2. ARPACK solvers

ARPACK is a library for solving large scale eigenvalue problems. The package is designed to compute a few eigenvalues and corresponding eigenvectors of a general n by n matrix A. It is most appropriate for large sparse or structured matrices A where structured means that a matrix-vector product w <- Av requires order n rather than the usual order n^2 floating point operations. Please see http://www.caam.rice.edu/software/ARPACK/ for details.

The eigenvalue calculation in ARPACK (in the simplest case) involves the calculation of the Av product where A is the matrix we work with and v is an arbitrary vector. A user-defined function of type igraph_arpack_function_t is expected to perform this product. If the product can be done efficiently, e.g. if the matrix is sparse, then ARPACK is usually able to calculate the eigenvalues very quickly.

In igraph, eigenvalue/eigenvector calculations usually involve the following steps:

  1. Initialization of an igraph_arpack_options_t data structure using igraph_arpack_options_init.

  2. Setting some options in the initialized igraph_arpack_options_t object.

  3. Defining a function of type igraph_arpack_function_t. The input of this function is a vector, and the output should be the output matrix multiplied by the input vector.

  4. Calling igraph_arpack_rssolve() (is the matrix is symmetric), or igraph_arpack_rnsolve().

The igraph_arpack_options_t object can be used multiple times.

If we have many eigenvalue problems to solve, then it might worth to create an igraph_arpack_storage_t object, and initialize it via igraph_arpack_storage_init(). This structure contains all memory needed for ARPACK (with the given upper limit regerding to the size of the eigenvalue problem). Then many problems can be solved using the same igraph_arpack_storage_t object, without always reallocating the required memory. The igraph_arpack_storage_t object needs to be destroyed by calling igraph_arpack_storage_destroy() on it, when it is not needed any more.

igraph does not contain all ARPACK routines, only the ones dealing with symmetric and non-symmetric eigenvalue problems using double precision real numbers.

3.1. Data structures

3.1.1. igraph_arpack_options_t — Options for ARPACK
3.1.2. igraph_arpack_storage_t — Storage for ARPACK
3.1.3. igraph_arpack_function_t — Type of the ARPACK callback function
3.1.4. igraph_arpack_options_init — Initialize ARPACK options
3.1.5. igraph_arpack_storage_init — Initialize ARPACK storage
3.1.6. igraph_arpack_storage_destroy — Deallocate ARPACK storage

3.1.1. igraph_arpack_options_t — Options for ARPACK

typedef struct igraph_arpack_options_t {
    /* INPUT */
    char bmat[1];          /* I-standard problem, G-generalized */
    int n;                 /* Dimension of the eigenproblem */
    char which[2];         /* LA, SA, LM, SM, BE */
    int nev;               /* Number of eigenvalues to be computed */
    igraph_real_t tol;     /* Stopping criterion */
    int ncv;               /* Number of columns in V */
    int ldv;               /* Leading dimension of V */
    int ishift;            /* 0-reverse comm., 1-exact with tridiagonal */
    int mxiter;            /* Maximum number of update iterations to take */
    int nb;                /* Block size on the recurrence, only 1 works */
    int mode;              /* The kind of problem to be solved (1-5)
                               1: A*x=l*x, A symmetric
                               2: A*x=l*M*x, A symm. M pos. def.
                               3: K*x = l*M*x, K symm., M pos. semidef.
                               4: K*x = l*KG*x, K s. pos. semidef. KG s. indef.
                               5: A*x = l*M*x, A symm., M symm. pos. semidef. */
    int start;             /* 0: random, 1: use the supplied vector */
    int lworkl;            /* Size of temporary storage, default is fine */
    igraph_real_t sigma;   /* The shift for modes 3,4,5 */
    igraph_real_t sigmai;  /* The imaginary part of shift for rnsolve */
    /* OUTPUT */
    int info;              /* What happened, see docs */
    int ierr;              /* What happened  in the dseupd call */
    int noiter;            /* The number of iterations taken */
    int nconv;
    int numop;             /* Number of OP*x operations */
    int numopb;            /* Number of B*x operations if BMAT='G' */
    int numreo;            /* Number of steps of re-orthogonalizations */
    /* INTERNAL */
    int iparam[11];
    int ipntr[14];
} igraph_arpack_options_t;

This data structure contains the options of thee ARPACK eigenvalue solver routines. It must be initialized by calling igraph_arpack_options_init() on it. Then it can be used for multiple ARPACK calls, as the ARPACK solvers do not modify it. Input options:

Values: 

bmat:

Character. Whether to solve a standard ('I') ot a generalized problem ('B').

n:

Dimension of the eigenproblem.

which:

Specifies which eigenvalues/vectors to compute. Possible values for symmetric matrices:

LA

Compute nev largest (algebraic) eigenvalues.

SA

Compute nev smallest (algebraic) eigenvalues.

LM

Compute nev largest (in magnitude) eigenvalues.

SM

Compute nev smallest (in magnitude) eigenvalues.

BE

Compute nev eigenvalues, half from each end of the spectrum. When nev is odd, compute one more from the high en than from the low end.

Possible values for non-symmetric matrices:

LM

Compute nev largest (in magnitude) eigenvalues.

SM

Compute nev smallest (in magnitude) eigenvalues.

LR

Compute nev eigenvalues of largest real part.

SR

Compute nev eigenvalues of smallest real part.

LI

Compute nev eigenvalues of largest imaginary part.

SI

Compute nev eigenvalues of smallest imaginary part.

nev:

The number of eigenvalues to be computed.

tol:

Stopping criterion: the relative accuracy of the Ritz value is considered acceptable if its error is less than tol times its estimated value. If this is set to zero then machine precision is used.

ncv:

Number of Lanczos vectors to be generated. Setting this to zero means that igraph_arpack_rssolve and igraph_arpack_rnsolve will determine a suitable value for ncv automatically.

ldv:

Numberic scalar. It should be set to zero in the current igraph implementation.

ishift:

Either zero or one. If zero then the shifts are provided by the user via reverse communication. If one then exact shifts with respect to the reduced tridiagonal matrix T. Please always set this to one.

mxiter:

Maximum number of Arnoldi update iterations allowed.

nb:

Blocksize to be used in the recurrence. Please always leave this on the default value, one.

mode:

The type of the eigenproblem to be solved. Possible values if the input matrix is symmetric:

  1. A*x=lambda*x, A is symmetric.

  2. A*x=lambda*M*x, A is symmetric, M is symmetric positive definite.

  3. K*x=lambda*M*x, K is symmetric, M is symmetric positive semi-definite.

  4. K*x=lambda*KG*x, K is symmetric positive semi-definite, KG is symmetric indefinite.

  5. A*x=lambda*M*x, A is symmetric, M is symmetric positive semi-definite. (Cayley transformed mode.)

Please note that only mode ==1 was tested and other values might not work properly. Possible values if the input matrix is not symmetric:

  1. A*x=lambda*x.

  2. A*x=lambda*M*x, M is symmetric positive definite.

  3. A*x=lambda*M*x, M is symmetric semi-definite.

  4. A*x=lambda*M*x, M is symmetric semi-definite.

Please note that only mode == 1 was tested and other values might not work properly.

start:

Whether to use the supplied starting vector (1), or use a random starting vector (0). The starting vector must be supplied in the first column of the vectors argument of the igraph_arpack_rssolve() of igraph_arpack_rnsolve() call.

Output options:

Values: 

info:

Error flag of ARPACK. Possible values:

0

Normal exit.

1

Maximum number of iterations taken.

3

No shifts could be applied during a cycle of the Implicitly restarted Arnoldi iteration. One possibility is to increase the size of ncv relative to nev.

ARPACK can return other error flags as well, but these are converted to igraph errors, see igraph_error_type_t.

ierr:

Error flag of the second ARPACK call (one eigenvalue computation usually involves two calls to ARPACK). This is always zero, as other error codes are converted to igraph errors.

noiter:

Number of Arnoldi iterations taken.

nconv:

Number of converged Ritz values. This represents the number of Ritz values that satisfy the convergence critetion.

numop:

Total number of matrix-vector multiplications.

numopb:

Not used currently.

numreo:

Total number of steps of re-orthogonalization.

Internal options:

Values: 

lworkl:

Do not modify this option.

sigma:

The shift for the shift-invert mode.

sigmai:

The imaginary part of the shift, for the non-symmetric or complex shift-invert mode.

iparam:

Do not modify this option.

ipntr:

Do not modify this option.

3.1.2. igraph_arpack_storage_t — Storage for ARPACK

typedef struct igraph_arpack_storage_t {
    int maxn, maxncv, maxldv;
    igraph_real_t *v;
    igraph_real_t *workl;
    igraph_real_t *workd;
    igraph_real_t *d;
    igraph_real_t *resid;
    igraph_real_t *ax;
    int *select;
    igraph_real_t *di;        /* These two only for non-symmetric problems */
    igraph_real_t *workev;
} igraph_arpack_storage_t;

Public members, do not modify them directly, these are considered to be read-only.

Values: 

maxn:

Maximum rank of matrix.

maxncv:

Maximum NCV.

maxldv:

Maximum LDV.

These members are considered to be private:

Values: 

workl:

Working memory.

workd:

Working memory.

d:

Memory for eigenvalues.

resid:

Memory for residuals.

ax:

Working memory.

select:

Working memory.

di:

Memory for eigenvalues, non-symmetric case only.

workev:

Working memory, non-symmetric case only.

3.1.3. igraph_arpack_function_t — Type of the ARPACK callback function

typedef int igraph_arpack_function_t(igraph_real_t *to, const igraph_real_t *from,
                                     int n, void *extra);

Arguments: 

to:

Pointer to an igraph_real_t, the result of the matrix-vector product is expected to be stored here.

from:

Pointer to an igraph_real_t, the input matrix should be multiplied by the vector stored here.

n:

The length of the vector (which is the same as the order of the input matrix).

extra:

Extra argument to the matrix-vector calculation function. This is coming from the igraph_arpack_rssolve() or igraph_arpack_rnsolve() function.

Returns: 

Error code, if not zero, then the ARPACK solver considers this as an error, stops and calls the igraph error handler.

3.1.4. igraph_arpack_options_init — Initialize ARPACK options

void igraph_arpack_options_init(igraph_arpack_options_t *o);

Initializes ARPACK options, set them to default values. You can always pass the initialized 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. 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.

Arguments: 

o:

The igraph_arpack_options_t object to initialize.

Time complexity: O(1).

3.1.5. igraph_arpack_storage_init — Initialize ARPACK storage

int igraph_arpack_storage_init(igraph_arpack_storage_t *s, long int maxn,
                               long int maxncv, long int maxldv,
                               igraph_bool_t symm);

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 storage argument of both igraph_arpack_rssolve() and igraph_arpack_rnsolve() to make memory allocated and deallocated automatically.

Don't forget to call the igraph_arpack_storage_destroy() function on the storage object if you don't need it any more.

Arguments: 

s:

The igraph_arpack_storage_t object to initialize.

maxn:

The maximum order of the matrices.

maxncv:

The maximum NCV parameter intended to use.

maxldv:

The maximum LDV parameter intended to use.

symm:

Whether symmetric or non-symmetric problems will be solved using this igraph_arpack_storage_t. (You cannot use the same storage both with symmetric and non-symmetric solvers.)

Returns: 

Error code.

Time complexity: O(maxncv*(maxldv+maxn)).

3.1.6. igraph_arpack_storage_destroy — Deallocate ARPACK storage

void igraph_arpack_storage_destroy(igraph_arpack_storage_t *s);

Arguments: 

s:

The igraph_arpack_storage_t object for which the memory will be deallocated.

Time complexity: operating system dependent.

3.2. ARPACK solvers

3.2.1. igraph_arpack_rssolve — ARPACK solver for symmetric matrices
3.2.2. igraph_arpack_rnsolve — ARPACK solver for non-symmetric matrices
3.2.3. igraph_arpack_unpack_complex — Make the result of the non-symmetric ARPACK solver more readable

3.2.1. igraph_arpack_rssolve — ARPACK solver for symmetric matrices

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);

This is the ARPACK solver for symmetric matrices. Please use igraph_arpack_rnsolve() for non-symmetric matrices.

Arguments: 

fun:

Pointer to an igraph_arpack_function_t object, the function that performs the matrix-vector multiplication.

extra:

An extra argument to be passed to fun.

options:

An igraph_arpack_options_t object.

storage:

An igraph_arpack_storage_t object, or a null pointer. In the latter case memory allocation and deallocation is performed automatically. Either this or the vectors argument must be non-null if the ARPACK iteration is started from a given starting vector. If both are given vectors take precedence.

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.

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 vectors argument must be non-null if the ARPACK iteration is started from a given starting vector. If both are given vectors take precedence.

Returns: 

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.

3.2.2. igraph_arpack_rnsolve — ARPACK solver for non-symmetric matrices

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);

Please always consider calling igraph_arpack_rssolve() if your matrix is symmetric, it is much faster. 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.

Arguments: 

fun:

Pointer to an igraph_arpack_function_t object, the function that performs the matrix-vector multiplication.

extra:

An extra argument to be passed to fun.

options:

An igraph_arpack_options_t object.

storage:

An igraph_arpack_storage_t object, or a null pointer. In the latter case memory allocation and deallocation is performed automatically.

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.

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 positive imaginary part only. Since in this case the eigenvector is also complex, it will occupy 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 negative imaginary part for a complex conjugate eigenvalue pair, the result vector will still store the eigenvector corresponding to the eigenvalue with the positive imaginary part (since this is how ARPACK works).

Returns: 

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.

3.2.3. igraph_arpack_unpack_complex — Make the result of the non-symmetric ARPACK solver more readable

int igraph_arpack_unpack_complex(igraph_matrix_t *vectors, igraph_matrix_t *values,
                                 long int nev);

This function works on the output of igraph_arpack_rnsolve and brushes it up a bit: it only keeps nev eigenvalues/vectors and every eigenvector is stored in two columns of the 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 igraph_arpack_rnsolve with its vectors and values argument and options->nev as nev. This will add the vectors for eigenvalues with a negative imaginary part and return all vectors as 2 columns, a real and imaginary part.

Arguments: 

vectors:

The eigenvector matrix, as returned by igraph_arpack_rnsolve. It will be resized, typically it will be larger.

values:

The eigenvalue matrix, as returned by igraph_arpack_rnsolve. It will be resized, typically extra, unneeded rows (=eigenvalues) will be removed.

nev:

The number of eigenvalues/vectors to keep. Can be less or equal than the number originally requested from ARPACK.

Returns: 

Error code.

Time complexity: linear in the number of elements in the vectors matrix.

← Chapter 29. Graph operators Chapter 31. Bipartite, i.e. two-mode graphs →