StarNEig Library  v0.1.3
A task-based library for solving dense nonsymmetric eigenvalue problems
sep_sm_full_chain.c
#include "validate.h"
#include <stdlib.h>
#include <stdio.h>
#include <time.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;
// allocate space for the eigenvalues and the eigenvalue 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 a default number of CPU cores and
// GPUs. The STARNEIG_HINT_SM flag indicates that the library should
// initialize itself for shared memory computations.
// reduce the full matrix matrix A to upper Hessenberg form
printf("Hessenberg reduction...\n");
starneig_SEP_SM_Hessenberg(n, A, ldA, Q, ldQ);
// reduce the upper Hessenberg matrix A to Schur form
printf("Schur reduction...\n");
starneig_SEP_SM_Schur(n, A, ldA, Q, ldQ, real, imag);
// select eigenvalues that have positive a real part
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);
// reorder the selected eigenvalues to the upper left corner of the matrix A
printf("Reordering...\n");
starneig_SEP_SM_ReorderSchur(n, select, A, ldA, Q, ldQ, real, imag);
// de-initialize the StarNEig library
// 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(real);
free(imag);
free(select);
return 0;
}