Надежны ли вычисления SVD LAPACKE_cgesdd и LAPACKE_cgesvd?

я использую LAPACKE_cgesdd а также LAPACKE_cgesvd вычислить особые значения матрицы. Обе процедуры имеют возможность вычислять только особые значения. Проблема в том, что в следующих четырех тестах:

1. Полный СВД с LAPACKE_cgesdd;
2. Полный СВД с LAPACKE_cgesvd;
3. Единственные значения только с LAPACKE_cgesdd;
4. Единственные значения только с 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);

Полный СВД с LAPACKE_cgesdd

+17,8592720031738
+11,4463796615601
+6,74482488632202

Полный СВД с LAPACKE_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