c++ - 避免 LAPACK 中的矩阵半矢量化

Question

我的问题的答案很可能是“否”，但也许有人对这个问题有聪明的解决方案？

问题来了。例如，lapack 函数 zheev计算复数 Hermitian 矩阵的特征值。问题是矩阵的所有 C++ 实现都存储行主矩阵或列主矩阵，而 zheev() 采用密集的上三角矩阵或下三角矩阵。

所以我的问题是:有什么方法可以避免将我的矩阵复制到仅存储下三角部分或上三角部分的新数组并在 lapack 中使用我当前的完整矩阵？

Answer 1

您在 zheev() 上链接的示例已经使用了未压缩的 LDA*N=N*N 矩阵。事实上，厄密矩阵不需要打包:您可能不必复制矩阵。注意:zheev() 修改矩阵 A!

LAPACK 处理矩阵的其他存储模式:参见 the naming scheme拉帕克。例如:

• zheev() : 内存占用 N*N 和存储类似于一般解压缩的 N*N 矩阵之一。根据参数 UPLO 的值，使用或忽略上三角部分。无论如何，矩阵可以像大小为 N*N 的一般解压缩矩阵一样填充。在这种情况下，参数 UPLO 的值不应改变获得的特征值。
• zhpev() : 打包格式。根据参数 UPLO 的值，存储上对角线项或下对角线项。矩阵存储的内存占用量为 (N*(N+1))/2。
• zhbev() :专用于乐队存储。

当您使用 C 或 C++ 时，这里是通过接口(interface) LAPACKE 使用 zheev() 的示例代码。它可以通过gcc main.c -o main -llapacke -llapack -lblas -lm -Wall编译。此外，此代码确保函数 zheev() 返回右特征向量，而不是左特征向量。左特征向量是右特征向量的复共轭，如解释的那样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;
}

在上面的代码中，执行了矩阵的拷贝，但这不是 LAPACKE_zheev() 所要求的。处理一个2000*2000的矩阵，上面代码的内存占用约为167MB。这是矩阵大小 (64MB) 的两倍多，因为执行了复制。但如果不进行复制的话，也不会超过两次。因此，LAPACKE_zheev() 不执行矩阵的任何复制。另请注意 LAPACKE_zheev()为工作数组分配一些空间。