sep_sm_eigenvectors.c

 
#include "validate.h"
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <starneig/starneig.h>
 
// a predicate function that selects all eigenvalues that have positive a real
// part
static int predicate(double real, double imag, void *arg)
{
    if (0.0 < real)
        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 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;
 
    double *X = NULL; int ldX = 0;
 
    // allocate space for the eigenvalues and the eigenvector selection vector
 
    double *real = malloc(n*sizeof(double));
    double *imag = malloc(n*sizeof(double));
    int *select = malloc(n*sizeof(int));
 
    // Initialize the StarNEig library using all available CPU cores and
    // GPUs. The STARNEIG_HINT_SM flag indicates that the library should
    // initialize itself for shared memory computations.
 
    starneig_node_init(STARNEIG_USE_ALL, STARNEIG_USE_ALL, STARNEIG_HINT_SM);
 
    // reduce the full matrix matrix A to Schur form (skip reordering)
 
    printf("Reduce...\n");
    starneig_SEP_SM_Reduce(
        n, A, ldA, Q, ldQ, real, imag, NULL, NULL, NULL, NULL);
 
    // select eigenvalues that have positive a real part and allocate space for
    // the eigenvectors
 
    int num_selected;
    starneig_SEP_SM_Select(n, A, ldA, &predicate, NULL, select, &num_selected);
    printf("Selected %d eigenvalues out of %d.\n", num_selected, n);
 
    ldX = ((n/8)+1)*8;
    X = malloc(num_selected*ldX*sizeof(double));
 
    // compute a selected set of eigenvectors
 
    printf("Eigenvectors...\n");
    starneig_SEP_SM_Eigenvectors(n, select, A, ldA, Q, ldQ, X, ldX);
 
    // de-initialize the StarNEig library
 
    starneig_node_finalize();
 
    // check residual || Q A Q^T - C ||_F / || C ||_F
 
    check_residual(n, ldQ, ldA, ldQ, ldC, Q, A, Q, C);
 
    // check residual || Q Q^T - I ||_F / || I ||_F
 
    check_orthogonality(n, ldQ, Q);
 
    // cleanup
 
    free(A);
    free(C);
    free(Q);
    free(X);
 
    free(real);
    free(imag);
    free(select);
 
    return 0;
}