Avoid matrix half-vectorization in LAPACK

Question

The answer to my question is most likely "No", but maybe someone has a smart solution to this problem?

Here's the problem. For example, the lapack function zheev calculates the eigenvalues of a complex Hermitian matrix. The problem is that all C++ implementations of matrices store either row-major or column-major matrices, while zheev() takes a dense upper or lower triangular matrix.

So my question is: Is there any way to avoid copying my matrix to a new array that stores only the lower or upper triangular part and use my current full matrix in lapack?

Answer 1

The example you linked on zheev() already makes use of an unpacked LDA*N=N*N matrix. Indeed, the hermitian matrix does not need to be packed: you may not have to make a copy of your matrix. Watch out: zheev() modifies the matrix A!

LAPACK handles other storage mode for matrices: see the naming scheme of LAPACK. For instance:

• zheev(): the memory footprint N*N and the storage are similar to the one of general unpacked N*N matrices. Depending on the value of the argument UPLO, the upper triangular part is used or ignored. Anyway, the matrix can be populated as if it were a general unpacked matrix of size N*N. In this case, the value of the argument UPLO should not change the obtained eigenvalues.
• zhpev(): packed format. Either the upper diagonal items or the lower diagonal items are stored, depending on the value of the argument UPLO. The memory footprint for matrix storage is (N*(N+1))/2.
• zhbev(): dedicated to band storage.

As you work with C or C++, here is a sample code using zheev() via the interface LAPACKE. It can be compiled by gcc main.c -o main -llapacke -llapack -lblas -lm -Wall. Moreover, this code ensures that the function zheev() returns the right eigenvectors, not the left ones. The left eigenvectors are the the complex conjugates of the right ones, as explained here.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <lapacke.h>


int main(void){

    int n=200;

    srand(time(NULL));

    // allocate the matrix, unpacked storage
    double complex** A=malloc(n*sizeof(double complex*));
    if(A==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    A[0]=malloc(n*n*sizeof(double complex));
    if(A[0]==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    int i;
    for(i=1;i<n;i++){
        A[i]=&A[0][n*i];
    }

    //populte the matrix, only the lower diagonal part
    int j;
    for(i=0;i<n;i++){
        for(j=0;j<i;j++){
            A[i][j]=rand()/((double)RAND_MAX)+rand()/((double)RAND_MAX)*I;

        }
        A[i][i]=rand()/((double)RAND_MAX)+42.;
    }



    //saving the matrix, because zheev() changes it.
    double complex** As=malloc(n*sizeof(double complex*));
    if(As==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    As[0]=malloc(n*n*sizeof(double complex));
    if(As[0]==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    for(i=1;i<n;i++){
        As[i]=&As[0][n*i];
    }
    for(i=0;i<n;i++){
        for(j=0;j<i;j++){
            As[i][j]=A[i][j];
        }
        As[i][i]=A[i][i];
    }
    //transpose part, conjugate
    for(i=0;i<n;i++){
        for(j=i+1;j<n;j++){
            As[i][j]=conj(A[j][i]);
        }
    }
    /*
for(i=0;i<n;i++){
for(j=0;j<n;j++){
printf("%g+I%g ",creal(As[i][j]),cimag(As[i][j]));
}
printf("\n");
}
     */

    //let's get the eigenvalues and the eigenvectors!
    double* w=malloc(n*sizeof(double));
    if(w==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    int lda = n;
    int info = LAPACKE_zheev(LAPACK_ROW_MAJOR, 'V', 'L', n, A[0], lda,  w);
    if(info<0){
        fprintf(stderr,"argument %d has illegal value\n",-info);
    }
    if(info>0){
        fprintf(stderr,"algorithm failed to converge... bad condition number ?\n");
    }

    //printing the eigenvalues...
    for(i=0;i<n;i++){
        printf("%d %g\n",i,w[i]);
    }

    //let's check that the column i of A is now a right eigenvectors corresponding to the eigenvalue w[i]...
    int l=42;

    double complex *res=malloc(n*sizeof(double complex));
    if(res==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    for(i=0;i<n;i++){
        res[i]=0;
        for(j=0;j<n;j++){
            res[i]+=As[i][j]*A[j][l];
        }
        res[i]-=w[l]*A[i][l];
    }
    double norm2res=0;
    double norm2vec=0;
    for(i=0;i<n;i++){
        norm2res+=creal(res[i])*creal(res[i])+cimag(res[i])*cimag(res[i]);
        norm2vec+=creal(A[i][l])*creal(A[i][l])+cimag(A[i][l])*cimag(A[i][l]);
    }
    printf("the norm of the eigenvector is %g\n",sqrt(norm2vec));
    printf("||Ax-\\lambda*x||_2/||x||_2 is about %g\n",sqrt(norm2res/norm2vec));
    //free the matrix
    free(A[0]);
    free(A);

    free(As[0]);
    free(As);

    free(w);
    free(res);
    return 0;
}

In the code above, a copy of the matrix is performed, but that is not required by LAPACKE_zheev(). Dealing with a matrix of 2000*2000, the memory footprint of the code above is about 167MB. That's more than twice of the size of the matrix (64MB) because a copy is performed. But it would be less than twice if the copy were not performed. Hence, LAPACKE_zheev() does not perform any copy of the matrix. Notice also that LAPACKE_zheev() allocates some space for the working array.