Using the cblas_chpr function

I'm having issue when I tried to use the cblas_chpr() function from the cblas library to calculate the correlation matrix of a float complex vector.

After downloading the Lapack v3.10.0 Library from netLib.org, I compiled it and copied the libcblas.a, liblapack.a, liblapacke.a, librefblas.a and libtmglib.a files to my project and make sure that the libraries are linked correctly.

According to the description, the cblas_chpr function calculates alpha * x * conjg(x') + A and stores the result in A.

The function is defined as:

void cblas_chpr(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
                const CBLAS_INDEX N, const float alpha, const void *X,
                const CBLAS_INDEX incX, void *A);

Where the parameters are:

layout - this is an emun and the two possible inputs are CblasRowMajor and CblasColMajor.
Uplo - this is an enum and the two possible inputs are CblasUpper and CblasLower.
N - The order of matrix A and the number of elements in vector x.
alpha - The scaling factor that vector X is multiplied by.
X - Vector X.
incX - Stride within X. For example, if incX is 7, every 7th element is used.
A - Matrix A. Overwritten by results on return.

The body of my function is as followed:

   /* Number of elements */
   int Ne = 10;


   /* Set the parameters */
   CBLAS_LAYOUT layout = CblasColMajor;   /* Layout is column major */
   CBLAS_UPLO Uplo = CblasUpper;          /* Upper triangle of the matrix */
   CBLAS_INDEX N = Ne;                    /* Number of elements in vector X */
   float alpha = 1.0;                     /* No scaling, alpha = 1.0 */

   /* The vector X */
   float complex * X = malloc(Ne * sizeof(* X));

   /* Set values of X - for illustration purpose only */
   for(int i = 0; i < Ne; i++)
   {
      X[i] = CMPLXF(i, i + 1.0);
   }

   CBLAS_INDEX incX = 1;                  /* Use data from every element */

   /* The correlation matrix is a Ne x Ne matrix */
   float complex ** A = malloc(Ne * sizeof(*A));

   for(int i = 0; i < Ne; i++)
   {
      A[i] = malloc(Ne * sizeof(*A[i]));
   }

   cblas_chpr(layout, Uplo, N, alpha, X, incX, A);

   float complex print_val = A[0][0];
   printf("%+.10f %+.10f", crealf(print_val), cimagf(print_val));

However, the program crashed with the No source available for "chpr_() at 0x55555555e70b" error.

I'm guessing that my input parameters are not correct. CBLAS is a wrapper for the Fortran BLAS Library.

Has anyone encountered this error before and know how to resolve it?

Solution

Answering my own question in case someone else ran into the same issue. A should be a 1D array of size N * (N + 1) / 2. Also, the value each element in array A must be initialize to zero. Otherwise, the result will be wrong. Read the description for the cblas_chpr() function for why this is the case.

/* Number of elements */
int Ne = 10;

/* Set the parameters */
CBLAS_LAYOUT layout = CblasColMajor;   /* Layout is column major */
CBLAS_UPLO Uplo = CblasUpper;          /* Upper triangle of the matrix */
CBLAS_INDEX N = Ne;                    /* Number of elements in vector X */
float alpha = 1.0;                     /* No scaling, alpha = 1.0 */

/* The vector X */
float complex * X = malloc(Ne * sizeof(* X));

/* Set values of X - for illustration purpose only */
for(int i = 0; i < Ne; i++)
{
  X[i] = CMPLXF(i, i + 1.0);
}

CBLAS_INDEX incX = 1;                  /* Use data from every element */

/* Initialize the array that store correlation matrix */
int size_A = Ne * (Ne + 1) / 2;
float complex * A = malloc(size_A * sizeof(*A));

for(int i = 0; i < size_A; i++)
{
    A[i] = 0.0;
}

cblas_chpr(layout, Uplo, N, alpha, X, incX, A);

/* Print the first value of the result */
float complex print_val = A[0][0];
printf("%+.10f %+.10f", crealf(print_val), cimagf(print_val));

free(X);
free(A);