Commit: c264b8f014af8f5cd8a6bddef813dd40359066da - MolDS (git)

Commit MetaInfo

• Revisión: c264b8f014af8f5cd8a6bddef813dd40359066da (tree)
• Tiempo: 2013-01-23 19:46:12
• Autor: Katsuhiko Nishimra <ktns.87@gmai...>
• Commiter: Katsuhiko Nishimra

Log Message

Implement the wrapper function for dsyrk. #30562 #30566

git-svn-id: https://svn.sourceforge.jp/svnroot/molds/trunk@1247 1136aad2-a195-0410-b898-f5ea1d11b9d8

Cambiar Resumen

• modified: src/wrappers/Blas.cpp
• modified: src/wrappers/Blas.h

Diferencia incremental

--- a/src/wrappers/Blas.cpp
+++ b/src/wrappers/Blas.cpp
@@ -375,4 +375,46 @@ void Blas::Dgemm(bool isColumnMajorMatrixA,
 #endif
 }
 
+// matrixC = matrixA*matrixA^T
+//    matrixA: n*k-matrix
+//    matrixC: n*n-matrix,symmetric (Use the upper triangular part, and copy it to the lower part.)
+void Blas::Dsyrk(molds_blas_int n, molds_blas_int k,
+                 double const* const* matrixA,
+                 double ** matrixC)const{
+   bool isMatrixATransposed = false;
+   bool isLowerTriangularPartMatrixCUsed = false;
+   double alpha = 1.0 , beta = 0.0;
+   this->Dsyrk(n, k, isMatrixATransposed, isLowerTriangularPartMatrixCUsed, alpha, matrixA, beta, matrixC);
+}
+
+// matrixC = alpha*matrixA*matrixA^T + beta*matrixC (isMatrixATransposed==false)
+//  or
+// matrixC = alpha*matrixA^T*matrixA + beta*matrixC (isMatrixATransposed==true)
+//    matrixA: n*k-matrix (isMatrixATransposed==false) or k*n-matrix (isMatrixATransposed==true)
+//    matrixC: n*n-matrix,symmetric (Use the upper triangular part, and copy it to the lower part.)
+void Blas::Dsyrk(molds_blas_int n, molds_blas_int k,
+                 bool isMatrixATransposed,
+                 bool isLowerTriangularPartMatrixCUsed,
+                 double alpha, double const* const* matrixA,
+                 double beta,  double ** matrixC)const{
+   double* c = &matrixC[0][0];
+   double* a = const_cast<double*>(&matrixA[0][0]);
+   CBLAS_UPLO uploC= isLowerTriangularPartMatrixCUsed ? CblasLower : CblasUpper;
+   CBLAS_TRANSPOSE transA= isMatrixATransposed ? CblasTrans : CblasNoTrans;
+   molds_blas_int lda = &matrixA[1][0] - &matrixA[0][0];
+   molds_blas_int ldc = &matrixC[1][0] - &matrixC[0][0];
+   cblas_dsyrk(CblasRowMajor, uploC, transA, n, k, alpha, a, lda, beta, c, ldc);
+#pragma omp parallel for schedule(auto)
+   for(molds_blas_int i=0;i<n;i++){
+      for(molds_blas_int j=i+1;j<n;j++){
+         if(isLowerTriangularPartMatrixCUsed){
+            matrixC[i][j] = matrixC[j][i];
+         }
+         else{
+            matrixC[j][i] = matrixC[i][j];
+         }
+      }
+   }
+}
+
 }

--- a/src/wrappers/Blas.h
+++ b/src/wrappers/Blas.h
@@ -85,6 +85,14 @@ public:
               double const* const* matrixB,
               double beta,
               double**             matrixC) const;
+              void Dsyrk(molds_blas_int n, molds_blas_int k,
+                         double const *const* matrixA,
+                         double**             matrixC)const;
+              void Dsyrk(molds_blas_int n, molds_blas_int k,
+                         bool isMatrixATransposed,
+                         bool isLowerTriangularPartMatrixCUsed,
+                         double alpha, double const* const* matrixA,
+                         double beta,  double**             matrixC)const;
 private:
    Blas();
    ~Blas();