patch.cc

patch.cc

diff --git a/src/shogun/distributions/Gaussian.cpp b/src/shogun/distributions/Gaussian.cpp
index 4955cc3d2..7c77440e5 100644
--- a/src/shogun/distributions/Gaussian.cpp
+++ b/src/shogun/distributions/Gaussian.cpp
@@ -10,14 +10,14 @@
  */
 #include <shogun/lib/config.h>
 
-#ifdef HAVE_LAPACK
-
+#include <shogun/mathematics/linalg/LinalgNamespace.h>
 #include <shogun/distributions/Gaussian.h>
 #include <shogun/mathematics/Math.h>
 #include <shogun/base/Parameter.h>
 #include <shogun/mathematics/lapack.h>
 
 using namespace shogun;
+using namespace linalg;
 
 CGaussian::CGaussian() : CDistribution(), m_constant(0), m_d(), m_u(), m_mean(), m_cov_type(FULL)
 {
@@ -166,12 +166,12 @@ float64_t CGaussian::update_params_em(float64_t* alpha_k, int32_t len)
                switch (cov_type)
                {
                        case FULL:
-                               //cblas_dger(CblasRowMajor, num_dim, num_dim, alpha_k[j], v.vector, 1, v.vector,
-                               //                       1, (double*) cov_sum.matrix, num_dim);
-                               SGMatrix<float64_t> feature_martix = SGVector<float64_t>::convert_to_matrix(v,num_dim, 1, false);
-                                SGMatrix<float64_t> temp_martix = SGMatrix<float64_t>::matrix_multiply(feature_martix, feature_martix,
-                                                false, true, alpha_k[j]);
-                                linalg::add(cov_sum, temp_martix, cov_sum);
+#ifdef HAVE_LAPACK
+                               cblas_dger(CblasRowMajor, num_dim, num_dim, alpha_k[j], v.vector, 1, v.vector,
+                                                        1, (double*) cov_sum.matrix, num_dim);
+#else
+                               linalg::cblas_dger<float64_t>(alpha_k[j], v, v, cov_sum);
+#endif
                                break;
                        case DIAG:
                                for (int32_t k=0; k<num_dim; k++)
@@ -226,8 +226,8 @@ float64_t CGaussian::compute_log_PDF(SGVector<float64_t> point)
 {
        ASSERT(m_mean.vector && m_d.vector)
        ASSERT(point.vlen == m_mean.vlen)
-       float64_t* difference=SG_MALLOC(float64_t, m_mean.vlen);
-       sg_memcpy(difference, point.vector, sizeof(float64_t)*m_mean.vlen);
+       SGVector<float64_t> difference(m_mean.vlen);
+       sg_memcpy(difference.vector, point.vector, sizeof(float64_t)*m_mean.vlen);
 
        for (int32_t i = 0; i < m_mean.vlen; i++)
                difference[i] -= m_mean.vector[i];
@@ -236,14 +236,15 @@ float64_t CGaussian::compute_log_PDF(SGVector<float64_t> point)
 
        if (m_cov_type==FULL)
        {
-               float64_t* temp_holder=SG_MALLOC(float64_t, m_d.vlen);
-               cblas_dgemv(CblasRowMajor, CblasNoTrans, m_d.vlen, m_d.vlen,
-                                       1, m_u.matrix, m_d.vlen, difference, 1, 0, temp_holder, 1);
-
+               SGVector<float64_t> temp_holder(m_d.vlen);
+               temp_holder.zero();
+               //cblas_dgemv(CblasRowMajor, CblasNoTrans, m_d.vlen, m_d.vlen,
+               //                      1, m_u.matrix, m_d.vlen, difference, 1, 0, temp_holder, 1);
+               linalg::cblas_dgemv<float64_t>(1.0, m_u, false, difference, 0.0, temp_holder);
+               
                for (int32_t i=0; i<m_d.vlen; i++)
                        answer+=temp_holder[i]*temp_holder[i]/m_d.vector[i];
 
-               SG_FREE(temp_holder);
        }
        else if (m_cov_type==DIAG)
        {
@@ -256,7 +257,6 @@ float64_t CGaussian::compute_log_PDF(SGVector<float64_t> point)
                        answer+=difference[i]*difference[i]/m_d.vector[0];
        }
 
-       SG_FREE(difference);
 
        return -0.5*answer;
 }
@@ -426,5 +426,3 @@ CGaussian* CGaussian::obtain_from_generic(CDistribution* distribution)
        SG_REF(casted);
        return casted;
 }
-
-#endif // HAVE_LAPACK
diff --git a/src/shogun/mathematics/linalg/LinalgNamespace.h b/src/shogun/mathematics/linalg/LinalgNamespace.h
index 374c3832d..9c39291ba 100644
--- a/src/shogun/mathematics/linalg/LinalgNamespace.h
+++ b/src/shogun/mathematics/linalg/LinalgNamespace.h
@@ -349,15 +349,46 @@ Container<T> add(Container<T>& a, Container<T>& b, T alpha=1, T beta=1)
        return result;
 }
 
+/**
+ * Performs the operation A = alpha * x * y' + A
+ *
+ * @param scaling factor for vector x
+ * @param vector x
+ * @param vector y
+ * @param Martix A
+ */
+template <typename T>
+void cblas_dger(T alpha, const SGVector<T> x, const SGVector<T> y, SGMatrix<T>& A)
+{
+       auto x_martix = SGVector<T>::convert_to_matrix(x, A.num_rows, 1, false);
+       auto y_martix = SGVector<T>::convert_to_matrix(y, A.num_cols, 1, false);
+
+       auto temp_martix = SGMatrix<T>::matrix_multiply(x_martix, y_martix, false, true, alpha);
+       add(A, temp_martix, A);
+}
+
+
+/**
+ * Performs the operation y = \alpha ax + \beta y
+ * This function multiplies a * x (after transposing a, if needed) 
+ * and multiplies the resulting matrix by alpha. It then multiplies vector y by beta. 
+ * It stores the sum of these two products in vector y
+ *
+ * @param scaling factor for vector ax
+ * @param vector a
+ * @param Whether to transpose matrix a
+ * @param vector x
+ * @param scaling factor for vector y
+ * @param vector y
+ */
 template <typename T>
+void cblas_dgemv(T alpha, SGMatrix<T> a, bool transpose, SGVector<T> x, T beta, SGVector<T>& y)
 {
+       auto temp_vector = matrix_prod(a, x, false);
+       add(y, temp_vector, y, alpha, beta);
 }
 
+
 /**
  * Compute the cholesky decomposition \f$A = L L^{*}\f$ or \f$A = U^{*} U\f$
  * of a Hermitian positive definite matrix