MAGMA 1.4.1 QR factorization incorrect results

Open discussion for MAGMA library (Matrix Algebra on GPU and Multicore Architectures)
Post Reply
pioneerty
Posts: 2
Joined: Sat Apr 05, 2014 3:27 pm

MAGMA 1.4.1 QR factorization incorrect results

Post by pioneerty » Sat Apr 05, 2014 4:48 pm

Hello everyone,

I tried to do QR factorization with MAGMA but several MAGMA functions all give incorrect results. I wonder if it's because I set certain parameters wrong or it's a bug in MAGMA.

My environment:
GPU: NVIDIA GeForce 8600M GS (compute capability 1.1) (it's ancient :-()
OS: Windows 8.1 Pro x64
MAGMA 1.4.1 compiled with MKL 11.1 in Visual Studio 2012

MAGMA's functions I tried:

magma_sgeqrf
magma_sgeqrf_gpu
magma_sgeqrf2_gpu
magma_sgeqrf3_gpu
magma_sgeqp3
magma_sgeqp3_gpu

Below is the code I used for testing (also attached). I started with a random 66-by-66 matrix generated by lapackf77_slarnv and intentionally made it rank deficient by duplicating the last column. I obtain reference results with lapackf77_sgeqrf and lapackf77_sgeqp3.

Code: Select all

#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <string.h>

#include "cuda_runtime.h"
#include "cublas_v2.h"

#include "magma.h"
#include "magma_lapack.h"

#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))

FILE *fdbg;

void sprint(FILE *fo, int m, int n, float *A, int lda)
{
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) 
			fprintf(fo, "%8.4f", A[i + j * lda]);
		fprintf(fo, "\n");
	}
}

int main()
{
	magma_init();

	fopen_s(&fdbg, "dbg_sgeqrf.txt", "w");

	int m = 66, n = 66;
	int lda = m;
	int mn = m * n;

	int idist = 2;
	int iseed[4] = {0, 0, 0, 1};
	float *hA;
	magma_smalloc_cpu(&hA, m * n);
	lapackf77_slarnv(&idist, iseed, &mn, hA);
	
	for (int i = 0; i < m; i++)
		for (int j = 0; j < n; j += 3)
			hA[i + j * lda] = hA[i + (n - 1) * lda];

	fprintf(fdbg, "A [%d, %d] = \n\n", m, n);
	sprint(fdbg, m, n, hA, lda);
	fprintf(fdbg, "\n");

	int ldda = lda;
	float *dA;
	magma_smalloc(&dA, m * n);
	magma_ssetmatrix(m, n, hA, lda, dA, ldda);

	// LAPACK SGEQRF
	{
		float *hA1;
		magma_smalloc_cpu(&hA1, m * n);
		lapackf77_slacpy(MagmaUpperLowerStr, &m, &n, hA, &lda, hA1, &lda);
		float *tau1;
		magma_smalloc_cpu(&tau1, MIN(m,n));
		float *work;
		int lwork = -1, info;
		magma_smalloc_cpu(&work, 1);
		lapackf77_sgeqrf(&m, &n, NULL, &lda, NULL, work, &lwork, &info);
		lwork = int(work[0]);
		magma_free_cpu(work);
		magma_smalloc_cpu(&work, lwork);

		lapackf77_sgeqrf(&m, &n, hA1, &lda, tau1, work, &lwork, &info);

		fprintf(fdbg, ">> LAPACK SGEQRF diag(R) [%d] = \n\n", n);
		for (int i = 0; i < n; i++) fprintf(fdbg, "%8.4f", hA1[i + i * lda]);
		fprintf(fdbg, "\n\n");

		magma_free_cpu(work);
		magma_free_cpu(tau1);
		magma_free_cpu(hA1);
	}

	// MAGMA SGEQRF (CPU)
	{
		float *hA2;
		magma_smalloc_pinned(&hA2, m * n);
		lapackf77_slacpy(MagmaUpperLowerStr, &m, &n, hA, &lda, hA2, &lda);
		float *tau2;
		magma_smalloc_cpu(&tau2, MIN(m,n));
		float *work;
		int lwork = -1, info;
		magma_smalloc_pinned(&work, 1);
		magma_sgeqrf(m, n, NULL, lda, NULL, work, lwork, &info);
		lwork = int(work[0]);
		magma_free_pinned(work);
		magma_smalloc_pinned(&work, lwork);

		magma_sgeqrf(m, n, hA2, lda, tau2, work, lwork, &info);
		cudaDeviceSynchronize();

		fprintf(fdbg, "@@ MAGMA SGEQRF (CPU) diag(R) [%d] = \n\n", n);
		for (int i = 0; i < n; i++) fprintf(fdbg, "%8.4f", hA2[i + i * lda]);
		fprintf(fdbg, "\n\n");

		magma_free_pinned(work);
		magma_free_cpu(tau2);
		magma_free_pinned(hA2);
	}

	// MAGMA SGEQRF_GPU
	{
		int ldda3 = (m + 31) / 32 * 32;
		float *dA3;
		magma_smalloc(&dA3, ldda3 * n);
		magma_scopymatrix(m, n, dA, lda, dA3, ldda3);
		float *tau3;
		magma_smalloc_cpu(&tau3, MIN(m,n));
		int nb = magma_get_sgeqrf_nb(m);
		int tSize = (2 * MIN(m,n) + (n + 31) / 32 * 32 ) * nb;
		float *dT3;
		magma_smalloc(&dT3, tSize);
		int info;

		magma_sgeqrf_gpu(m, n, dA3, ldda3, tau3, dT3, &info);
		cudaDeviceSynchronize();

		float *hRdiag3;
		magma_smalloc_pinned(&hRdiag3, n);
		magma_sgetvector(n, dA3, ldda3 + 1, hRdiag3, 1);
		fprintf(fdbg, "@@ MAGMA SGEQF_GPU diag(R) [%d] = \n\n", n);
		for (int i = 0; i < n; i++) fprintf(fdbg, "%8.4f", hRdiag3[i]);
		fprintf(fdbg, "\n\n");
		magma_free_pinned(hRdiag3);

		magma_free(dT3);
		magma_free_cpu(tau3);
		magma_free(dA3);
	}

	// MAGMA SGEQRF2_GPU
	{
		int ldda4 = lda;
		float *dA4;
		magma_smalloc(&dA4, ldda4 * n);
		magma_scopymatrix(m, n, dA, ldda, dA4, ldda4);
		float *tau4;
		magma_smalloc_cpu(&tau4, MIN(m,n));
		int info;

		magma_sgeqrf2_gpu(m, n, dA4, ldda4, tau4, &info);
		cudaDeviceSynchronize();

		float *hRdiag4;
		magma_smalloc_pinned(&hRdiag4, n);
		magma_sgetvector(n, dA4, ldda4 + 1, hRdiag4, 1);
		fprintf(fdbg, "@@ MAGMA SGEQF2_GPU diag(R) [%d] = \n\n", n);
		for (int i = 0; i < n; i++) fprintf(fdbg, "%8.4f", hRdiag4[i]);
		fprintf(fdbg, "\n\n");
		magma_free_pinned(hRdiag4);

		magma_free_cpu(tau4);
		magma_free(dA4);
	}

	// MAGMA SGEQRF3_GPU
	{
		int ldda5 = (m + 31) / 32 * 32;
		float *dA5;
		magma_smalloc(&dA5, ldda5 * n);
		magma_scopymatrix(m, n, dA, ldda, dA5, ldda5);
		float *tau5;
		magma_smalloc_cpu(&tau5, MIN(m,n));
		int nb = magma_get_sgeqrf_nb(n);
		int tSize = (2 * MIN(m,n) + (n + 31) / 32 * 32 ) * nb;
		float *dT5;
		magma_smalloc(&dT5, tSize);
		int info;

		magma_sgeqrf3_gpu(m, n, dA5, ldda5, tau5, dT5, &info);
		cudaDeviceSynchronize();

		float *hRdiag5;
		magma_smalloc_pinned(&hRdiag5, n);
		magma_sgetvector(n, dA5, ldda5 + 1, hRdiag5, 1);
		fprintf(fdbg, "@@ MAGMA SGEQF3_GPU diag(R) [%d] = \n\n", n);
		for (int i = 0; i < n; i++) fprintf(fdbg, "%8.4f", hRdiag5[i]);
		fprintf(fdbg, "\n\n");
		magma_free_pinned(hRdiag5);

		magma_free(dT5);
		magma_free_cpu(tau5);
		magma_free(dA5);
	}

	// LAPACK SGEQP3
	{
		float *hA6;
		magma_smalloc_cpu(&hA6, m * n);
		lapackf77_slacpy(MagmaUpperLowerStr, &m, &n, hA, &lda, hA6, &lda);
		int *jpvt6;
		magma_malloc_cpu((void **)&jpvt6, sizeof(int) * n);
		memset((void *)jpvt6, 0, sizeof(int) * n);
		float *tau6;
		magma_smalloc_cpu(&tau6, MIN(m,n));
		float *work;
		int lwork = -1, info;
		magma_smalloc_cpu(&work, 1);
		lapackf77_sgeqp3(&m, &n, NULL, &lda, NULL, NULL, work, &lwork, &info);
		lwork = int(work[0]);
		magma_free_cpu(work);
		magma_smalloc_cpu(&work, lwork);

		lapackf77_sgeqp3(&m, &n, hA6, &lda, jpvt6, tau6, work, &lwork, &info);

		fprintf(fdbg, ">> LAPACK SGEQP3 diag(R) [%d] = \n\n", n);
		for (int i = 0; i < n; i++) fprintf(fdbg, "%8.4f", hA6[i + i * lda]);
		fprintf(fdbg, "\n\n");

		magma_free_cpu(work);
		magma_free_cpu(tau6);
		magma_free_cpu(jpvt6);
		magma_free_cpu(hA6);
	}

	// MAGMA SGEQP3 (CPU)
	{
		float *hA7;
		magma_smalloc_pinned(&hA7, m * n);
		lapackf77_slacpy(MagmaUpperLowerStr, &m, &n, hA, &lda, hA7, &lda);
		int *jpvt7;
		magma_malloc_cpu((void **)&jpvt7, sizeof(int) * n);
		memset((void *)jpvt7, 0, sizeof(int) * n);
		float *tau7;
		magma_smalloc_cpu(&tau7, MIN(m,n));
		float *work;
		int lwork = -1, info;
		magma_smalloc_pinned(&work, 1);
		magma_sgeqp3(m, n, NULL, lda, NULL, NULL, work, lwork, &info);
		lwork = int(work[0]);
		magma_free_pinned(work);
		magma_smalloc_pinned(&work, lwork);

		magma_sgeqp3(m, n, hA7, lda, jpvt7, tau7, work, lwork, &info);
		cudaDeviceSynchronize();

		fprintf(fdbg, "@@ MAGMA SGEQP3 (CPU) diag(R) [%d] = \n\n", n);
		for (int i = 0; i < n; i++) fprintf(fdbg, "%8.4f", hA7[i + i * lda]);
		fprintf(fdbg, "\n\n");

		magma_free_pinned(work);
		magma_free_cpu(tau7);
		magma_free_cpu(jpvt7);
		magma_free_pinned(hA7);
	}

	// MAGMA SGEQP3_GPU
	{
		int ldda8 = ldda;
		float *dA8;
		magma_smalloc(&dA8, ldda8 * n);
		magma_scopymatrix(m, n, dA, ldda, dA8, ldda8);
		int *jpvt8;
		magma_malloc_cpu((void **)&jpvt8, sizeof(int) * n);
		memset((void *)jpvt8, 0, sizeof(int) * n);
		float *dtau8;
		magma_smalloc(&dtau8, MIN(m,n));
		float *dwork;
		int lwork = MAX(2 * MIN(m,n) + (n + 1) * magma_get_sgeqp3_nb(m), n * n + n), info;
		magma_smalloc(&dwork, lwork);

		magma_sgeqp3_gpu(m, n, dA8, ldda8, jpvt8, dtau8, dwork, lwork, &info);
		cudaDeviceSynchronize();

		float *hRdiag8;
		magma_smalloc_pinned(&hRdiag8, n);
		magma_sgetvector(n, dA8, ldda8 + 1, hRdiag8, 1);
		fprintf(fdbg, "@@ MAGMA SGEQP3_GPU diag(R) [%d] = \n\n", n);
		for (int i = 0; i < n; i++) fprintf(fdbg, "%8.4f", hRdiag8[i]);
		fprintf(fdbg, "\n\n");
		magma_free_pinned(hRdiag8);

		magma_free(dwork);
		magma_free(dtau8);
		magma_free(dA8);
	}

	magma_free(dA);
	magma_free_cpu(hA);
	fclose(fdbg);

	magma_finalize();

	printf("DONE.");
	getchar();

	return EXIT_SUCCESS;
}
I output the diagonal elements of the resulting R matrix. Below is what I got.

Code: Select all

>> LAPACK SGEQRF diag(R) [66] = 

 -4.6936 -4.5387  4.9423 -0.0000  4.5067  4.7854  0.0000 -4.6959  4.6888 -0.0000  4.2137 -4.7521 -0.0000 -4.3230  4.1481 -0.0000  4.3407 -4.1273  0.0000 -3.8516 -3.6816 -0.0000 -3.9246 -3.4808  0.0000  3.9632 -3.5478  0.0000 -3.2878  3.6462  0.0000  3.7737 -3.0029  0.0000 -2.9937 -3.3786  0.0000  2.8499 -3.5509  0.0000  2.8493  3.0615  0.0000 -1.9662  3.2754  0.0000  2.4939  2.3461  0.0000 -2.4883  2.5237  0.0000  1.7641  2.1774  0.0000 -2.3135 -0.9234  0.0000  0.9351 -1.7379  0.0000 -0.7272 -0.7449 -0.0000  1.2691  0.0000

@@ MAGMA SGEQRF (CPU) diag(R) [66] = 

 -4.6936 -4.5387  4.9423 -0.0000  4.5067  4.7854  0.0000 -4.6959  4.6888 -0.0000  4.2137 -4.7521 -0.0000 -4.3230  4.1481 -0.0000  4.3407 -4.1273  0.0000 -3.8516 -3.6816 -0.0000 -3.9246 -3.4808  0.0000  3.9632 -3.5478  0.0000 -3.2878  3.6462  0.0000 -3.6621 -3.0308 -0.0000  2.9702 -3.3875 -0.0000  2.9556 -3.4710 -0.0000  2.8130  3.2044  0.0000  2.0256  3.1639 -0.0000  2.5169 -2.4163  0.0000 -2.2987  2.6618  0.0000  1.9230 -2.2114  0.0000 -1.6346 -1.2347  0.0000  1.2038  1.4904  0.0000  1.3782 -1.3146 -0.0000  0.5518  0.0000

@@ MAGMA SGEQF_GPU diag(R) [66] = 

  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.5518  0.0000

@@ MAGMA SGEQF2_GPU diag(R) [66] = 

 -4.6936 -4.5387  4.9423 -0.0000  4.5067  4.7854  0.0000 -4.6959  4.6888 -0.0000  4.2137 -4.7521 -0.0000 -4.3230  4.1481 -0.0000  4.3407 -4.1273  0.0000 -3.8516 -3.6816 -0.0000 -3.9246 -3.4808  0.0000  3.9632 -3.5478  0.0000 -3.2878  3.6462  0.0000 -3.6621 -3.0308 -0.0000  2.9702 -3.3875 -0.0000  2.9556 -3.4710 -0.0000  2.8130  3.2044  0.0000  2.0256  3.1639 -0.0000  2.5169 -2.4163  0.0000 -2.2987  2.6618  0.0000  1.9230 -2.2114  0.0000 -1.6346 -1.2347  0.0000  1.2038  1.4904  0.0000  1.3782 -1.3146 -0.0000  0.5518  0.0000

@@ MAGMA SGEQF3_GPU diag(R) [66] = 

  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.5518  0.0000

>> LAPACK SGEQP3 diag(R) [66] = 

 -5.2746 -5.1853  5.1076  5.0561 -5.0209 -5.0157 -4.9575 -4.8133  4.7857 -4.6638 -4.6296 -4.5563 -4.5463 -4.5169  4.4441 -4.4033 -4.3977  4.3234  4.1847 -4.1579 -4.0977 -4.0358 -3.9496  3.9434  3.7999  3.7574  3.6619  3.6303 -3.5383  3.4396 -3.4219  3.3419 -3.3332  3.1291  3.0927  3.0133 -3.0078 -2.8707  2.8486  2.7572  2.5879 -2.4765 -2.2840  2.1537 -0.0000  0.0000 -0.0000  0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

@@ MAGMA SGEQP3 (CPU) diag(R) [66] = 

 -5.2746 -5.1853  5.1076  5.0561 -5.0209 -5.0157 -4.9575 -4.8133  4.7857 -4.6638 -4.6296 -4.5563 -4.5463 -4.5169  4.4441 -4.4033 -4.3977  4.3234  4.1847 -4.1579 -4.0977 -4.0358 -3.9496  3.9579 -3.8270  3.7560  3.6750 -0.0000 -0.0000  0.0000 -3.4460 -3.2853 -3.2604 -3.1821  3.0965  3.0949  3.3624 -2.8824 -2.6928 -2.6222 -3.0478  2.5121  2.4515  2.3511 -3.7577  0.0000  0.0000 -0.0000  0.0000 -0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  1.3152 -1.2409  1.0329  0.8041 -0.7194  0.3682 -0.2283 -0.1823 -0.0486

@@ MAGMA SGEQP3_GPU diag(R) [66] = 

 -5.2746 -5.1853  5.1076  5.0561 -5.0209 -5.0157 -4.9575 -4.8133  4.7857 -4.6638 -4.6296 -4.5563 -4.5463 -4.5169  4.4441 -4.4033 -4.3977  4.3234  4.1847 -4.1579 -4.0977 -4.0358 -3.9496  3.9434  3.7999  3.7574  3.6619  3.6303 -3.5383  3.4396 -3.4219  3.3419 -3.3332  3.1291  3.0927  3.0133 -3.0078 -2.8707  2.8486  2.7572  2.5879 -2.4765 -2.2840  2.1537 -0.0000 -0.0000 -0.0000 -5.1853  5.1076  5.0561 -5.0209 -5.0157 -4.9575 -4.8133  4.7857 -4.6638 -4.6296 -4.5563 -4.5463 -4.5169  4.4441 -4.4033 -4.3977  4.3234  4.1847 -4.1579
The magma_sgeqp3_gpu function wrote a few lines of error message to console:

Code: Select all

 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
 -- recompute dnorms --
CUDA runtime error: invalid argument (11) in magma_slaqps2_gpu at D:/magma-1.4.1/magmablas/slaqps2_gpu.cu:284
I tried some other test cases, and I only checked the diagonal element of the resulting R matrix.

[*] magma_sgeqrf_gpu and magma_sgeqrf3_gpu will write 1's on the diagonal of the R matrix once the input matrix size is bigger than 32-by-32 and non-divisible by 32. Say the input matrix is 66-by-66. Then the first 64 diagonal elements of R will be 1s and the last two diagonal elements of R will be correct numbers. With matrix size smaller than 32-by-32, the two functions seem to do fine.

[*] magma_sgeqrf and magma_sgeqrf2_gpu seem to do fine with full-rank matrices, or rank-deficient matrices smaller than 32-by-32.

[*] With matrices bigger than 63-by-63 magma_sgeqp3 has the first few diagonal elements of R correct, then the numbers go wrong.

[*] magma_sgeqp3_gpu gives incorrect results once the matrix is rank-deficient or nearly-rank-deficient. It fail to give correct result on a Hilbert matrix of order 2, with a "recompute dnorms" message and no CUDA runtime error, and it fails on a Hilbert matrix of order 3, with CUDA runtime error messages.

In addition, I tried to run the MAGMA testing samples of these QR factorization functions. I copied the release-built executables and MKL dlls to a folder and run the exes. I couldn't replicate what I got above with the testing samples, but I got a few other failing cases with non-square matrices.

This is for magma_sgeqrf_gpu

Code: Select all

C:\Users\pioneerty\Desktop\testing>testing_sgeqrf_gpu -N 3,2 -c2 -l --version 1
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqrf_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||Ax-b||_F/(N*||A||_F*||x|
|_F)
================================================================================
====
    3     2      1.#J (   0.00)      0.00 (   0.00)   1.50e-001  failed

C:\Users\pioneerty\Desktop\testing>testing_sgeqrf_gpu -N 3,2 -c2 --version 1
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqrf_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||Ax-b||_F/(N*||A||_F*||x|
|_F)
================================================================================
====
    3     2     ---   (  ---  )      0.00 (   0.00)   5.91e-008
I don't know why the "-l" parameter made such a difference.

This is for magma_sgeqrf3_gpu

Code: Select all

C:\Users\pioneerty\Desktop\testing>testing_sgeqrf_gpu -N 3,2 -c2 -l --version 3
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqrf_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||Ax-b||_F/(N*||A||_F*||x|
|_F)
================================================================================
====
    3     2      1.#J (   0.00)      0.00 (   0.00)   1.50e-001  failed

C:\Users\pioneerty\Desktop\testing>testing_sgeqrf_gpu -N 3,2 -c2 --version 3
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqrf_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||Ax-b||_F/(N*||A||_F*||x|
|_F)
================================================================================
====
    3     2     ---   (  ---  )      0.00 (   0.00)   5.77e-008
This is for magma_sgeqrf2_gpu

Code: Select all

C:\Users\pioneerty\Desktop\testing>testing_sgeqrf_gpu -N 3,2 -c2 -l --version 2
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqrf_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||R||_F / ||A||_F
=======================================================================
    3     2      1.#J (   0.00)      0.00 (   0.00)   0.00e+000

C:\Users\pioneerty\Desktop\testing>testing_sgeqrf_gpu -N 3,2 -c -l --version 2
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqrf_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||R-Q'A||_1 / (M*||A||_1*e
ps) ||I-Q'Q||_1 / (M*eps)
================================================================================
=========================
    3     2      1.#J (   0.00)      0.00 (   0.00)   6.07e-001
     1.08e+000
This is for magma_sgeqp3_gpu

Code: Select all

C:\Users\pioneerty\Desktop\testing>testing_sgeqp3_gpu -N 2,2 -c -l
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqp3_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||A*P - Q*R||_F
=====================================================================
    2     2      0.00 (   0.00)      0.00 (   0.00)   4.43e-002

C:\Users\pioneerty\Desktop\testing>testing_sgeqp3_gpu -N 2,2 -c2 -l
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqp3_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||A*P - Q*R||_F
=====================================================================
    2     2      0.00 (   0.00)      0.00 (   0.00)   4.43e-002

C:\Users\pioneerty\Desktop\testing>testing_sgeqp3_gpu -N 3,2 -c -l
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqp3_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||A*P - Q*R||_F
=====================================================================
    3     2      0.00 (   0.00)      0.00 (   0.01)   8.16e-008

C:\Users\pioneerty\Desktop\testing>testing_sgeqp3_gpu -N 3,2 -c2 -l
MAGMA 1.4.1 , compiled for CUDA capability >= 1.0
device 0: GeForce 8600M GS, 1000.0 MHz clock, 256.0 MB memory, capability 1.1
Usage: testing_sgeqp3_gpu [options] [-h|--help]

  M     N     CPU GFlop/s (sec)   GPU GFlop/s (sec)   ||A*P - Q*R||_F
=====================================================================
    3     2      0.00 (   0.00)      0.00 (   0.00)   8.16e-008
Thanks in advance! Any help is appreciated.
Attachments
test_sgeqrf.cpp
code testing magma's various qr factorization
(7.72 KiB) Downloaded 101 times
Top
mgates3
Posts: 918
Joined: Fri Jan 06, 2012 2:13 pm

Re: MAGMA 1.4.1 QR factorization incorrect results

Post by mgates3 » Mon Apr 07, 2014 3:50 pm

Use dgeqrf2_gpu for LAPACK-complaint output.

dgeqrf_gpu and dgeqrf3_gpu store 1s on the diagonal and 0s above the diagonal of each block of V. Thus, the upper triangle is not the full R matrix as in LAPACK. The diagonal blocks of R or R^{-1} are stored in dT. This makes using QR in zgels faster.

I don't know if this was the cause or not, but it is possible for magma_dgeqp3 to produce a different, but correct, factorization than lapack dgeqp3. If two columns have similar norms during the factorization, slight differences in rounding could cause one or the other column to be chosen as a pivot. See if the JPVT vectors are the same or not. That's why we check A*P - Q*R, instead of directly comparing with lapack's results.

We will look more closely into these results to see if we can replicate them and what bugs they might expose in MAGMA. Thanks.

-mark
Top
pioneerty
Posts: 2
Joined: Sat Apr 05, 2014 3:27 pm

Re: MAGMA 1.4.1 QR factorization incorrect results

Post by pioneerty » Thu Apr 10, 2014 10:51 pm

Thanks a lot for the comments!

I was confused by the documentation. The description in sgeqrf_gpu.cpp does say that "the corresponding parts of the upper triangular R are inverted and stored separately in dT." However, the description of argument dA still says that it stores the R matrix on exit. I took it for granted. I guess it's the same for sgeqrf3_gpu.

Nevertheless, how can I recover the R matrix or at least the diagonal of R after magma_sgeqrf_gpu and magma_sgeqrf3_gpu? These two QR versions seem very attractive as they can make subsequent application of Q much faster. On the other hand, is there or will there be a version of sgeqp3 in MAGMA that makes application of Q faster?

For magma_sgeqrf and magma_sgeqrf2_gpu, I tried to check whether the returned results meet Q * R = A. It turned out that MAGMA produces a different yet correct QR factorization as compared to MKL's. It was my silly mistake to neglect this possibility.

I checked the JPVT arrays returned by magma_sgeqp3 and magma_sgeqp3_gpu, but they differed from the one returned by MKL's sgeqp3. I also checked the resulting factorization of magma_sgeqp3 and magma_sgeqp3_gpu by computing || A * P - Q * R || / (|| A || * M). It seems that A * P = Q * R was not met. Below is my output (updated testing code attached):

Code: Select all

>> LAPACK SGEQP3 diag(R) [66] = 

 -5.2746 -5.1853  5.1076  5.0561 -5.0209 -5.0157 -4.9575 -4.8133  4.7857 -4.6638 -4.6296 -4.5563 -4.5463 -4.5169  4.4441 -4.4033 -4.3977  4.3234  4.1847 -4.1579 -4.0977 -4.0358 -3.9496  3.9434  3.7999  3.7574  3.6619  3.6303 -3.5383  3.4396 -3.4219  3.3419 -3.3332  3.1291  3.0927  3.0133 -3.0078 -2.8707  2.8486  2.7572  2.5879 -2.4765 -2.2840  2.1537 -0.0000  0.0000 -0.0000  0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

   LAPACK SGEQP3 JPVT [66] = 

      12      56      63      51      32      30      42      54       8       3       6       9      65      14      17      39      18      21      36      26      62      22      45      15      33      48      57      23      53      11       5      20      35      24      60      41      50      38       2      59      44      47      29      27      64      66      13      31      49      46       4       7      43       1      55      16      19      58      40      25      61      34      10      28      52      37

  ||A * P - Q * R|| / (||A|| * M) = 4.444211e-009

@@ MAGMA SGEQP3 (CPU) diag(R) [66] = 

 -5.2746 -5.1853  5.1076  5.0561 -5.0209 -5.0157 -4.9575 -4.8133  4.7857 -4.6638 -4.6296 -4.5563 -4.5463 -4.5169  4.4441 -4.4033 -4.3977  4.3234  4.1847 -4.1579 -4.0977 -4.0358 -3.9496  3.9579 -3.8270  3.7560  3.6750 -0.0000 -0.0000  0.0000 -3.4460 -3.2853 -3.2604 -3.1821  3.0965  3.0949  3.3624 -2.8824 -2.6928 -2.6222 -3.0478  2.5121  2.4515  2.3511 -3.7577  0.0000  0.0000 -0.0000  0.0000 -0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  1.3152 -1.2409  1.0329  0.8041 -0.7194  0.3682 -0.2283 -0.1823 -0.0486

   MAGMA SGEQP3 (CPU) JPVT [66] = 

      12      56      63      51      32      30      42      54       8       3       6       9      65      14      17      39      18      21      36      26      62      66      45      15      48      33      57      28      31      25      23      53      35      24      59      29       2      41      60      50      20       7      47      44       5      43      37      34      49      40      52      55      46      58      61      64      22      13      38       4      27      11       1      19      10      16

  ||A * P - Q * R|| / (||A|| * M) = 5.016793e-003

@@ MAGMA SGEQP3_GPU diag(R) [66] = 

 -5.2746 -5.1853  5.1076  5.0561 -5.0209 -5.0157 -4.9575 -4.8133  4.7857 -4.6638 -4.6296 -4.5563 -4.5463 -4.5169  4.4441 -4.4033 -4.3977  4.3234  4.1847 -4.1579 -4.0977 -4.0358 -3.9496  3.9434  3.7999  3.7574  3.6619  3.6303 -3.5383  3.4396 -3.4219  3.3419 -3.3332  3.1291  3.0927  3.0133 -3.0078 -2.8707  2.8486  2.7572  2.5879 -2.4765 -2.2840  2.1537 -0.0000 -0.0000 -0.0000 -5.1853  5.1076  5.0561 -5.0209 -5.0157 -4.9575 -4.8133  4.7857 -4.6638 -4.6296 -4.5563 -4.5463 -4.5169  4.4441 -4.4033 -4.3977  4.3234  4.1847 -4.1579

   MAGMA SGEQP3_GPU JPVT [66] = 

      12      56      63      51      32      30      42      54       8       3       6       9      65      14      17      39      18      21      36      26      62      22      45      15      33      48      57      23      53      11       5      20      35      24      60      41      50      38       2      59      44      47      29      27      28      46       7      31      49      37       4      52      43       1      55      16      19      58      40      25      61      34      10      64      13      66

  ||A * P - Q * R|| / (||A|| * M) = 8.036694e-003
One more question, does MAGMA's sgeqp3 or sgeqp3_gpu guarantee that the diagonal of R decreases in absolute value? It's not documented but I'm afraid I take it for granted again.

Thank you!
Attachments
test_sgeqrf_rev.cpp
qr test w/ check
(14.81 KiB) Downloaded 111 times
Top
mgates3
Posts: 918
Joined: Fri Jan 06, 2012 2:13 pm

Re: MAGMA 1.4.1 QR factorization incorrect results

Post by mgates3 » Tue Apr 15, 2014 12:44 pm

Yes, sorry, the documentation for geqrf_gpu and geqrf3_gpu is not accurate in its description of dA. The upper triangle of the diagonal blocks are set to the identity, to allow multiplying by the V block columns (stored below the diagonal) much easier and faster on the GPU.

For magma_dgeqrf3, the diagonal blocks of R are stored in the dT workspace.
For magma_dgeqrf, the diagonal blocks of R are inverted and stored in the dT workspace. You can't easily exactly reconstruct R. You could re-invert the diagonal blocks, incurring some additional rounding errors.

I ran your code and observed the same problem with magma_sgeqp3. However, our tester doesn't exhibit the problem. I haven't yet been able to figure out the difference, but will keep investigating.

Yes, I believe that after geqp3, the diagonal entries should be decreasing in magnitude. It's a consequence of the current algorithm, which is the same as in LAPACK. I don't think that it is inherent in the problem, though -- other rank-revealing QR factorizations may not have this property.

-mark
Top
Post Reply

Return to “User discussion”