【问题标题】:Avoid matrix half-vectorization in LAPACK避免 LAPACK 中的矩阵半向量化
【发布时间】:2017-02-15 04:31:45
【问题描述】:

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

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

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

【问题讨论】:

  • 您是否有迹象表明这是一个问题? zheev 或 LAPACK 中任何有趣的东西的运行时间可能会支配转换所需的时间。
  • @AmiTavory 你可能是对的,但是,仅仅因为 lapack 有这个奇怪的接口,就不得不对我的矩阵进行矢量化和去矢量化,这很烦人。此外,内存成本也不容忽视。我处理边长约为 2^16 的矩阵。所以我想避免浪费内存。
  • 在您链接的示例中,zheev() 使用 LDA.N=N.N 矩阵。实际上,矩阵不需要打包:您可能不必复制矩阵。注意:zheev() 修改 A!您可能对zhpev()(专用于打包存储)感兴趣,或者您可能对zhbev()(专用于带存储,如果它适合您的矩阵)感兴趣。
  • @francis 非常感谢!实际上,我在示例中注意到了这一点,但无法解释。我们可以说每个函数都有一个包含非密集元素的副本吗?
  • 不客气! LAPACK 的每个功能都涉及到一个特定的矩阵存储。请参阅 LAPACK 的the naming convention。在ZHEEV 中,Z 指的是双精度复数,HE 指的是内存中大小为 N.N 的 Hermitian 矩阵。根据参数UPLO,使用或忽略上三角部分。无论如何,矩阵可以被填充,就好像它是一个大小为N.N 的一般解包矩阵。在这种情况下,参数 UPLO 的值不应改变获得的特征值。最后EV指的是特征值的计算。

标签: c++ c optimization matrix lapack


【解决方案1】:

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

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

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

当您使用 C 或 C++ 时,这里是通过接口 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() 为工作数组分配了一些空间。

【讨论】:

    猜你喜欢
    • 2013-05-30
    • 1970-01-01
    • 1970-01-01
    • 2013-02-19
    • 1970-01-01
    • 2014-06-25
    • 2023-01-11
    • 2019-04-25
    • 2013-01-14
    相关资源
    最近更新 更多