1

我正在计算密集非对称矩阵 A 的特征值。为此，我使用xGEHRD和xHSEQR Lapack 例程来首先计算 A 的上 Hessenberg 形式，然后计算所获得矩阵的唯一特征值。

这两个例程都需要参数 LWORK，并且都提供了一种机制来计算其最佳值。我相信这个参数与缓冲技术的内部阻塞有关，但我不知道它是如何确定的。

使用查询机制获取 LWORK 的最优值，工作流程应该是这样的：

int LWORK = -1;
float* OPT_LWORK = (float*) malloc( sizeof(float));

sgehrd_ (..., OPT_LWORK ,&LWORK, ...) // query optimal value for sgehrd

LWORK = (int) OPT_WORK
float* WORK = (float*) malloc( (int) sizeof(float) * LWORK );

sgehrd_ (..., WORK ,&LWORK, ...) // calculate Hessenberg

int LWORK = -1;
shseqr_ (..., OPT_LWORK ,&LWORK, ...) // query optimal value for shseqr

LWORK = (int) OPT_WORK
float* WORK = // possibly realloc with the new LWORK value

shseqr_ (..., WORK ,&LWORK, ...) // finally obtain eigenvalues

我已经进行了一些测试，始终为 WORK 数组的维度获得相同的最佳值。如果值相同，我可以大大简化我的代码（不需要 realloc 并且只需要一次调用来确定 LWORK 的值，更少的错误检查......）。

我的问题是，对于相同的矩阵和相同的 ILO 和 IHI 值，我可以假设这两个例程的值将相等吗？

4

1 回答 1

1

查看sgehrd.f，似乎WORKfor 例程的最佳大小sgehrd()是N*NB，其中NB是

NB = MIN( NBMAX, ILAENV( 1, 'SGEHRD', ' ', N, ILO, IHI, -1 ) )

哪里NBMAX=64。因此，最优LWORK取决于N和。ILOIHI

查看shseqr.f，最佳长度的计算LWORK更复杂：slaqr0()调用例程......但文件中的文档指出：

如果 LWORK = -1，则 SHSEQR 执行工作区查询。在这种情况下，SHSEQR 检查输入参数并估计给定 N、ILO 和 IHI 值的最佳工作空间大小。估算值在 WORK(1) 中返回。XERBLA 不会发出与 LWORK 相关的错误消息。既不访问 H 也不访问 Z。

的最佳长度对于和WORK可能不同。这是一个例子：sgehrd()shseqr()

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

extern void sgehrd_(int *n,int* ilo, int* ihi, float* a,int* lda,float* tau,float* work, int* lwork,int* info);

extern void shseqr_(char* job,char* compz,int *n,int* ilo, int* ihi,float* h,int* ldh,float* wr,float* wi,float* z,int* ldz,float* work, int* lwork,int* info);

int main()
{
    int n=10;
    int ilo=n;
    int ihi=n;
    float* a=malloc(sizeof(float)*n*n);
    int lda=n;
    float* tau=malloc(sizeof(float)*(n-1));
    int info;

    char job='S';
    char compz='I';
    float* h=malloc(sizeof(float)*n*n);
    int ldh=n;
    float* wr=malloc(sizeof(float)*(n));
    float* wi=malloc(sizeof(float)*(n));
    float* z=malloc(sizeof(float)*n*n);
    int ldz=n;

    int LWORK = -1;
    float* OPT_LWORK =(float*) malloc( sizeof(float));

    sgehrd_ (&n,&ilo, &ihi, a,&lda,tau,OPT_LWORK ,&LWORK,&info); // query optimal value for sgehrd
    if(info!=0){printf("sgehrd lwork=-1 : failed\n");}

    LWORK = (int) OPT_LWORK[0];
    printf("sgehrd,length of optimal work : %d \n",LWORK);

    float* WORK = (float*) malloc( (int) sizeof(float) * LWORK );

    sgehrd_ (&n,&ilo, &ihi, a,&lda,tau,WORK ,&LWORK,&info);// calculate Hessenberg
    if(info!=0){printf("sgehrd execution : failed\n");}

    LWORK = -1;
    shseqr_ (&job,&compz,&n,&ilo, &ihi, h,&ldh, wr, wi,z,&ldz, OPT_LWORK ,&LWORK, &info); // query optimal value for shseqr
    if(info!=0){printf("shgeqr lwork=-1 : failed\n");}

    LWORK = (int) OPT_LWORK[0];
    printf("shseqr,length of optimal work : %d \n",LWORK);

    WORK = realloc(WORK,(int) sizeof(float) * LWORK );// possibly realloc with the new LWORK value

    shseqr_ (&job,&compz,&n,&ilo, &ihi, h,&ldh, wr, wi,z,&ldz, WORK ,&LWORK, &info);  // finally obtain eigenvalues
    if(info!=0){printf("shgeqr execution : failed\n");}

    free(OPT_LWORK);
    free(WORK);

    free(a);
    free(tau);
    free(h);
    free(wr);
    free(wi);
    free(z);

}

编译gcc main.c -o main -llapack -lblas

我的输出是：

sgehrd,length of optimal work : 320 
shseqr,length of optimal work : 10
于 2014-11-16T17:02:50.590 回答