我正在实现一个迭代算法,该算法使用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/

10-11 19:02