2016-07-14 32 views
0

Ich arbeite an C, mit GNU-Bibliothek für wissenschaftliche Datenverarbeitung. Im Wesentlichen müssen, ich das Äquivalent des folgenden MATLAB-Code tun:Elementweise Produkt zwischen einem Vektor und einer Matrix mit GNU Blas Subroutinen

x=x.*(A*x); 

wobei x eine gsl_vector ist, und A ein gsl_matrix.

verwaltet I (A * x) mit dem folgenden Befehl zu tun:

gsl_blas_dgemv(CblasNoTrans, 1.0, A, x, 1.0, res); 

wo res ist eine andere gsl_vector, die das Ergebnis speichert. Wenn die Matrix A die Größe m * m hat und der Vektor x die Größe m * 1 hat, hat der Vektor res die Größe m * 1.

Jetzt müssen wir noch das elementweise Produkt der Vektoren x und res (die Ergebnisse sollten ein Vektor sein). Leider stecke ich darauf und kann die Funktion nicht finden, die das tut.

Wenn mir jemand dabei helfen kann, wäre ich sehr dankbar. Außerdem weiß jemand, ob es eine bessere Dokumentation von GNU als https://www.gnu.org/software/gsl/manual/html_node/GSL-BLAS-Interface.html#GSL-BLAS-Interface gibt, die mich bisher verwirren.

Schließlich, würde ich in der Zeit Leistung verlieren, wenn ich diesen Schritt tun, indem Sie einfach eine for-Schleife (die Größe des Vektors ist um 11000 und dieser Schritt wird 500-5000 mal wiederholt)?

for (i = 0; i < m; i++) 
    gsl_vector_set(res, i, gsl_vector_get(x, i) * gsl_vector_get(res, i)); 

Vielen Dank!

Antwort

2

Die Funktion ist Sie implementieren können:

gsl_vector_mul(res, x) 

Ich habe Intels MKL verwendet, und ich mag die Dokumentation auf ihrer Website für diese BLAS-Routinen.

+0

Tatsächlich scheint das genau so zu funktionieren, wie ich es wollte. – TheRevanchist

0

Die For-Schleife ist in Ordnung, wenn GSL gut gestaltet ist. Zum Beispiel gsl_vector_set() und gsl_vector_get() können inlined sein. Sie können die Laufzeit mit gsl_blas_daxpy vergleichen. Die For-Schleife ist gut optimiert, wenn das Timing-Ergebnis ähnlich ist.

Auf der anderen Seite können Sie eine viel bessere Matrix Bibliothek Eigen versuchen wollen, mit denen Sie Ihren Betrieb mit dem Code ähnlich wie diese

x = x.array() * (A * x).array(); 
+0

Eigen funktioniert nicht in C, es funktioniert nur in C++. – TheRevanchist

+0

ok. Ich habe es verpasst. – kangshiyin