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 operationssuch as vector-vector, matrix-vector and matrix-matrix product.Please seehttp://www.netlib.org/blas/ for details and a referenceimplementation in Fortran. igraph contains some wrapper functionsthat can be used to call BLAS routines in a somewhat moreuser-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 drivenby igraph's internal requirements. The wrapper functions usuallysubstitute double-precision floating point arrays used by BLAS withigraph_vector_t andigraph_matrix_t instances and alsoremove those parameters (such as the number of rows/columns) thatcan 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. |
Returns:
Error code. |
Time complexity: O(n) where n is the length of the vectors.
Example 29.1. Fileexamples/simple/blas.c
#include <igraph.h>intmain(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 tothedgemv function in BLAS.dgemv performs the operationy = alpha*A*x + beta*y, wherex andy are vectorsandA 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. Fileexamples/simple/blas.c
#include <igraph.h>intmain(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 tothedgemm function in BLAS.dgemm calculatesalpha*a*b + beta*c, where a, b and c are matrices, of which a and bcan 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 ofsize k × m.
Returns:
|
Example 29.3. Fileexamples/simple/blas_dgemm.c
#include <igraph.h>intmain(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 tothedgemv function in BLAS.dgemv performs the operationy = alpha*A*x + beta*y, where x and y are vectors and A is anappropriately 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 solvingsystems of simultaneous linear equations, least-squares solutionsof linear systems of equations, eigenvalue problems, and singularvalue problems. The associated matrix factorizations (LU, Cholesky,QR, SVD, Schur, generalized Schur) are also provided, as arerelated computations such as reordering of the Schur factorizationsand estimating condition numbers. Dense and banded matrices arehandled, but not general sparse matrices. In all areas, similarfunctionality is provided for real and complex matrices, in bothsingle and double precision.
igraph provides an interface to a very limited set of LAPACKfunctions, using the regular igraph data structures.
See more about LAPACK athttp://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 * Uwhere P is a permutation matrix, L is lower triangular with unitdiagonal elements (lower trapezoidal if m > n), and U is uppertriangular (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 = Bwith a general N-by-N matrix A using the LU factorizationcomputed byigraph_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 linearequations A * X = B, where A is an N-by-N matrix and X and B areN-by-NRHS matrices.
The LU decomposition with partial pivoting and rowinterchanges is used to factor A as A = P * L * U,where P is a permutation matrix, L is unit lower triangular, and U isupper triangular. The factored form of A is then used to solve thesystem 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. Fileexamples/simple/igraph_lapack_dgesv.c
#include <igraph.h>#include <stdio.h>#define DIM 10voidigraph_print_warning(const char *reason,const char *file, int line) {IGRAPH_UNUSED(file);IGRAPH_UNUSED(line);printf("Warning: %s\n", reason);}intmain(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 eigenvaluesand, optionally, eigenvectors of a real symmetric matrixA.Eigenvalues and eigenvectors can be selected by specifying eithera 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 |
| 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 |
| An initialized vector, the eigenvalues are stored here, unless it is a null pointer. It will be resized as needed. |
| An initialized matrix. A set of orthonormal 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. Fileexamples/simple/igraph_lapack_dsyevr.c
#include <igraph.h>intmain(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 realnonsymmetric matrix A, the eigenvalues and, optionally, the leftand/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)^Hwhere u(j)^H denotes the conjugate transpose of u(j).
The computed eigenvectors are normalized to have Euclidean normequal 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. Fileexamples/simple/igraph_lapack_dgeev.c
#include <igraph.h>#include <stdio.h>intmain(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 leftand/or right eigenvectors of a nonsymmetric N-by-N real matrix.
Optionally also, it computes a balancing transformation to improvethe conditioning of the eigenvalues and eigenvectors (ilo,ihi,scale, andabnrm), reciprocal condition numbers for theeigenvalues (rconde), and reciprocal condition numbers for theright 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)^Hwhere u(j)^H denotes the conjugate transpose of u(j).
The computed eigenvectors are normalized to have Euclidean normequal to 1 and largest component real.
Balancing a matrix means permuting the rows and columns to make itmore nearly upper triangular, and applying a diagonal similaritytransformation D * A * D^(-1), where D is a diagonal matrix, tomake its rows and columns closer in norm and the condition numbersof its eigenvalues and eigenvectors smaller. The computedreciprocal 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 furtherexplanation of balancing, see section 4.10.2 of the LAPACKUsers' Guide. Note that the eigenvectors obtained for the balancedmatrix are backtransformed to those ofA.
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. Fileexamples/simple/igraph_lapack_dgeevx.c
#include <igraph.h>#include <stdio.h>voidmatrix_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);}intmain(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 correspondingeigenvectors of a generaln byn matrixA. It ismost appropriate for large sparse or structured matricesA wherestructured means that a matrix-vector productw <- Av requiresordern rather than the usual ordern^2 floating pointoperations. Please seehttps://github.com/opencollab/arpack-ng for details.
The eigenvalue calculation in ARPACK (in the simplestcase) involves the calculation of theAv product whereAis the matrix we work with andv is an arbitrary vector. Auser-defined function of typeigraph_arpack_function_tis expected to perform this product. If the product can be doneefficiently, e.g. if the matrix is sparse, then ARPACK is usuallyable to calculate the eigenvalues very quickly.
In igraph, eigenvalue/eigenvector calculations usuallyinvolve the following steps:
Initialization of anigraph_arpack_options_t data structure usingigraph_arpack_options_init.
Setting some options in the initializedigraph_arpack_options_t object.
Defining a function of typeigraph_arpack_function_t. The input of this function is a vector, and the output should be the output matrix multiplied by the input vector.
Callingigraph_arpack_rssolve() (is the matrix is symmetric), origraph_arpack_rnsolve().
Theigraph_arpack_options_t object can be used multipletimes.
If we have many eigenvalue problems to solve, then it might worthto create anigraph_arpack_storage_t object, and initialize itviaigraph_arpack_storage_init(). This structure contains allmemory needed for ARPACK (with the given upper limit regerding tothe size of the eigenvalue problem). Then many problems can besolved using the sameigraph_arpack_storage_t object, withoutalways reallocating the required memory.Theigraph_arpack_storage_t object needs to be destroyed bycallingigraph_arpack_storage_destroy() on it, when it is notneeded any more.
igraph does not contain allARPACK routines, only the ones dealing with symmetric andnon-symmetric eigenvalue problems using double precision realnumbers.
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 eigenvaluesolver routines. It must be initialized by callingigraph_arpack_options_init() on it. Then it can be used formultiple 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 consideredto 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 initializedigraph_arpack_options_tobject to built-in igraph functions without any modification. Thebuilt-in igraph functions modify the options to perform theircalculation, e.g.igraph_pagerank() always searches for theeigenvalue with the largest magnitude, regardless of the suppliedvalue.
If you want to implement your own function involving eigenvaluecalculation using ARPACK, however, you will likely need to set upthe 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 eigenvaluecalculations using ARPACK, and want to spare the memoryallocation/deallocation between each two runs. Otherwise it is safeto supply a null pointer as thestorage argument of bothigraph_arpack_rssolve() andigraph_arpack_rnsolve() to makememory allocated and deallocated automatically.
Don't forget to call theigraph_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 useigraph_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-vectormultiplication. Usually a small number of iterations is enough, soif the matrix is sparse and the matrix-vector multiplication can bedone in O(n) time (the number of vertices), then the eigenvaluesare 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 callingigraph_arpack_rssolve() if yourmatrix 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 algebraicsolution 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 withpositive imaginary part only. Since in this case the eigenvector is also complex, it will occupytwo columns in the eigenvector matrix (the real and the imaginary parts, in this order). Caveat: if the eigenvalue vector returns only the eigenvalue with thenegative imaginary part for a complex conjugate eigenvalue pair, the result vector willstill 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-vectormultiplication. Usually a small number of iterations is enough, soif the matrix is sparse and the matrix-vector multiplication can bedone in O(n) time (the number of vertices), then the eigenvaluesare 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 ofigraph_arpack_rnsolve andbrushes it up a bit: it only keepsnev eigenvalues/vectors andevery eigenvector is stored in two columns of thevectorsmatrix.
The output of the non-symmetric ARPACK solver is somewhat hard toparse, 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 moreeigenvalues than requested. The common use of this function is tocall it directly afterigraph_arpack_rnsolve with itsvectors andvalues argument andoptions->nev asnev.This will add the vectors for eigenvalues with a negative imaginarypart 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 thevectorsmatrix.
| ← Chapter 28. Graph operators | Chapter 30. Bipartite, i.e. two-mode graphs → |