Are LAPACKE_cgesdd and LAPACKE_cgesvd SVD calculations reliable?

Question

I'm using LAPACKE_cgesdd and LAPACKE_cgesvd to compute the singular values of a matrix. Both the routines have the option to compute the singular values only. The problem I have is that, in the following four test cases:

1. Full SVD with LAPACKE_cgesdd;
2. Full SVD with LAPACKE_cgesvd;
3. Singular values only with LAPACKE_cgesdd;
4. Singular values only with LAPACKE_cgesvd

I receive different singular values. In particular:

Test, 3 x 4 matrix

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);

Full SVD with LAPACKE_cgesdd

17.8592720031738 11.4463796615601 6.74482488632202

Full SVD with LAPACKE_cgesvd

17.8651084899902 11.3695945739746 6.83876800537109

Singular values only with LAPACKE_cgesdd

17.8592758178711 11.4463806152344 6.74482440948486

Singular values only with LAPACKE_cgesvd

17.8705902099609 11.5145053863525 6.82878828048706

As it can be seen, even for the same routine, the results can change since the third significant digit when switching from full SVD to singular values only.

My questions are:

Is this reasonable?

Am I doing something wrong?

Thank you in advance for any help.

This is the code I'm using:

#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);
}

compiled as

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