避免 LAPACK 中的矩阵半矢量化

Question

我的问题的答案很可能是 "No"，但也许有人对此问题有巧妙的解决方案？

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

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

Answer 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() returns 右特征向量 ，而不是左特征向量。左特征向量是右特征向量的复共轭，如 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() 为工作数组分配了一些 space。

避免 LAPACK 中的矩阵半矢量化

Avoid matrix half-vectorization in LAPACK

c

c++

optimization

matrix

lapack