我正在实现一个迭代算法,该算法使用LAPACK进行PSD投影(这并不重要,关键是我一遍又一遍地调用此函数):
void useLAPACK(vector<double>& x, int N){
/* Locals */
int n = N, il, iu, m, lda = N, ldz = N, info, lwork, liwork;
double abstol;
double vl,vu;
int iwkopt;
int* iwork;
double wkopt;
double* work;
/* Local arrays */
int isuppz[N];
double w[N], z[N*N];
/* Negative abstol means using the default value */
abstol = -1.0;
/* Set il, iu to compute NSELECT smallest eigenvalues */
vl = 0;
vu = 1.79769e+308;
/* Query and allocate the optimal workspace */
lwork = -1;
liwork = -1;
dsyevr_( (char*)"Vectors", (char*)"V", (char*)"Upper", &n, &x[0], &lda, &vl, &vu, &il, &iu,
&abstol, &m, w, z, &ldz, isuppz, &wkopt, &lwork, &iwkopt, &liwork,
&info );
lwork = (int)wkopt;
work = (double*)malloc( lwork*sizeof(double) );
liwork = iwkopt;
iwork = (int*)malloc( liwork*sizeof(int) );
/* Solve eigenproblem */
dsyevr_( (char*)"Vectors", (char*)"V", (char*)"Upper", &n, &x[0], &lda, &vl, &vu, &il, &iu,
&abstol, &m, w, z, &ldz, isuppz, work, &lwork, iwork, &liwork,
&info );
/* Check for convergence */
if( info > 0 ) {
printf( "The dsyevr (useLAPACK) failed to compute eigenvalues.\n" );
exit( 1 );
}
/* Print the number of eigenvalues found */
//printf( "\n The total number of eigenvalues found:%2i\n", m );
//print_matrix( "Selected eigenvalues", 1, m, w, 1 );
//print_matrix( "Selected eigenvectors (stored columnwise)", n, m, z, ldz );
//Eigenvectors are returned as stacked columns (in total m)
//Outer sum calculation is fastest.
for(int i = 0; i < N*N; ++i) x[i] = 0;
double lambda;
double vrow1,vrow2;
for(int col = 0; col < m; ++col) {
lambda = w[col];
for (int row1 = 0; row1 < N; ++row1) {
vrow1 = z[N*col+row1];
for(int row2 = 0; row2 < N; ++row2){
vrow2 = z[N*col+row2];
x[row1*N+row2] += lambda*vrow1*vrow2;
}
}
}
free( (void*)iwork );
free( (void*)work );
}
现在我的时间测量表明,第一次通话大约需要4毫秒,但随后增加到100毫秒。这段代码对此有很好的解释吗? x每次都是相同的 vector 。
最佳答案
我想我已经解决了问题。我的算法从零矩阵开始,然后正特征值的数量或多或少为正一半为负。 dsyevr仅使用这些参数计算正特征值和相应的特征 vector 。我想如果全部为零,则不必真正计算任何特征 vector ,从而使算法更快。感谢您提供所有答案,并对缺少的信息感到抱歉。
关于c++ - LAPACK函数在第一次迭代后变慢,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/58647174/