From fde49f246b2dbed746bddfad3212797d1697d511 Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Thu, 30 Jan 2020 21:32:35 -0600 Subject: [PATCH] integrated Rupali's changes - replaced OpenCV matrix math with LAPACK --- src/main.cu | 92 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------- 1 file changed, 67 insertions(+), 25 deletions(-) diff --git a/src/main.cu b/src/main.cu index b6415cb..ffc0896 100644 --- a/src/main.cu +++ b/src/main.cu @@ -50,6 +50,43 @@ typedef struct { }gnome; gnome gnom; +//computing matrix determinant using LU decomposition +template +T mtxdeterminant(T* M, size_t r, size_t c) { + //----DETERMINANT using LU decomposition using LAPACK + /* DGETRF computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges. + The factorization has the form + A = P * L * U + where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n).*/ + int* ipiv = (int*)malloc(sizeof(int) * r); + memset(ipiv, 0, r * sizeof(int)); + LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)r, (int)c, M, (int)r, ipiv); + + //determinant of U + T product = 1; + for (size_t i = 0; i < r; i++) { + for (size_t j = 0; j < r; j++) { + if (i == j) { + product = product * M[i * r + j]; + } + } + } + + //determinant of P + int j; + double detp = 1.; + for (j = 0; j < r; j++) { + if (j + 1 != ipiv[j]) { + // j+1 : following feedback of ead : ipiv is from Fortran, hence starts at 1. + // hey ! This is a transpose ! + detp = -detp; + } + } + T detSw = product * detp * 1; //det(U)*det(P)*det(L) + if (ipiv != NULL) std::free(ipiv); + return detSw; + +} void gpuComputeEignS( size_t g, size_t fea){ //eigen value computation will return r = (nC-1) eigen vectors so new projected data will have dimension of r rather than f @@ -87,7 +124,7 @@ void gpuComputeEignS( size_t g, size_t fea){ int *IPIV = (int*) malloc(sizeof(int) * f); //computing inverse of matrix Sw memset(IPIV, 0, f * sizeof(int)); - LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)f, (int)f, gSw_a, (int)f, IPIV); + LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)f, (int)f, gSw_a, (int)f, IPIV); // DGETRI computes the inverse of a matrix using the LU factorization computed by DGETRF. LAPACKE_sgetri(LAPACK_COL_MAJOR, (int)f, gSw_a, (int)f, IPIV); @@ -146,23 +183,22 @@ void gpuComputeEignS( size_t g, size_t fea){ mtxMul(tempgSw, &gnom.lda[g * r * f], &gnom.Sw[g * f * f], r, f, f,f); float* nSw = (float*)calloc(r * r, sizeof(float)); mtxMultranspose(nSw, tempgSw, &gnom.lda[g * r * f], r, f, r, f); - if(debug){ - std::cout<<"From Eigen function: projected Sb sn Sw"<