Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Angelo_L_1
Beginner
115 Views

Intel MKL dgesvd does not seem to return correct singular vectors

I'm using Intel MKL dgesvd in the code below. When comparing the results with Matlab, the singular values are ok, while the singular vectors are not.

The C++ code I'm using

#include <stdio.h>
#include <stdlib.h>
#include <sstream>
#include <string.h>
#include <iostream>
#include <fstream>
#include <time.h>
#include <mkl.h>
#include "mkl_lapacke.h"

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

using namespace std;

/**********************************/
/* RANDOM MATRIX ENTRY GENERATION */
/**********************************/
double generateRandomMatrixEntry(int N) { return 2000. * ((double)rand() / (double)(RAND_MAX - 0.2) * (1. / (double)N)) + 100.*((double)rand() / (double)(RAND_MAX - 0.2) * (1. / (double)N)) + 24.; }

/********/
/* MAIN */
/********/
int main(int argc, char **argv) {

        int Nrows, Ncols, K, numExecutions;

    char *jobu = "S", *jobvt = "S";     // --- "N" means not to calculate the corresponding singular vectors

    TimingCPU timer;

    if (argc != 5) { printf("usage: svdMKLDouble Nrows Ncols K numExecutions\n"); exit(0); }
        { std::stringstream ss(argv[1]); ss >> Nrows; }
        { std::stringstream ss(argv[2]); ss >> Ncols; }
        { std::stringstream ss(argv[3]); ss >> K; }    
        { std::stringstream ss(argv[4]); ss >> numExecutions; }

        int lda = Nrows, ldu = Nrows, ldvt = Ncols;

        double TotalTime = 0., AverageTime = 0.;

    double *A       = (double*)mkl_malloc(Nrows * Ncols * K * sizeof(double), 64);  // --- Input matrices
    double *A_QUERY = (double*)mkl_malloc(Nrows * Ncols *     sizeof(double), 64);  // --- Input matrix for MKL query
    double *U       = (double*)mkl_malloc(Nrows * Nrows * K * sizeof(double), 64);  // --- Matrices of LEFT singular vectors
    double *S       = (double*)mkl_malloc(        Ncols * K * sizeof(double), 64);  // --- Singular values
    double *V       = (double*)mkl_malloc(Ncols * Ncols * K * sizeof(double), 64);  // --- Matrices of RIGHT singular vectors

    srand(time(NULL));
    //srand(0);
    int N = 5;

    /*************************************************/
    /* QUERY AND ALLOCATION OF THE OPTIMAL WORKSPACE */
    /*************************************************/
        for (int h = 0; h < Nrows * Ncols; h++) A_QUERY = generateRandomMatrixEntry(N);

    int info, lwork = -1;          
    double wkopt;
        dgesvd(jobu, jobvt, &Nrows, &Ncols, A_QUERY, &lda, S, U, &ldu, V, &ldvt, &wkopt, &lwork, &info);
        lwork = (int)wkopt;
        double *work = (double *)malloc(lwork * K * sizeof(double));                

    /**************/
    /* DUMMY CALL */
    /**************/
    for (int i = 0; i < K; i++)
            dgesvd(jobu, jobvt, &Nrows, &Ncols, A + (i * Nrows * Ncols), &lda, S + (i * Ncols), U + (i * Nrows * Nrows), &ldu, V + (i * Ncols * Ncols), &ldvt, work + i * lwork, &lwork, &info);

    /*********************/
    /* TIME MEASUREMENTS */
    /*********************/
    for (int k = 0; k < numExecutions; k++) {

        //srand(k);

            for (int h = 0; h < Nrows * Ncols * K; h++) { A = generateRandomMatrixEntry(N); }

            timer.StartCounter();
        #pragma omp parallel for
        for (int i = 0; i < K; i++) {
            dgesvd(jobu, jobvt, &Nrows, &Ncols, A + (i * Nrows * Ncols), &lda, S + (i * Ncols), U + (i * Nrows * Nrows), &ldu, V + (i * Ncols * Ncols), &ldvt, work + i * lwork, &lwork, &info);
            // --- Convergence check            
            if (info > 0) { printf("Convergence failure.\n"); exit(1); }
        }
        #pragma omp barrier   
            TotalTime += timer.GetCounter(); // --- TotalTime in ms

    }

    AverageTime = TotalTime / numExecutions;
    std::cout << std::scientific << Nrows << " " << Ncols << " " << K << " " << numExecutions << " " << AverageTime << " " << std::endl;

    /*****************************/
    /* SAVINGS FOR RESULT CHECKS */
    /*****************************/
    // --- Saving the last matrix
    saveCPUrealtxt<double>(A + ((K - 1) * Nrows * Ncols), "Source_Matrix.txt", Nrows * Ncols);  
    // --- Saving the SVD
    saveCPUrealtxt<double>(U + ((K - 1) * Nrows * Nrows), "U.txt", Nrows * Nrows);  
    saveCPUrealtxt<double>(V + ((K - 1) * Ncols * Ncols), "V.txt", Ncols * Ncols);  
    saveCPUrealtxt<double>(S + ((K - 1) * Ncols),         "S.txt", Ncols        );  

    return 0;
}

The compilation command

icpc -lpthread -lmkl_sequential -lmkl_core -lmkl_intel_lp64 -lrt -openmp svdMKLDouble.cpp TimingCPU.cpp InputOutput.cpp -O3 -o svdMKLDouble

Code usage as...

./svdMKLDouble 8 8 1 1

The Matlab code I'm using as reference

clear all
close all
clc

load Source_Matrix.txt
load U.txt
load V.txt
load S.txt

Source_Matrix = reshape(Source_Matrix, 8, 8);
U             = reshape(U, 8, 8);
V             = reshape(V, 8, 8);

[UU, SS, VV]  = svd(Source_Matrix);

100 * sqrt(sum(abs(diag(SS) - S).^2) / sum(abs(diag(SS)).^2))

100 * sqrt(sum(sum((abs(Source_Matrix - UU * SS * VV.').^2))) / sum(sum((abs(Source_Matrix).^2))))

100 * sqrt(sum(sum((abs(Source_Matrix - U  * diag(S)  * V).^2))) / sum(sum((abs(Source_Matrix).^2)))) 

The comparison results

ans =

    0.0512        ---> Singular values correctly computed


ans =

   5.1387e-13     ---> Matlab SVD works correctly


ans =

  132.3127        ---> Intel MKL singular values do not seem to be correct

The question: What I'm doing wrong?

0 Kudos
0 Replies
Reply