c++ - LAPACKE_cgesdd 和 LAPACKE_cgesvd SVD 计算是否可靠？

Question

我正在使用LAPACKE_cgesddandLAPACKE_cgesvd来计算矩阵的奇异值。这两个例程都可以选择仅计算奇异值。我遇到的问题是，在以下四个测试用例中：

完整的 SVD 与LAPACKE_cgesdd;
完整的 SVD 与LAPACKE_cgesvd;
奇异值仅与LAPACKE_cgesdd;
奇异值仅与LAPACKE_cgesvd

我收到不同的奇异值。尤其是：

测试，3 x 4矩阵

a[0].real(5.91); a[0].imag(-5.96);
a[1].real(7.09); a[1].imag(2.72);
a[2].real(7.78); a[2].imag(-4.06);
a[3].real(-0.79); a[3].imag(-7.21);
a[4].real(-3.15); a[4].imag(-4.08);
a[5].real(-1.89); a[5].imag(3.27);
a[6].real(4.57); a[6].imag(-2.07);
a[7].real(-3.88); a[7].imag(-3.30);
a[8].real(-4.89); a[8].imag(4.20);
a[9].real(4.10); a[9].imag(-6.70);
a[10].real(3.28); a[10].imag(-3.84);
a[11].real(3.84); a[11].imag(1.19);

完整的 SVDLAPACKE_cgesdd

17.8592720031738 11.4463796615601 6.74482488632202

完整的 SVDLAPACKE_cgesvd

17.8651084899902 11.3695945739746 6.83876800537109

奇异值仅与LAPACKE_cgesdd

17.8592758178711 11.4463806152344 6.74482440948486

奇异值仅与LAPACKE_cgesvd

17.8705902099609 11.5145053863525 6.82878828048706

可以看出，即使对于相同的例程，当仅从完整 SVD 切换到奇异值时，结果也会从第三位有效数字开始发生变化。

我的问题是：

这合理吗？

难道我做错了什么？

预先感谢您的任何帮助。

这是我正在使用的代码：

#include <stdlib.h>
#include <stdio.h>
#include <algorithm>    // std::min
#include <time.h>
#include <complex>
#include <mkl.h>
#include "mkl_lapacke.h"

#include "TimingCPU.h"
#include "InputOutput.h"

//#define FULL_SVD
#define PRINT_MATRIX
#define PRINT_SINGULAR_VALUES
//#define PRINT_LEFT_SINGULAR_VECTORS
//#define PRINT_RIGHT_SINGULAR_VECTORS
#define SAVE_MATRIX
#define SAVE_SINGULAR_VALUES
//#define SAVE_LEFT_SINGULAR_VECTORS
//#define SAVE_RIGHT_SINGULAR_VECTORS

#define GESDD
//#define GESVD

/*************************************************************/
/* PRINT A SINGLE PRECISION COMPLEX MATRIX STORED COLUMNWISE */
/*************************************************************/
void print_matrix_col(char const *desc, MKL_INT Ncols, MKL_INT Nrows, std::complex<float>* a, MKL_INT LDA) {
    printf( "\n %s\n[", desc);
    for(int i = 0; i < Ncols; i++) {
        for(int j = 0; j < Nrows; j++)
            printf("(%6.2f,%6.2f)", a[i * LDA + j].real(), a[i * LDA + j].imag());
        printf( "\n" );
    }
}

/**********************************************************/
/* PRINT A SINGLE PRECISION COMPLEX MATRIX STORED ROWWISE */
/**********************************************************/
void print_matrix_row(char const *desc, int Nrows, int Ncols, std::complex<float>* a, int LDA) {
    printf( "\n %s\n", desc);
    for (int i = 0; i < Ncols; i++) {
        for (int j = 0; j < Ncols; j++)
            printf("%2.10f + 1i * %2.10f ", a[i * LDA + j].real(), a[i * LDA + j].imag());
        printf( ";\n" );
    }
}

/****************************************/
/* PRINT A SINGLE PRECISION REAL MATRIX */
/****************************************/
void print_rmatrix(char const *desc, MKL_INT m, MKL_INT n, float* a, MKL_INT lda ) {
    MKL_INT i, j;
    printf( "\n %s\n", desc );
    for( i = 0; i < m; i++ ) {
        for( j = 0; j < n; j++ ) printf( " %6.2f", a[i*lda+j] );
        printf( "\n" );
    }
}

/********/
/* MAIN */
/********/
int main() {

    const int Nrows = 3;            // --- Number of rows
    const int Ncols = 4;            // --- Number of columns
    const int LDA   = Ncols;
    const int LDU   = Nrows;
    const int LDVT  = Ncols;    

    const int numRuns   = 20;       // --- Number of runs for timing

    TimingCPU timerCPU; 

    // --- Allocating space and initializing the input matrix
    std::complex<float> *a = (std::complex<float> *)malloc(Nrows * Ncols * sizeof(std::complex<float>));
    srand(time(NULL));
//  for (int k = 0; k < Nrows * Ncols; k++) {
//      a[k].real((float)rand() / (float)(RAND_MAX));   
//      a[k].imag((float)rand() / (float)(RAND_MAX));   
//  }
    a[0].real(5.91); a[0].imag(-5.96);
    a[1].real(7.09); a[1].imag(2.72);
    a[2].real(7.78); a[2].imag(-4.06);
    a[3].real(-0.79); a[3].imag(-7.21);
    a[4].real(-3.15); a[4].imag(-4.08);
    a[5].real(-1.89); a[5].imag(3.27);
    a[6].real(4.57); a[6].imag(-2.07);
    a[7].real(-3.88); a[7].imag(-3.30);
    a[8].real(-4.89); a[8].imag(4.20);
    a[9].real(4.10); a[9].imag(-6.70);
    a[10].real(3.28); a[10].imag(-3.84);
    a[11].real(3.84); a[11].imag(1.19);

    // --- Allocating space for the singular vector matrices
    #ifdef FULL_SVD 
    std::complex<float> *u  = (std::complex<float> *)malloc(Nrows * LDU  * sizeof(std::complex<float>));
    std::complex<float> *vt = (std::complex<float> *)malloc(Ncols * LDVT * sizeof(std::complex<float>));
    #endif

    // --- Allocating space for the singular values 
    float *s = (float *)malloc(std::min(Nrows, Ncols) * sizeof(float));

    #ifdef GESVD
    float *superb = (float *)malloc((std::min(Nrows, Ncols) - 1) * sizeof(float));
    #endif

    // --- Print and/or save input matrix
    #ifdef PRINT_MATRIX 
    print_matrix_row("Matrix (stored rowwise)", Ncols, Nrows, a, LDA);
    #endif
    #ifdef SAVE_MATRIX
    saveCPUcomplextxt(a, "/home/angelo/Project/SVD/MKL/a.txt", Ncols * Nrows);
    #endif

    // --- Compute singular values
    MKL_INT info;
    float   timing = 0.f;   
    for (int k = 0; k < numRuns; k++) {
        timerCPU.StartCounter();
        // --- The content of the input matrix a is destroyed on output
        #if defined(FULL_SVD) && defined(GESDD)
        printf("Running Full SVD - GESDD\n");
        MKL_INT info = LAPACKE_cgesdd(LAPACK_ROW_MAJOR, 'A', Nrows, Ncols, (MKL_Complex8 *)a, LDA, s, (MKL_Complex8 *)u, LDU, (MKL_Complex8 *)vt, LDVT);
        #endif      
        #if !defined(FULL_SVD) && defined(GESDD)        
        printf("Running singular values only - GESDD\n");
        MKL_INT info = LAPACKE_cgesdd(LAPACK_ROW_MAJOR, 'N', Nrows, Ncols, (MKL_Complex8 *)a, LDA, s, NULL, LDU, NULL, LDVT);
        #endif
        #if defined(FULL_SVD) && defined(GESVD)     
        printf("Running Full SVD - GESVD\n");
        MKL_INT info = LAPACKE_cgesvd(LAPACK_ROW_MAJOR, 'A', 'A', Nrows, Ncols, (MKL_Complex8 *)a, LDA, s,
         (MKL_Complex8 *)u, LDU, (MKL_Complex8 *)vt, LDVT, superb);
        #endif
        #if !defined(FULL_SVD) && defined(GESVD)
        printf("Running singular values only - GESVD\n");
        MKL_INT info = LAPACKE_cgesvd(LAPACK_ROW_MAJOR, 'N', 'N', Nrows, Ncols, (MKL_Complex8 *)a, LDA, s,
         NULL, LDU, NULL, LDVT, superb);
        #endif
        if(info > 0) { // --- Check for convergence
            printf( "The algorithm computing SVD failed to converge.\n" );
            exit(1);
        }
        timing = timing + timerCPU.GetCounter();
    }
    printf("Timing = %f\n", timing / numRuns);

    // --- Print and/or save singular values
    #ifdef PRINT_SINGULAR_VALUES
    print_rmatrix("Singular values", 1, Ncols, s, 1);
    #endif
    #ifdef SAVE_SINGULAR_VALUES 
    saveCPUrealtxt(s, "/home/angelo/Project/SVD/MKL/s.txt", std::min(Nrows, Ncols));
    #endif

    // --- Print left singular vectors
    #ifdef PRINT_LEFT_SINGULAR_VECTORS  
    print_matrix_col("Left singular vectors (stored columnwise)", Ncols, Ncols, u, LDU);
    #endif
    #if defined(FULL_SVD) && defined(SAVE_LEFT_SINGULAR_VECTORS) 
    saveCPUcomplextxt(u, "/home/angelo/Project/SVD/MKL/u.txt", Nrows * LDU);
    #endif

    // --- Print right singular vectors
    #ifdef PRINT_RIGHT_SINGULAR_VECTORS 
    print_matrix_col("Right singular vectors (stored rowwise)", Ncols, Nrows, vt, LDVT);
    #endif
    #if defined(FULL_SVD) && defined(SAVE_RIGHT_SINGULAR_VECTORS) 
    saveCPUcomplextxt(vt, "/home/angelo/Project/SVD/MKL/vt.txt", Ncols * LDVT);
    #endif

    exit(0);
}

编译为

g++ -DMKL_ILP64 -fopenmp -m64 -I$MKLROOT/include svdMKLComplexSingle.cpp TimingCPU.cpp InputOutput.cpp -L$MKLROOT/lib/intel64 -lmkl_intel_ilp64 -lmkl_core -lmkl_sequential -lpthread -lm -fno-exceptions

c++ - LAPACKE_cgesdd 和 LAPACKE_cgesvd SVD 计算是否可靠？

0 回答 0

Related

Reference