Discussion :

[alexgille]

Bonjour à tous,

j'ai une résolution de système d'équations non carré du type Ax = b à faire sur un GPU. Pour cela j'utilise cublas sous unix mais je n'arrive pas aux résultats que je devrais avoir.

Quelqu'un s'y connais un peu ? je peux poster mon code ici ?

Merci beaucoup !

[kiluak]

Tu ne perds rien à nous montrer (si ton code ne fait pas 500000 lignes de préférence )

[alexgille]

Le code étant un peu gros, je ne poste que la partie qui me pose problème.
Le but de mon programme est de résoudre un système d'équation du type A.x=b en "out of core". C-à-d que la matrice A est stockée dans un fichier et est lue bloc par bloc (pas chargée entièrement en mémoire).Je fais pareil pour le vecteur b.

Le seul problème est que Cublas préconise l'utilisation un stockage "column-major", or en C le stockage est "row-major". Cependant lorsque je fait mon calcul AiT.Ai (Ai étant un bloc de A), j'obtiens de bons résultats.

Mais l'utilisation de la fonction cublasDtrsm() ne renvois que des valeurs fausses.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
/* Chargement d'un bloc de la matrice en mémoire */
		read_bloc(file_a, file_b, Ai, bi, NC, BLOC, &lines_read);
		printf("Chargement de %d / %d lignes depuis '%s' et '%s' ... OK\n", (lines_read), NL, argv[2], argv[3]);
		for(i=0; i < BLOC*NC; ++i) {
			printf("%lf ", *(Ai+i));
 
			if((i+1)%NC == 0) printf("\n");
		}		
		for(i=0; i < BLOC; ++i) {
			printf("%lf\n", *(bi+i));
 
		/* Mise en mémoire GPU du vecteur b */
		printf("Mise en GPU du vecteur b ... ");
		if(cublasSetVector(NC, sizeof(double), bi, 1, dev_bi, 1) != CUBLAS_STATUS_SUCCESS) {
			printf("fail!\n");
			cublasShutdown();
			return 1;
		}
		printf("OK.\n");
 
		/* Mise en mémoire GPU de la matrice Ai */
		printf("Mise en GPU de la matrice ... ");
		if(cublasSetMatrix(BLOC, NC, sizeof(double), Ai, BLOC, dev_Ai, BLOC) != CUBLAS_STATUS_SUCCESS) {
			printf("fail!\n");
			cublasShutdown();
			return 1;
		}
		printf("OK.\n");
 
		/* Calcul de (AiT . Ai) */
		printf("Calcul de AiT.Ai ... ");
		for(i=0; i < NC; ++i) {
			for(j=0; j < NC; ++j) {
				*(AiTAi+i*NC+j) = cublasDdot(BLOC, dev_Ai+i, NC, dev_Ai+j, NC);
			}
		}
		if(cublasGetError() != CUBLAS_STATUS_SUCCESS) {
			getMessage(cublasGetError());
			cublasShutdown();
			return 1;
		}
		printf("OK.\n");
 
		/* Calcul de (AiT . b) */
		printf("Calcul de AiT.b ... ");
		for(i=0; i < BLOC; ++i) {
			*(AiTbi+i) = cublasDdot(BLOC, dev_Ai+i, NC, dev_bi, 1);
		}
		if(cublasGetError() != CUBLAS_STATUS_SUCCESS) {
			getMessage(cublasGetError());
			cublasShutdown();
			return 1;
		}
		printf("OK.\n");
 
		/* Mise en mémoire GPU de la matrice AiTAi */
		printf("Mise en GPU de la matrice AiTAi ... ");
		if(cublasSetMatrix(NC, NC, sizeof(double), AiTAi, NC, dev_AiTAi, NC) != CUBLAS_STATUS_SUCCESS) {
			printf("fail!\n");
			cublasShutdown();
			return 1;
		}
		printf("OK.\n");
 
		/* Mise en mémoire GPU de la matrice AiTbi */
		printf("Mise en GPU du vecteur AiTbi ... ");
		if(cublasSetVector(NC, sizeof(double), AiTbi, 1, dev_AiTbi, 1) != CUBLAS_STATUS_SUCCESS) {
			printf("fail!\n");
			cublasShutdown();
			return 1;
		}
		printf("OK.\n");
 
		/* Résolution du système */
		cublasDtrsm('L', 'U', 'N', 'N', BLOC, 1, (double)1, dev_Ai, BLOC, dev_bi, BLOC);
 
		/* Récupération de la matrice */
		printf("Récupération du résultat en GPU ... ");
		if(cublasGetMatrix(NC, 1, sizeof(double), dev_AiTbi, NC, bi, NC) != CUBLAS_STATUS_SUCCESS) {
			printf("fail!\n");
			getMessage(cublasGetError());
			cublasShutdown();
			return 1;
		}
		printf("OK.\n");

BLOC: Longueur en ligne d'un bloc.
NC: Nombre de colonnes de A.
NL: Nombre de ligne total de la matrice A.

Merci beaucoup !

[alexgille]

Voici le lien vers le forum de nVidia qui explique toute la résolution du système et comment parvenir à cette résolution.

http://forums.nvidia.com/index.php?s...&#entry1242937