?geev

Interface Definition

C interface:

void sgeev_( const char *jobvl, const char *jobvr, const int *n, float *a, const int *lda, float *wr, float *wi, float *vl, const int *ldvl, float *vr, const int *ldvr, float *work, int *lwork, int *info);

void dgeev_(const char *jobvl, const char *jobvr, const int *n, double *a, const int *lda, double *wr, double *wi, double *vl, const int *ldvl, double *vr, const int *ldvr, double *work, int *lwork, int *info);

void cgeev_( const char *jobvl, const char *jobvr, const int *n, float _Complex *a, const int *lda, float _Complex *w, float _Complex *vl, const int *ldvl, float _Complex *vr, const int *ldvr, float _Complex *work, int *lwork, float *rwork, int *info);

void zgeev_( const char *jobvl, const char *jobvr, const int *n, double _Complex *a, const int *lda, double _Complex *w, double _Complex *vl, const int *ldvl, double _Complex *vr, const int *ldvr, double _Complex *work, int *lwork, double *rwork, int *info);

Fortran interface:

SGEEV(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info);

DGEEV(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info);

CGEEV(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info);

ZGEEV(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info);

Parameters

Dependencies

#include "klapack.h"

Examples

C interface:

#include <stdio.h>
#include <stdlib.h>
#include "lapack.h"

int main() {
    int i, j;
    int n = 3; // Matrix dimensions.
    int lda = n;
    int ldvl = n;
    int ldvr = n;
    int lwork;
    int info;

    // Define matrix A (3 x 3).
    float a[9] = {
        1.0, 2.0, 3.0,
        4.0, 5.0, 6.0,
        7.0, 8.0, 9.0
    };

    // Output the original matrix.
    printf("Original matrix A:\n");
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            printf("%8.4f ", a[i*n + j]);
        }
        printf("\n");
    }
    printf("\n");

    // Allocate memory for storing eigenvalues and eigenvectors.
    float *wr = (float*)malloc(n * sizeof(float));  // Real part.
    float *wi = (float*)malloc(n * sizeof(float));  // Imaginary part.
    float *vl = (float*)malloc(n * n * sizeof(float));  // Left eigenvector.
    float *vr = (float*)malloc(n * n * sizeof(float));  // Right eigenvector.

    // Query the optimal workspace size.
    lwork = -1;
    float query_work;
    sgeev_("V", "V", &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, &query_work, &lwork, &info,1,1);

    // Allocate a workspace.
    lwork = (int)query_work;
    float *work = (float*)malloc(lwork * sizeof(float));

    // Compute the eigenvalues and eigenvectors.
    sgeev_("V", "V", &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info,1,1);

    // Check whether the operation is successful.
    if (info != 0) {
        printf("SGEEV failed with error code %d\n", info);
        return 1;
    }

    // Output the eigenvalues.
    printf("Eigenvalues:\n");
    for (i = 0; i < n; i++) {
        if (wi[i] == 0.0) {
            printf("%8.4f\n", wr[i]);
        } else {
            printf("%8.4f + %8.4fi\n", wr[i], wi[i]);
            printf("%8.4f - %8.4fi\n", wr[i], -wi[i]);
            i++;  // Skip the conjugate pairs.
        }
    }
    printf("\n");

    // Output the right eigenvectors.
    printf("Right eigenvectors:\n");
    for (i = 0; i < n; i++) {
        if (wi[i] == 0.0) {
            printf("Eigenvector for eigenvalue %8.4f:\n", wr[i]);
            for (j = 0; j < n; j++) {
                printf("%8.4f ", vr[j*n + i]);
            }
            printf("\n");
        } else {
            printf("Eigenvector pair for eigenvalue %8.4f + %8.4fi:\n", wr[i], wi[i]);
            for (j = 0; j < n; j++) {
                printf("%8.4f + %8.4fi ", vr[j*n + i], vr[j*n + i+1]);
            }
            printf("\n");
            for (j = 0; j < n; j++) {
                printf("%8.4f - %8.4fi ", vr[j*n + i], -vr[j*n + i+1]);
            }
            printf("\n");
            i++;  // Skip the conjugate pairs.
        }
    }

    // Free memory.
    free(wr);
    free(wi);
    free(vl);
    free(vr);
    free(work);

    return 0;
}

Fortran interface:

     ! Define constants and variables.
    integer, parameter :: n = 3          ! Matrix dimension.
    character(1) :: jobvl = 'N'         ! Do not compute the left eigenvectors.
    character(1) :: jobvr = 'V'         ! Compute the right eigenvectors.
    integer :: lda = n, ldvl = n, ldvr = n
    integer :: lwork, info, i, j,k
    double precision :: a(n, n)         ! Input matrix (column-major order).
    double precision :: wr(n), wi(n)    ! Real and imaginary parts of the eigenvalues.
    double precision :: vr(n, n), vl(n, n)        ! Eigenvectors (column-major order).
    double precision, allocatable :: work(:)  ! Workspace array.

    ! Initialize the test matrix (column-major order).
    ! Matrix example:
    ! [ 1.0  4.0  7.0 ]
    ! [ 2.0  5.0  8.0 ]
    ! [ 3.0  6.0  9.0 ]
    a = reshape([1.0d0, 2.0d0, 3.0d0, &
                 4.0d0, 5.0d0, 6.0d0, &
                 7.0d0, 8.0d0, 9.0d0], [n, n])

    ! === Step 1: Query the optimal workspace size. ===
    lwork = -1
    allocate(work(1))  ! Temporarily allocate the minimum space.
    call dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)

    ! Check whether the query is successful.
    if (info /= 0) then
        print *, "Workspace query failed. Error code: ", info
        stop
    end if

    ! Allocate the workspace based on the query result.
    lwork = int(work(1))
    deallocate(work)
    allocate(work(lwork))

    ! === Step 2: Compute the eigenvalues and eigenvectors. ===
    call dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)

    ! Check whether the computation is successful.
    if (info /= 0) then
        print *, "Eigenvalue computation failed. Error code: ", info
        stop
    end if

    k = 1
    ! === Output results ===
    print *, "=== Eigenvalues ==="
    do i = 1, n
        if (wi(k) == 0.0d0) then
            ! Real eigenvalues.
            write(*, '(A, I1, A, F12.6)') "λ", k, " = ", wr(k)
        else
            ! Complex conjugate pairs (one pair displayed).
            write(*, '(A, I1, A, F12.6, A, F12.6, A)') "λ", k, " = ", wr(k), " ± ", wi(k), "i"
            k = k + 1  ! Skip the next conjugate imaginary part.
        end if
        k = k + 1
    end do

    print *, ""
    print *, "=== Right eigenvectors (column-major order) ==="
    do j = 1, n
        write(*, '(A, I1, A)', advance='no') "Vector", j, ": "
        do i = 1, n
            write(*, '(F12.6, 2X)', advance='no') vr(i, j)
        end do
        print *, ""
    end do

    ! Free memory.
    deallocate(work)