我在计算稠密非对称矩阵a的特征值。为此,我使用xGEHRD和xHSEQRLapack例程,首先计算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
我做了一些测试,得到了工作数组维数的最优值。如果值相同,我可以简化代码很多(不需要realloc,只需要一个调用来确定LWORK的值,更少的错误检查…)。
我的问题是,对于ILO和IHI的相同矩阵和相同值,我能假设两个例程的值将相等吗?
最佳答案
从sgehrd.f来看,常规WORK
的最佳sgehrd()
大小是N*NB
,其中NB
是
NB = MIN( NBMAX, ILAENV( 1, 'SGEHRD', ' ', N, ILO, IHI, -1 ) )
其中
NBMAX=64
。因此,最佳LWORK
取决于N
、ILO
和IHI
。从shseqr.f来看,计算最佳长度
LWORK
更为复杂:例程slaqr0()
被称为。。。但文件中的文件说:如果LWORK=-1,则SHSEQR执行工作区查询。
在这种情况下,SHSEQR检查输入参数和
估计给定的最佳工作空间大小
N值、ILO值和IHI值。估计数已退回
工作中(1)。没有与LWORK相关的错误消息是
由XERBLA发行。H和Z都不可访问。
WORK
和sgehrd()
的最佳长度可能不同。下面是一个例子:#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);
}
编译人
shseqr()
我的输出是:
sgehrd,length of optimal work : 320
shseqr,length of optimal work : 10