gep_sm_full_chain.c

#include "validate.h"
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <starneig/starneig.h>

// a predicate function that selects all finate eigenvalues that have positive
// a real part
static int predicate(double real, double imag, double beta, void *arg)
{
    if (0.0 < real && beta != 0.0)
        return 1;
    return 0;
}

int main()
{
    const int n = 3000; // matrix dimension

    srand((unsigned) time(NULL));

    // generate a full random matrix A and a copy C

    int ldA = ((n/8)+1)*8, ldC = ((n/8)+1)*8;
    double *A = malloc(n*ldA*sizeof(double));
    double *C = malloc(n*ldC*sizeof(double));
    for (int j = 0; j < n; j++)
        for (int i = 0; i < n; i++)
            A[j*ldA+i] = C[j*ldC+i] = 2.0*rand()/RAND_MAX - 1.0;

    // generate a full random matrix B and a copy D

    int ldB = ((n/8)+1)*8, ldD = ((n/8)+1)*8;
    double *B = malloc(n*ldB*sizeof(double));
    double *D = malloc(n*ldD*sizeof(double));
    for (int j = 0; j < n; j++)
        for (int i = 0; i < n; i++)
            B[j*ldB+i] = D[j*ldD+i] = 2.0*rand()/RAND_MAX - 1.0;

    // generate an identity matrix Q

    int ldQ = ((n/8)+1)*8;
    double *Q = malloc(n*ldA*sizeof(double));
    for (int j = 0; j < n; j++)
        for (int i = 0; i < n; i++)
            Q[j*ldQ+i] = i == j ? 1.0 : 0.0;

    // generate an identity matrix Z

    int ldZ = ((n/8)+1)*8;
    double *Z = malloc(n*ldZ*sizeof(double));
    for (int j = 0; j < n; j++)
        for (int i = 0; i < n; i++)
            Z[j*ldZ+i] = i == j ? 1.0 : 0.0;

    // allocate space for the eigenvalues and the eigenvalue selection vector

    double *real = malloc(n*sizeof(double));
    double *imag = malloc(n*sizeof(double));
    double *beta = malloc(n*sizeof(double));
    int *select = malloc(n*sizeof(int));

    // Initialize the StarNEig library using a default number of CPU cores and
    // GPUs. The STARNEIG_HINT_SM flag indicates that the library should
    // initialize itself for shared memory computations and the
    // STARNEIG_AWAKE_WORKERS indicates that the library should should keep
    // StarPU worker threads awake between interface function calls.

    starneig_node_init(-1, -1, STARNEIG_HINT_SM | STARNEIG_AWAKE_WORKERS);

    // reduce the dense-dense matrix pencil (A,B) to Hessenberg-triangular form

    printf("Hessenberg-triangular reduction...\n");
    starneig_GEP_SM_HessenbergTriangular(n, A, ldA, B, ldB, Q, ldQ, Z, ldZ);

    // reduce the Hessenberg-triangular matrix pencil (A,B) to generalized Schur
    // form

    printf("Schur reduction...\n");
    starneig_GEP_SM_Schur(n, A, ldA, B, ldB, Q, ldQ, Z, ldZ, real, imag, beta);

    // select eigenvalues that have positive a real part

    int num_selected;
    starneig_GEP_SM_Select(
        n, A, ldA, B, ldB, &predicate, NULL, select, &num_selected);
    printf("Selected %d eigenvalues out of %d.\n", num_selected, n);

    // reorder selected eigenvalues to the upper left corner of the generalized
    // Schur form (A,B)

    printf("Reordering...\n");
    starneig_GEP_SM_ReorderSchur(
        n, select, A, ldA, B, ldB, Q, ldQ, Z, ldZ, real, imag, beta);

    // de-initialize the StarNEig library

    starneig_node_finalize();

    // check residual || Q A Z^T - C ||_F / || C ||_F

    check_residual(n, ldQ, ldA, ldZ, ldC, Q, A, Z, C);

    // check residual || Q B Z^T - D ||_F / || D ||_F

    check_residual(n, ldQ, ldB, ldZ, ldD, Q, B, Z, D);

    // check residual || Q Q^T - I ||_F / || I ||_F

    check_orthogonality(n, ldQ, Q);

    // check residual || Z Z^T - I ||_F / || I ||_F

    check_orthogonality(n, ldZ, Z);

    // cleanup

    free(A);
    free(C);
    free(B);
    free(D);
    free(Q);
    free(Z);

    free(real);
    free(imag);
    free(beta);
    free(select);

    return 0;
}