testBlas.cpp 3.97 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
// ===================================================================================
// Ce LOGICIEL "ScalFmm" est couvert par le copyright Inria 20xx-2012.
// Inria détient tous les droits de propriété sur le LOGICIEL, et souhaite que
// la communauté scientifique l'utilise afin de le tester et de l'évaluer.
// Inria donne gracieusement le droit d'utiliser ce LOGICIEL. Toute utilisation
// dans un but lucratif ou à des fins commerciales est interdite sauf autorisation
// expresse et préalable d'Inria.
// Toute utilisation hors des limites précisées ci-dessus et réalisée sans l'accord
// expresse préalable d'Inria constituerait donc le délit de contrefaçon.
// Le LOGICIEL étant un produit en cours de développement, Inria ne saurait assurer
// aucune responsabilité et notamment en aucune manière et en aucun cas, être tenu
// de répondre d'éventuels dommages directs ou indirects subits par l'utilisateur.
// Tout utilisateur du LOGICIEL s'engage à communiquer à Inria ses remarques
// relatives à l'usage du LOGICIEL
// ===================================================================================

// ==== CMAKE =====
// @FUSE_BLAS
// ================

21 22 23
#include <iostream>
#include <stdlib.h>

24
#include "../../Src/Utils/FBlas.hpp"
25

26 27
#include "../../Src/Utils/FParameterNames.hpp"

28 29

template <class FReal>
BRAMAS Berenger's avatar
BRAMAS Berenger committed
30
FReal FRandom() { return (FReal(drand48())); }
31 32 33 34 35

/**
 * Test functionality of C - interfaced BLAS functions
 */

36
int main(int argc, char** argv)
37
{
38 39
    FHelpDescribeAndExit(argc, argv, "Simply ensure that blas are compuling and running.");

40
    typedef double FReal;
41 42 43 44 45
	const unsigned int m = 4, n = 4; // to be able to test both, transpose and not transpose operations
	FReal* A = new FReal [m * n]; // matrix: column major ordering
	FReal* x = new FReal [n];
	FReal* y = new FReal [m];
	
46
	const FReal d = FRandom<FReal>();
47 48
	
	for (unsigned int j=0; j<n; ++j) {
49
		x[j] = FRandom<FReal>();
50
		for (unsigned int i=0; i<m; ++i) {
51
			A[j*m + i] = FRandom<FReal>();
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 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
		}
	}

//	std::cout << "A = " << std::endl;
//	for (unsigned int i=0; i<m; ++i) {
//		for (unsigned int j=0; j<n; ++j) 
//			std::cout << A[j*m + i] << " ";
//		std::cout << std::endl;
//	}
//	std::cout << std::endl;	
//	
//	std::cout << "x = " << std::endl;
//	for (unsigned int i=0; i<m; ++i)
//		std::cout << x[i] << std::endl;
//	std::cout << std::endl;	


	
	FReal* z = new FReal [m];
	

	// y = d Ax ////////////////////////////////////
	// cblas
	FBlas::gemv(m, n, d, A, x, y);
	// direct
	for (unsigned int i=0; i<m; ++i) {
		z[i] = FReal(0.);
		for (unsigned int j=0; j<n; ++j) 
			z[i] += A[j*m + i] * x[j];
		z[i] *= d;
	}
	// compare
	std::cout << "\ny = d Ax (zeros are correct)" << std::endl;
	for (unsigned int i=0; i<m; ++i)
		std::cout << z[i] - y[i] << std::endl;


	// y = d A^Tx ////////////////////////////////////
	// cblas
	FBlas::gemtv(m, n, d, A, x, y);
	// direct
	for (unsigned int i=0; i<m; ++i) {
		z[i] = FReal(0.);
		for (unsigned int j=0; j<n; ++j) 
			z[i] += A[i*m + j] * x[j];
		z[i] *= d;
	}
	// compare
	std::cout << "\ny = d A^Tx (zeros are correct)" << std::endl;
	for (unsigned int i=0; i<m; ++i)
		std::cout << z[i] - y[i] << std::endl;


	// y += d Ax ////////////////////////////////////
	// cblas
	FBlas::gemva(m, n, d, A, x, y);
	// direct
	for (unsigned int i=0; i<m; ++i) {
		FReal _z = FReal(0.);
		for (unsigned int j=0; j<n; ++j) 
			_z += A[j*m + i] * x[j];
		z[i] += _z * d;
	}
	// compare
	std::cout << "\ny += d Ax (zeros are correct)" << std::endl;
	for (unsigned int i=0; i<m; ++i)
		std::cout << z[i] - y[i] << std::endl;


	// y += d A^Tx ////////////////////////////////////
	// cblas
	FBlas::gemtva(m, n, d, A, x, y);
	// direct
	for (unsigned int i=0; i<m; ++i) {
		FReal _z = FReal(0.);
		for (unsigned int j=0; j<n; ++j) 
			_z += A[i*m + j] * x[j];
		z[i] += _z * d;
	}
	// compare
	std::cout << "\ny += d A^Tx (zeros are correct)" << std::endl;
	for (unsigned int i=0; i<m; ++i)
		std::cout << z[i] - y[i] << std::endl;



	delete [] A;
	delete [] x;
	delete [] y;
	delete [] z;

	
	return 0;
}