For using the igraph C library
igraph_blas_ddot — Dot product of two vectors.igraph_blas_dnrm2 — Euclidean norm of a vector.igraph_blas_dgemv — Matrix-vector multiplication using BLAS, vector version.igraph_blas_dgemm — Matrix-matrix multiplication using BLAS.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.
igraph_error_t igraph_blas_ddot(const igraph_vector_t *v1, const igraph_vector_t *v2,
igraph_real_t *res);
Arguments:
|
The first vector. |
|
The second vector. |
|
Pointer to a real, the result will be stored here. |
Time complexity: O(n) where n is the length of the vectors.
Example 29.1. File examples/simple/blas.c
#include <igraph.h> int main(void) { 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; }
igraph_real_t igraph_blas_dnrm2(const igraph_vector_t *v);
Arguments:
|
The vector. |
Returns:
Real value, the norm of |
Time complexity: O(n) where n is the length of the vector.
igraph_error_t 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:
|
Whether to transpose the matrix |
|
The constant |
|
The matrix |
|
The vector |
|
The constant |
|
The vector |
Time complexity: O(nk) if the matrix is of size n x k
Returns:
|
See also:
|
Example 29.2. File examples/simple/blas.c
#include <igraph.h> int main(void) { 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; }
igraph_error_t igraph_blas_dgemm(igraph_bool_t transpose_a, igraph_bool_t transpose_b,
igraph_real_t alpha, const igraph_matrix_t *a, const igraph_matrix_t *b,
igraph_real_t beta, igraph_matrix_t *c);
This function is a somewhat more user-friendly interface to
the dgemm function in BLAS. dgemm calculates
alpha*a*b + beta*c, where a, b and c are matrices, of which a and b
can be transposed.
Arguments:
|
whether to transpose the matrix |
|
whether to transpose the matrix |
|
the constant |
|
the matrix |
|
the matrix |
|
the constant |
|
the matrix |
Time complexity: O(n m k) where matrix a is of size n × k, and matrix b is of size k × m.
Returns:
|
Example 29.3. File examples/simple/blas_dgemm.c
#include <igraph.h> int main(void) { igraph_matrix_t a, b, c; igraph_matrix_init(&a, 2, 2); MATRIX(a, 0, 0) = 1; MATRIX(a, 0, 1) = 2; MATRIX(a, 1, 0) = 3; MATRIX(a, 1, 1) = 4; igraph_matrix_init(&b, 2, 2); MATRIX(b, 0, 0) = 5; MATRIX(b, 0, 1) = 6; MATRIX(b, 1, 0) = 7; MATRIX(b, 1, 1) = 8; igraph_matrix_init(&c, 2, 2); igraph_blas_dgemm(1, 1, 0.5, &a, &b, 0, &c); igraph_matrix_printf(&c, "%g"); igraph_matrix_destroy(&a); igraph_matrix_destroy(&b); igraph_matrix_destroy(&c); return 0; }
igraph_error_t 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:
|
whether to transpose the matrix |
|
the constant |
|
the matrix |
|
the vector |
|
the constant |
|
the vector |
Time complexity: O(nk) if the matrix is of size n x k
Returns:
|
See also:
|
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/
igraph_error_t 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:
|
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. |
|
An integer vector, the pivot indices are stored here,
unless it is a null pointer. Row |
|
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.
igraph_error_t 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:
|
Boolean, whether to transpose the input matrix. |
|
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. |
|
An integer vector, the pivot indices from
|
|
The right hand side matrix must be given here. The solution will also be placed here. |
Returns:
Error code. |
Time complexity: TODO.
igraph_error_t 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:
|
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. |
|
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). |
|
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. |
|
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 29.4. File examples/simple/igraph_lapack_dgesv.c
#include <igraph.h> #include <stdio.h> #define DIM 10 void igraph_print_warning(const char *reason, const char *file, int line) { IGRAPH_UNUSED(file); IGRAPH_UNUSED(line); printf("Warning: %s\n", reason); } int main(void) { 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_init_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-11) { 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; }
igraph_error_t 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:
|
Matrix, on entry it contains the symmetric input matrix. Only the leading N-by-N upper triangular part is used for the computation. |
|
Constant that gives which eigenvalues (and possibly
the corresponding eigenvectors) to calculate. Possible
values are |
|
If |
|
If |
|
An upper bound for the number of eigenvalues in
the (vl,vu] interval, if |
|
The index of the smallest eigenvalue to return, if |
|
The index of the largets eigenvalue to return, if |
|
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. |
|
An initialized vector, the eigenvalues are stored here, unless it is a null pointer. It will be resized as needed. |
|
An initialized matrix, the eigenvectors are stored in its columns, unless it is a null pointer. It will be resized as needed. |
|
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
|
Returns:
Error code. |
Time complexity: TODO.
Example 29.5. File examples/simple/igraph_lapack_dsyevr.c
#include <igraph.h> int main(void) { igraph_matrix_t A; igraph_matrix_t vectors; igraph_vector_t values; igraph_matrix_init(&A, 2, 2); igraph_matrix_init(&vectors, 0, 0); igraph_vector_init(&values, 0); MATRIX(A, 0, 0) = 2.0; MATRIX(A, 0, 1) = -1.0; MATRIX(A, 1, 0) = -1.0; MATRIX(A, 1, 1) = 3.0; printf("Take a subset:\n"); igraph_lapack_dsyevr(&A, IGRAPH_LAPACK_DSYEV_SELECT, /*vl=*/ 0, /*vu=*/ 0, /*vestimate=*/ 0, /*il=*/ 1, /*iu=*/ 1, /*abstol=*/ 1e-10, &values, &vectors, /*support=*/ 0); printf("eigenvalues:\n"); igraph_vector_print(&values); printf("eigenvectors:\n"); igraph_matrix_print(&vectors); printf("\nTake a subset based on an interval:\n"); igraph_lapack_dsyevr(&A, IGRAPH_LAPACK_DSYEV_INTERVAL, /*vl*/ 3, /*vu*/ 4, /*vestimate=*/ 1, /*il=*/ 0, /*iu=*/ 0, /*abstol=*/ 1e-10, &values, &vectors, /*support=*/ 0); printf("eigenvalues:\n"); igraph_vector_print(&values); printf("eigenvectors:\n"); igraph_matrix_print(&vectors); igraph_vector_destroy(&values); igraph_matrix_destroy(&vectors); igraph_matrix_destroy(&A); return 0; }
igraph_error_t 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:
|
matrix. On entry it contains the N-by-N input matrix. |
|
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. |
|
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. |
|
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. |
|
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. |
|
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 |
Returns:
Error code. |
Time complexity: TODO.
Example 29.6. File examples/simple/igraph_lapack_dgeev.c
#include <igraph.h> #include <stdio.h> int main(void) { igraph_matrix_t A; igraph_matrix_t vectors_left, vectors_right; igraph_vector_t values_real, values_imag; igraph_vector_complex_t values; int info = 1; igraph_matrix_init(&A, 2, 2); 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); MATRIX(A, 0, 0) = 1.0; MATRIX(A, 0, 1) = 1.0; MATRIX(A, 1, 0) = -1.0; MATRIX(A, 1, 1) = 1.0; igraph_lapack_dgeev(&A, &values_real, &values_imag, &vectors_left, &vectors_right, &info); igraph_vector_complex_create(&values, &values_real, &values_imag); printf("eigenvalues:\n"); igraph_vector_complex_print(&values); printf("left eigenvectors:\n"); igraph_matrix_print(&vectors_left); printf("right eigenvectors:\n"); igraph_matrix_print(&vectors_right); igraph_vector_destroy(&values_imag); igraph_vector_destroy(&values_real); igraph_vector_complex_destroy(&values); igraph_matrix_destroy(&vectors_right); igraph_matrix_destroy(&vectors_left); igraph_matrix_destroy(&A); return 0; }
igraph_error_t 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. Note that the eigenvectors obtained for the balanced
matrix are backtransformed to those of A.
Arguments:
|
Indicates whether the input matrix should be balanced. Possible values:
|
||||||||
|
The input matrix, must be square. |
||||||||
|
An initialized vector, or a |
||||||||
|
An initialized vector, or a |
||||||||
|
An initialized matrix or a |
||||||||
|
An initialized matrix or a |
||||||||
|
|
||||||||
|
if not NULL, |
||||||||
|
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
The order in which the interchanges are made is N to |
||||||||
|
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.) |
||||||||
|
An initialized vector or a NULL pointer. If not a null pointer, then the reciprocal condition numbers of the eigenvalues are stored here. |
||||||||
|
An initialized vector or a NULL pointer. If not a null pointer, then the reciprocal condition numbers of the right eigenvectors are stored here. |
||||||||
|
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 |
Returns:
Error code. |
Time complexity: TODO
Example 29.7. File examples/simple/igraph_lapack_dgeevx.c
#include <igraph.h> #include <stdio.h> void matrix_to_complex_vectors(igraph_vector_complex_t *c1, igraph_vector_complex_t *c2, igraph_matrix_t *m) { IGRAPH_REAL(VECTOR(*c1)[0]) = MATRIX(*m, 0, 0); IGRAPH_REAL(VECTOR(*c1)[1]) = MATRIX(*m, 1, 0); IGRAPH_IMAG(VECTOR(*c1)[0]) = MATRIX(*m, 0, 1); IGRAPH_IMAG(VECTOR(*c1)[1]) = MATRIX(*m, 1, 1); IGRAPH_REAL(VECTOR(*c2)[0]) = MATRIX(*m, 0, 0); IGRAPH_REAL(VECTOR(*c2)[1]) = MATRIX(*m, 1, 0); IGRAPH_IMAG(VECTOR(*c2)[0]) = -MATRIX(*m, 0, 1); IGRAPH_IMAG(VECTOR(*c2)[1]) = -MATRIX(*m, 1, 1); } int main(void) { igraph_matrix_t A; igraph_matrix_t vectors_left, vectors_right; igraph_vector_t values_real, values_imag; igraph_vector_complex_t values; igraph_vector_complex_t eigenvector1; igraph_vector_complex_t eigenvector2; int info = 1; igraph_real_t abnrm; igraph_matrix_init(&A, 2, 2); 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); igraph_vector_complex_init(&eigenvector1, 2); igraph_vector_complex_init(&eigenvector2, 2); MATRIX(A, 0, 0) = 1.0; MATRIX(A, 0, 1) = 1.0; MATRIX(A, 1, 0) = -1.0; MATRIX(A, 1, 1) = 1.0; igraph_lapack_dgeevx(IGRAPH_LAPACK_DGEEVX_BALANCE_BOTH, &A, &values_real, &values_imag, &vectors_left, &vectors_right, NULL, NULL, /*scale=*/ NULL, &abnrm, /*rconde=*/ NULL, /*rcondv=*/ NULL, &info); igraph_vector_complex_create(&values, &values_real, &values_imag); printf("eigenvalues:\n"); igraph_vector_complex_print(&values); printf("\nleft eigenvectors:\n"); /*matrix_to_complex_vectors only works because we have two complex conjugate eigenvalues */ matrix_to_complex_vectors(&eigenvector1, &eigenvector2, &vectors_left); igraph_vector_complex_print(&eigenvector1); igraph_vector_complex_print(&eigenvector2); printf("\nright eigenvectors:\n"); matrix_to_complex_vectors(&eigenvector1, &eigenvector2, &vectors_right); igraph_vector_complex_print(&eigenvector1); igraph_vector_complex_print(&eigenvector2); printf("\nOne-norm of the balanced matrix:\n%g\n", abnrm); igraph_vector_destroy(&values_imag); igraph_vector_destroy(&values_real); igraph_vector_complex_destroy(&values); igraph_vector_complex_destroy(&eigenvector1); igraph_vector_complex_destroy(&eigenvector2); igraph_matrix_destroy(&vectors_right); igraph_matrix_destroy(&vectors_left); igraph_matrix_destroy(&A); return 0; }
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:
Initialization of an igraph_arpack_options_t data
structure using igraph_arpack_options_init.
Setting some options in the initialized igraph_arpack_options_t object.
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.
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.
igraph_arpack_options_t — Options for ARPACK.igraph_arpack_storage_t — Storage for ARPACK.igraph_arpack_function_t — Type of the ARPACK callback function.igraph_arpack_options_init — Initialize ARPACK options.igraph_arpack_storage_init — Initialize ARPACK storage.igraph_arpack_storage_destroy — Deallocate ARPACK storage.
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 the 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:
|
Character. Whether to solve a standard ('I') ot a generalized problem ('B'). |
||||||||||||||||||||||
|
Dimension of the eigenproblem. |
||||||||||||||||||||||
|
Specifies which eigenvalues/vectors to compute. Possible values for symmetric matrices:
Possible values for non-symmetric matrices:
|
||||||||||||||||||||||
|
The number of eigenvalues to be computed. |
||||||||||||||||||||||
|
Stopping criterion: the relative accuracy
of the Ritz value is considered acceptable if its error is less
than |
||||||||||||||||||||||
|
Number of Lanczos vectors to be generated. Setting this
to zero means that |
||||||||||||||||||||||
|
Numberic scalar. It should be set to zero in the current igraph implementation. |
||||||||||||||||||||||
|
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 |
||||||||||||||||||||||
|
Maximum number of Arnoldi update iterations allowed. |
||||||||||||||||||||||
|
Blocksize to be used in the recurrence. Please always leave this on the default value, one. |
||||||||||||||||||||||
|
The type of the eigenproblem to be solved. Possible values if the input matrix is symmetric:
Please note that only
Please note that only |
||||||||||||||||||||||
|
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 |
Output options:
Values:
|
Error flag of ARPACK. Possible values:
ARPACK can return other error flags as well, but these are
converted to igraph errors, see |
||||||
|
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. |
||||||
|
Number of Arnoldi iterations taken. |
||||||
|
Number of converged Ritz values. This represents the number of Ritz values that satisfy the convergence critetion. |
||||||
|
Total number of matrix-vector multiplications. |
||||||
|
Not used currently. |
||||||
|
Total number of steps of re-orthogonalization. |
Internal options:
Values:
|
Do not modify this option. |
|
The shift for the shift-invert mode. |
|
The imaginary part of the shift, for the non-symmetric or complex shift-invert mode. |
|
Do not modify this option. |
|
Do not modify this option. |
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;
/* The following two are only used for non-symmetric problems: */
igraph_real_t *di;
igraph_real_t *workev;
} igraph_arpack_storage_t;
Public members, do not modify them directly, these are considered to be read-only.
Values:
|
Maximum rank of matrix. |
|
Maximum NCV. |
|
Maximum LDV. |
These members are considered to be private:
Values:
|
Working memory. |
|
Working memory. |
|
Memory for eigenvalues. |
|
Memory for residuals. |
|
Working memory. |
|
Working memory. |
|
Memory for eigenvalues, non-symmetric case only. |
|
Working memory, non-symmetric case only. |
typedef igraph_error_t igraph_arpack_function_t(igraph_real_t *to, const igraph_real_t *from,
int n, void *extra);
Arguments:
|
Pointer to an |
|
Pointer to an |
|
The length of the vector (which is the same as the order of the input matrix). |
|
Extra argument to the matrix-vector calculation
function. This is coming from the |
Returns:
Error code. If not |
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:
|
The |
Time complexity: O(1).
igraph_error_t igraph_arpack_storage_init(igraph_arpack_storage_t *s, igraph_integer_t maxn,
igraph_integer_t maxncv, igraph_integer_t 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:
|
The |
|
The maximum order of the matrices. |
|
The maximum NCV parameter intended to use. |
|
The maximum LDV parameter intended to use. |
|
Whether symmetric or non-symmetric problems will be
solved using this |
Returns:
Error code. |
Time complexity: O(maxncv*(maxldv+maxn)).
void igraph_arpack_storage_destroy(igraph_arpack_storage_t *s);
Arguments:
|
The |
Time complexity: operating system dependent.
igraph_error_t 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:
|
Pointer to an |
|
An extra argument to be passed to |
|
An |
|
An |
|
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. |
|
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 |
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.
igraph_error_t 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:
|
Pointer to an |
|
An extra argument to be passed to |
|
An |
|
An |
|
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. |
|
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.
igraph_error_t igraph_arpack_unpack_complex(igraph_matrix_t *vectors, igraph_matrix_t *values,
igraph_integer_t 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:
|
The eigenvector matrix, as returned by |
|
The eigenvalue matrix, as returned by |
|
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 28. Graph operators | Chapter 30. Bipartite, i.e. two-mode graphs → |