您所在位置：主题：调用CLAPACK中函数dsygvx的一个问题
librightly
[专家分：0] 发布于 2010-06-14 20:56:00
HI
我在调用CLapack库的dsygvx计算广义特征值时，遇到如下问题：
矩阵A为
4.00000 0.000000 2.00000 0.000000 1.00000 0.000000 2.00000 0.000000
0.000000 4.00000 0.000000 2.00000 0.000000 1.00000 0.000000 2.00000
2.00000 0.000000 4.00000 0.000000 2.00000 0.000000 1.00000 0.000000
0.000000 2.00000 0.000000 4.00000 0.000000 2.00000 0.000000 1.00000
1.00000 0.000000 2.00000 0.000000 4.00000 0.000000 2.00000 0.000000
0.000000 1.00000 0.000000 2.00000 0.000000 4.00000 0.000000 2.00000
2.00000 0.000000 1.00000 0.000000 2.00000 0.000000 4.00000 0.000000
0.000000 2.00000 0.000000 1.00000 0.000000 2.00000 0.000000 4.00000
矩阵B为8阶的单位阵
代码如下
////////
#include <iostream>
#include <fstream>
#include <iomanip>
#include <time.h>
#include <math.h>
extern "C"
{
#include "f2c.h"
#include "clapack.h"
//#include "mkl_lapack.h"
void* malloc(size_t n);
}
using namespace std;
void print_matrix( char* desc, int m, int n, double* a, int lda );
extern const char* KFilename="f:\\K1.txt";
extern const char* MFilename="f:\\M1.txt";
extern const char* OutFilename="f:\\outf.txt";
ifstream Kinf;
ifstream Minf;
ofstream outf;
const int N=8;
int main()
{
Kinf.open(KFilename,ios_base::in);
if ( !Kinf )
{
cerr<<"File could not be opened"<<endl;
exit( 1 );
}
double K[N*N];
int i=0;
while ( Kinf >> K[i] )
i++;
Kinf.close();
Minf.open(MFilename,ios_base::in);
if ( !Minf )
{
cerr<<"File could not be opened"<<endl;
exit( 1 );
}
double M[N*N];
i=0;
while ( Minf >> M[i] )
i++;
Minf.close();
integer n=N;
outf.open(OutFilename,ios_base::out);
outf.close();
print_matrix("K",n,n,K,n);
print_matrix("M",n,n,M,n);
integer lda=max(1,n);
doublereal* a=(doublereal*)malloc( lda*n*sizeof(doublereal) );
a=K;
integer ldb=max(1,n);
doublereal* b=(doublereal*)malloc( ldb*n*sizeof(doublereal) );
b=M;
integer itype=1;
char jobz='V';
char range= 'I';
char uplo='U';
integer lwork=max(1,8*n);
doublereal* work= (doublereal*)malloc( sizeof(doublereal) * max(1,lwork));
integer* iwork= (integer*)malloc( sizeof(integer) * 5*n);
integer info;
integer il=1;
integer iu=2;
if( range=='I' )
std::cin>>il>>iu;
doublereal* abstol=(doublereal*)malloc( sizeof(doublereal) );
integer m=iu-il+1;
doublereal* w = (doublereal*)malloc( sizeof(doublereal) * n) ;
integer ldz=max(1,n);
doublereal* z = (doublereal*)malloc( sizeof(doublereal) * ldz*m) ;
integer* ifail=(integer*)malloc( sizeof(integer) * n) ;
doublereal vl=0.0 ;
doublereal vu=1.0 ;
if( range=='V' )
std::cin>>vl>>vu;
clock_t start,finish;
start=clock();
dsygvx_(&itype,&jobz,&range,&uplo,&n,a,&lda,b,&ldb,&vl,&vu,&il,&iu,abstol,&m,w,z,&ldz,work,&lwork,iwork,ifail,&info);
finish=clock();
double duration=double(finish-start)/CLOCKS_PER_SEC;
outf.open(OutFilename,ios_base::out|ios_base::app);
outf<<duration<<" seconds has spent"<<endl;
outf.close();
cout<<info<<' '<<*ifail<<endl;
//cout<<dec<<work[1];
print_matrix( "Eigenvalues", 1, m, w, 1 );
print_matrix( "Eigenvector", n, m, z,ldz);
print_matrix( "aa", n, n, a,n);
print_matrix( "bb", n, n, b,n);
free(work);
free(w);
free(z);
int aaaa;
cin>>aaaa;
return 0;
}
/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
outf.open(OutFilename,ios_base::out|ios_base::app);
outf<<desc<<endl;
//printf( "\n %s\n", desc );
for( int i = 0; i <m; i++ ) {
for( int j = 0; j < n; j++ )
{
outf.precision(6);
outf<<setiosflags(ios::right|ios::showpoint)<<setw(15)<< a[i+j*lda];//printf( " %6.4f", a[i+j*lda] );
}
outf<<endl;//printf( "\n" );
}
outf.close();
}
//////////////
结果如下
特征值：
1.00000 1.00000 3.00000 3.00000 3.00000 3.00000 9.00000 9.00000
特征向量：
-0.0286397 0.499179 0.0888613 0.681537 0.142184 -0.0859920 -2.85212e-017 -0.500000
0.499179 0.0286397 0.139710 -0.0568828 -0.212256 -0.657414 0.500000 -5.55112e-017
0.0286397 -0.499179 -0.118747 0.173301 -0.653249 0.170681 -2.29833e-017 -0.500000
-0.499179 -0.0286397 -0.677115 0.0473130 0.0894260 -0.176863 0.500000 -2.77556e-016
-0.0286397 0.499179 -0.0888613 -0.681537 -0.142184 0.0859920 -2.85212e-017 -0.500000
0.499179 0.0286397 -0.139710 0.0568828 0.212256 0.657414 0.500000 2.77556e-016
0.0286397 -0.499179 0.118747 -0.173301 0.653249 -0.170681 -6.05718e-017 -0.500000
-0.499179 -0.0286397 0.677115 -0.0473130 -0.0894260 0.176863 0.500000 0.000000
变量a的输出值为：
3.00000 -4.30211e-016 -0.302776 -0.160828 0.0325749 0.373545 0.000000 0.000000
0.000000 3.00000 4.00297e-016 -0.573368 -0.0474143 -0.280159 0.190744 0.666667
2.00000 0.000000 2.18511 1.13585 -0.310681 0.186773 0.000000 0.000000
0.000000 2.00000 0.000000 3.52853 -3.00861 -0.000000 -0.762974 0.333333
1.00000 0.000000 2.00000 0.000000 7.28636 7.44760e-016 0.000000 0.000000
0.000000 1.00000 0.000000 2.00000 0.000000 2.33333 -1.88562 0.666667
2.00000 0.000000 1.00000 0.000000 2.00000 0.000000 6.66667 -3.00000
0.000000 2.00000 0.000000 1.00000 0.000000 2.00000 0.000000 4.00000
/////////////////
哪位达人可以帮忙解释下为什么结果有这么多的极小数，变量a貌似也不对
是什么原因造成的啊？
非常感谢！
我在调用CLapack库的dsygvx计算广义特征值时，遇到如下问题：
矩阵A为
4.00000 0.000000 2.00000 0.000000 1.00000 0.000000 2.00000 0.000000
0.000000 4.00000 0.000000 2.00000 0.000000 1.00000 0.000000 2.00000
2.00000 0.000000 4.00000 0.000000 2.00000 0.000000 1.00000 0.000000
0.000000 2.00000 0.000000 4.00000 0.000000 2.00000 0.000000 1.00000
1.00000 0.000000 2.00000 0.000000 4.00000 0.000000 2.00000 0.000000
0.000000 1.00000 0.000000 2.00000 0.000000 4.00000 0.000000 2.00000
2.00000 0.000000 1.00000 0.000000 2.00000 0.000000 4.00000 0.000000
0.000000 2.00000 0.000000 1.00000 0.000000 2.00000 0.000000 4.00000
矩阵B为8阶的单位阵
代码如下
////////
#include <iostream>
#include <fstream>
#include <iomanip>
#include <time.h>
#include <math.h>
extern "C"
{
#include "f2c.h"
#include "clapack.h"
//#include "mkl_lapack.h"
void* malloc(size_t n);
}
using namespace std;
void print_matrix( char* desc, int m, int n, double* a, int lda );
extern const char* KFilename="f:\\K1.txt";
extern const char* MFilename="f:\\M1.txt";
extern const char* OutFilename="f:\\outf.txt";
ifstream Kinf;
ifstream Minf;
ofstream outf;
const int N=8;
int main()
{
Kinf.open(KFilename,ios_base::in);
if ( !Kinf )
{
cerr<<"File could not be opened"<<endl;
exit( 1 );
}
double K[N*N];
int i=0;
while ( Kinf >> K[i] )
i++;
Kinf.close();
Minf.open(MFilename,ios_base::in);
if ( !Minf )
{
cerr<<"File could not be opened"<<endl;
exit( 1 );
}
double M[N*N];
i=0;
while ( Minf >> M[i] )
i++;
Minf.close();
integer n=N;
outf.open(OutFilename,ios_base::out);
outf.close();
print_matrix("K",n,n,K,n);
print_matrix("M",n,n,M,n);
integer lda=max(1,n);
doublereal* a=(doublereal*)malloc( lda*n*sizeof(doublereal) );
a=K;
integer ldb=max(1,n);
doublereal* b=(doublereal*)malloc( ldb*n*sizeof(doublereal) );
b=M;
integer itype=1;
char jobz='V';
char range= 'I';
char uplo='U';
integer lwork=max(1,8*n);
doublereal* work= (doublereal*)malloc( sizeof(doublereal) * max(1,lwork));
integer* iwork= (integer*)malloc( sizeof(integer) * 5*n);
integer info;
integer il=1;
integer iu=2;
if( range=='I' )
std::cin>>il>>iu;
doublereal* abstol=(doublereal*)malloc( sizeof(doublereal) );
integer m=iu-il+1;
doublereal* w = (doublereal*)malloc( sizeof(doublereal) * n) ;
integer ldz=max(1,n);
doublereal* z = (doublereal*)malloc( sizeof(doublereal) * ldz*m) ;
integer* ifail=(integer*)malloc( sizeof(integer) * n) ;
doublereal vl=0.0 ;
doublereal vu=1.0 ;
if( range=='V' )
std::cin>>vl>>vu;
clock_t start,finish;
start=clock();
dsygvx_(&itype,&jobz,&range,&uplo,&n,a,&lda,b,&ldb,&vl,&vu,&il,&iu,abstol,&m,w,z,&ldz,work,&lwork,iwork,ifail,&info);
finish=clock();
double duration=double(finish-start)/CLOCKS_PER_SEC;
outf.open(OutFilename,ios_base::out|ios_base::app);
outf<<duration<<" seconds has spent"<<endl;
outf.close();
cout<<info<<' '<<*ifail<<endl;
//cout<<dec<<work[1];
print_matrix( "Eigenvalues", 1, m, w, 1 );
print_matrix( "Eigenvector", n, m, z,ldz);
print_matrix( "aa", n, n, a,n);
print_matrix( "bb", n, n, b,n);
free(work);
free(w);
free(z);
int aaaa;
cin>>aaaa;
return 0;
}
/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
outf.open(OutFilename,ios_base::out|ios_base::app);
outf<<desc<<endl;
//printf( "\n %s\n", desc );
for( int i = 0; i <m; i++ ) {
for( int j = 0; j < n; j++ )
{
outf.precision(6);
outf<<setiosflags(ios::right|ios::showpoint)<<setw(15)<< a[i+j*lda];//printf( " %6.4f", a[i+j*lda] );
}
outf<<endl;//printf( "\n" );
}
outf.close();
}
//////////////
结果如下
特征值：
1.00000 1.00000 3.00000 3.00000 3.00000 3.00000 9.00000 9.00000
特征向量：
-0.0286397 0.499179 0.0888613 0.681537 0.142184 -0.0859920 -2.85212e-017 -0.500000
0.499179 0.0286397 0.139710 -0.0568828 -0.212256 -0.657414 0.500000 -5.55112e-017
0.0286397 -0.499179 -0.118747 0.173301 -0.653249 0.170681 -2.29833e-017 -0.500000
-0.499179 -0.0286397 -0.677115 0.0473130 0.0894260 -0.176863 0.500000 -2.77556e-016
-0.0286397 0.499179 -0.0888613 -0.681537 -0.142184 0.0859920 -2.85212e-017 -0.500000
0.499179 0.0286397 -0.139710 0.0568828 0.212256 0.657414 0.500000 2.77556e-016
0.0286397 -0.499179 0.118747 -0.173301 0.653249 -0.170681 -6.05718e-017 -0.500000
-0.499179 -0.0286397 0.677115 -0.0473130 -0.0894260 0.176863 0.500000 0.000000
变量a的输出值为：
3.00000 -4.30211e-016 -0.302776 -0.160828 0.0325749 0.373545 0.000000 0.000000
0.000000 3.00000 4.00297e-016 -0.573368 -0.0474143 -0.280159 0.190744 0.666667
2.00000 0.000000 2.18511 1.13585 -0.310681 0.186773 0.000000 0.000000
0.000000 2.00000 0.000000 3.52853 -3.00861 -0.000000 -0.762974 0.333333
1.00000 0.000000 2.00000 0.000000 7.28636 7.44760e-016 0.000000 0.000000
0.000000 1.00000 0.000000 2.00000 0.000000 2.33333 -1.88562 0.666667
2.00000 0.000000 1.00000 0.000000 2.00000 0.000000 6.66667 -3.00000
0.000000 2.00000 0.000000 1.00000 0.000000 2.00000 0.000000 4.00000
/////////////////
哪位达人可以帮忙解释下为什么结果有这么多的极小数，变量a貌似也不对
是什么原因造成的啊？
非常感谢！