Revisión | 1d8136b2244e23760d755d77c85ff26785508a8e (tree) |
---|---|
Tiempo | 2013-01-22 20:38:52 |
Autor | Katsuhiko Nishimra <ktns.87@gmai...> |
Commiter | Katsuhiko Nishimra |
WIP on gdiisl: 7d2b4b1 tmp BFGS output SCF log
@@ -375,4 +375,46 @@ void Blas::Dgemm(bool isColumnMajorMatrixA, | ||
375 | 375 | #endif |
376 | 376 | } |
377 | 377 | |
378 | +// matrixC = matrixA*matrixA^T | |
379 | +// matrixA: n*k-matrix | |
380 | +// matrixC: n*n-matrix,symmetric (Use the upper triangular part, and copy it to the lower part.) | |
381 | +void Blas::Dsyrk(molds_blas_int n, molds_blas_int k, | |
382 | + double const** matrixA, | |
383 | + double ** matrixA)const{ | |
384 | + molds_blas_int lda = n; | |
385 | + molds_blas_int ldc = n; | |
386 | + bool isMatrixATransposed = false; | |
387 | + bool isLowerTriangularPartMatrixCUsed = false; | |
388 | + this->Dsyrk(n, k, isMatrixATransposed, isLowerTriangularPartMatrixCUsed, alpha, matrixA, lda, matrixC, ldc); | |
389 | +} | |
390 | + | |
391 | +// matrixC = alpha*matrixA*matrixA^T + beta*matrixC (isMatrixATransposed==false) | |
392 | +// or | |
393 | +// matrixC = alpha*matrixA^T*matrixA + beta*matrixC (isMatrixATransposed==true) | |
394 | +// matrixA: n*k-matrix (isMatrixATransposed==false) or k*n-matrix (isMatrixATransposed==true) | |
395 | +// matrixC: n*n-matrix,symmetric (Use the upper triangular part, and copy it to the lower part.) | |
396 | +void Blas::Dsyr(molds_blas_int n, molds_blas_int k, | |
397 | + bool isMatrixATransposed, | |
398 | + bool isLowerTriangularPartMatrixCUsed, | |
399 | + double alpha, double const* vectorX, molds_blas_int lda, | |
400 | + double beta, double ** matrixC, molds_blas_int ldc)const{ | |
401 | + double* c = &matrixC[0][0]; | |
402 | +#ifdef __INTEL_COMPILER | |
403 | + char uploA='L'; | |
404 | + molds_blas_int lda = n; | |
405 | + dsyr(&uploA, &n, &alpha, vectorX, &incrementX, a, &lda); | |
406 | +#else | |
407 | + double* x = const_cast<double*>(&vectorX[0]); | |
408 | + CBLAS_UPLO uploA=CblasUpper; | |
409 | + int lda = n; | |
410 | + cblas_dsyr(CblasRowMajor, uploA, n, alpha, x, incrementX, a, lda); | |
411 | +#endif | |
412 | +#pragma omp parallel for schedule(auto) | |
413 | + for(molds_blas_int i=0;i<n;i++){ | |
414 | + for(molds_blas_int j=i+1;j<n;j++){ | |
415 | + matrixA[j][i] = matrixA[i][j]; | |
416 | + } | |
417 | + } | |
418 | +} | |
419 | + | |
378 | 420 | } |