提交: a151270c4cbebcc9ddea54b99445edfa16260f85 - MolDS (git)

Commit MetaInfo

• 修订版: a151270c4cbebcc9ddea54b99445edfa16260f85 (tree)
• 时间: 2012-12-29 14:23:50
• 作者: Katsuhiko Nishimra <ktns.87@gmai...>
• Commiter: Katsuhiko Nishimra

Log Message

Integrate Blas/Lapack_GNU.cpp into Blas/Lapack.cpp. #30409

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

更改概述

• modified: src/Makefile_GNU
• modified: src/wrappers/Blas.cpp
• modified: src/wrappers/Lapack.cpp

差异

--- a/src/Makefile_GNU
+++ b/src/Makefile_GNU
@@ -28,9 +28,9 @@ OPENBLAS_LIBS = -lopenblas
 EXENAME = MolDS.out
 DEPFILE = obj/objfile.dep
 
-ALL_CPP_FILES = Main.cpp base/MolDSException.cpp wrappers/Blas_GNU.cpp wrappers/Lapack_GNU.cpp base/Utilities.cpp base/Enums.cpp base/MathUtilities.cpp base/MallocerFreer.cpp base/EularAngle.cpp base/Parameters.cpp base/atoms/Atom.cpp base/atoms/Hatom.cpp base/atoms/Liatom.cpp base/atoms/Catom.cpp base/atoms/Natom.cpp base/atoms/Oatom.cpp base/atoms/Satom.cpp base/factories/AtomFactory.cpp base/Molecule.cpp base/InputParser.cpp base/GTOExpansionSTO.cpp base/RealSphericalHarmonicsIndex.cpp base/loggers/MOLogger.cpp base/loggers/DensityLogger.cpp base/loggers/HoleDensityLogger.cpp base/loggers/ParticleDensityLogger.cpp cndo/Cndo2.cpp indo/Indo.cpp zindo/ZindoS.cpp mndo/Mndo.cpp am1/Am1.cpp am1/Am1D.cpp pm3/Pm3.cpp pm3/Pm3D.cpp pm3/Pm3Pddg.cpp base/factories/ElectronicStructureFactory.cpp md/MD.cpp mc/MC.cpp rpmd/RPMD.cpp nasco/NASCO.cpp optimization/Optimizer.cpp optimization/ConjugateGradient.cpp optimization/SteepestDescent.cpp optimization/BFGS.cpp base/factories/OptimizerFactory.cpp  base/MolDS.cpp
+ALL_CPP_FILES = Main.cpp base/MolDSException.cpp wrappers/Blas.cpp wrappers/Lapack.cpp base/Utilities.cpp base/Enums.cpp base/MathUtilities.cpp base/MallocerFreer.cpp base/EularAngle.cpp base/Parameters.cpp base/atoms/Atom.cpp base/atoms/Hatom.cpp base/atoms/Liatom.cpp base/atoms/Catom.cpp base/atoms/Natom.cpp base/atoms/Oatom.cpp base/atoms/Satom.cpp base/factories/AtomFactory.cpp base/Molecule.cpp base/InputParser.cpp base/GTOExpansionSTO.cpp base/RealSphericalHarmonicsIndex.cpp base/loggers/MOLogger.cpp base/loggers/DensityLogger.cpp base/loggers/HoleDensityLogger.cpp base/loggers/ParticleDensityLogger.cpp cndo/Cndo2.cpp indo/Indo.cpp zindo/ZindoS.cpp mndo/Mndo.cpp am1/Am1.cpp am1/Am1D.cpp pm3/Pm3.cpp pm3/Pm3D.cpp pm3/Pm3Pddg.cpp base/factories/ElectronicStructureFactory.cpp md/MD.cpp mc/MC.cpp rpmd/RPMD.cpp nasco/NASCO.cpp optimization/Optimizer.cpp optimization/ConjugateGradient.cpp optimization/SteepestDescent.cpp optimization/BFGS.cpp base/factories/OptimizerFactory.cpp  base/MolDS.cpp
 ALL_HEAD_FILES = base/PrintController.h base/MolDSException.h base/Uncopyable.h wrappers/Blas.h wrappers/Lapack.h base/Utilities.h base/Enums.h base/MathUtilities.h base/MallocerFreer.h base/EularAngle.h base/Parameters.h base/atoms/Atom.h base/atoms/Hatom.h base/atoms/Liatom.h base/atoms/Catom.h base/atoms/Natom.h base/atoms/Oatom.h base/atoms/Satom.h base/factories/AtomFactory.h base/Molecule.h base/InputParser.h base/GTOExpansionSTO.h base/RealSphericalHarmonicsIndex.h base/loggers/MOLogger.h base/loggers/DensityLogger.h base/loggers/HoleDensityLogger.h base/loggers/ParticleDensityLogger.h base/ElectronicStructure.h cndo/Cndo2.h cndo/ReducedOverlapAOsParameters.h indo/Indo.h zindo/ZindoS.h mndo/Mndo.h am1/Am1.h am1/Am1D.h pm3/Pm3.h pm3/Pm3D.h pm3/Pm3Pddg.h base/factories/ElectronicStructureFactory.h md/MD.h mc/MC.h rpmd/RPMD.h nasco/NASCO.h optimization/Optimizer.h optimization/ConjugateGradient.h optimization/SteepestDescent.h optimization/BFGS.h base/factories/OptimizerFactory.h base/MolDS.h
-ALL_OBJ_FILES = obj/Main.o obj/MolDSException.o obj/Blas_GNU.o obj/Lapack_GNU.o obj/Utilities.o obj/Enums.o obj/MathUtilities.o obj/MallocerFreer.o obj/EularAngle.o obj/Parameters.o obj/Atom.o obj/Hatom.o obj/Liatom.o obj/Catom.o obj/Natom.o obj/Oatom.o obj/Satom.o obj/AtomFactory.o obj/Molecule.o obj/InputParser.o obj/GTOExpansionSTO.o obj/RealSphericalHarmonicsIndex.o obj/MOLogger.o obj/DensityLogger.o obj/HoleDensityLogger.o obj/ParticleDensityLogger.o obj/Cndo2.o obj/Indo.o obj/ZindoS.o obj/Mndo.o obj/Am1.o obj/Am1D.o obj/Pm3.o obj/Pm3D.o obj/Pm3Pddg.o obj/ElectronicStructureFactory.o obj/MD.o obj/MC.o obj/RPMD.o obj/NASCO.o obj/Optimizer.o obj/ConjugateGradient.o obj/SteepestDescent.o obj/BFGS.o obj/OptimizerFactory.o obj/MolDS.o
+ALL_OBJ_FILES = obj/Main.o obj/MolDSException.o obj/Blas.o obj/Lapack.o obj/Utilities.o obj/Enums.o obj/MathUtilities.o obj/MallocerFreer.o obj/EularAngle.o obj/Parameters.o obj/Atom.o obj/Hatom.o obj/Liatom.o obj/Catom.o obj/Natom.o obj/Oatom.o obj/Satom.o obj/AtomFactory.o obj/Molecule.o obj/InputParser.o obj/GTOExpansionSTO.o obj/RealSphericalHarmonicsIndex.o obj/MOLogger.o obj/DensityLogger.o obj/HoleDensityLogger.o obj/ParticleDensityLogger.o obj/Cndo2.o obj/Indo.o obj/ZindoS.o obj/Mndo.o obj/Am1.o obj/Am1D.o obj/Pm3.o obj/Pm3D.o obj/Pm3Pddg.o obj/ElectronicStructureFactory.o obj/MD.o obj/MC.o obj/RPMD.o obj/NASCO.o obj/Optimizer.o obj/ConjugateGradient.o obj/SteepestDescent.o obj/BFGS.o obj/OptimizerFactory.o obj/MolDS.o
 
 $(EXENAME): $(DEPFILE) $(ALL_OBJ_FILES)
 	$(CC) -o $@ -Wl,-rpath=$(OPENBLAS_LIB_DIR) $(ALL_OBJ_FILES) -L$(OPENBLAS_LIB_DIR) $(LIBS) $(OPENBLAS_LIBS)

--- a/src/wrappers/Blas.cpp
+++ b/src/wrappers/Blas.cpp
@@ -25,7 +25,11 @@
 #include<string>
 #include<stdexcept>
 #include<boost/format.hpp>
+#ifdef __INTEL_COMPILER
 #include"mkl.h"
+#else
+#include"cblas.h"
+#endif
 #include"../base/PrintController.h"
 #include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
@@ -75,7 +79,12 @@ void Blas::Dcopy(molds_blas_int n,
 void Blas::Dcopy(molds_blas_int n,
                  double const* vectorX, molds_blas_int incrementX,
                  double*       vectorY, molds_blas_int incrementY) const{
+#ifdef __INTEL_COMPILER
    dcopy(&n, vectorX, &incrementX, vectorY, &incrementY);
+#else
+   double* x = const_cast<double*>(&vectorX[0]);
+   cblas_dcopy(n, x, incrementX, vectorY, incrementY);
+#endif
 }
 
 // vectorY = alpha*vectorX + vectorY
@@ -95,7 +104,12 @@ void Blas::Daxpy(molds_blas_int n, double alpha,
 void Blas::Daxpy(molds_blas_int n, double alpha,
            double const* vectorX, molds_blas_int incrementX,
            double*       vectorY, molds_blas_int incrementY) const{
+#ifdef __INTEL_COMPILER
    daxpy(&n, &alpha, vectorX, &incrementX, vectorY, &incrementY);
+#else
+   double* x = const_cast<double*>(&vectorX[0]);
+   cblas_daxpy(n, alpha, x, incrementX, vectorY, incrementY);
+#endif
 }
 
 // returns vectorX^T*vectorY
@@ -115,7 +129,13 @@ double Blas::Ddot(molds_blas_int n,
 double Blas::Ddot(molds_blas_int n,
             double const* vectorX, molds_blas_int incrementX,
             double const* vectorY, molds_blas_int incrementY)const{
+#ifdef __INTEL_COMPILER
    return ddot(&n, vectorX, &incrementX, vectorY, &incrementY);
+#else
+   double* x=const_cast<double*>(vectorX),
+         * y=const_cast<double*>(vectorY);
+   return cblas_ddot(n, x, incrementX, y, incrementY);
+#endif
 }
 
 // vectorY = matrixA*vectorX
@@ -147,6 +167,7 @@ void Blas::Dgemv(bool isColumnMajorMatrixA,
                  double beta,
                  double* vectorY,
                  molds_blas_int incrementY) const{
+#ifdef __INTEL_COMPILER
    double const* a = &matrixA[0][0];
    char transA;
    if(isColumnMajorMatrixA){
@@ -158,6 +179,20 @@ void Blas::Dgemv(bool isColumnMajorMatrixA,
    }
    molds_blas_int lda = m;
    dgemv(&transA, &m, &n, &alpha, a, &lda, vectorX, &incrementX, &beta, vectorY, &incrementY);
+#else
+   double* a = const_cast<double*>(&matrixA[0][0]);
+   double* x = const_cast<double*>(&vectorX[0]);
+   CBLAS_TRANSPOSE transA;
+   if(isColumnMajorMatrixA){
+      transA = CblasNoTrans;
+   }
+   else{
+      transA = CblasTrans;
+      swap(m,n);
+   }
+   int lda = m;
+   cblas_dgemv(CblasColMajor, transA, m, n, alpha, a, lda, x, incrementX, beta, vectorY, incrementY);
+#endif
 }
 
 // vectorY = matrixA*vectorX
@@ -183,10 +218,18 @@ void Blas::Dsymv(molds_blas_int n, double alpha,
            double const* vectorX, molds_blas_int incrementX,
            double beta,
            double*       vectorY, molds_blas_int incrementY) const{
+#ifdef __INTEL_COMPILER
    double const* a = &matrixA[0][0];
    char uploA='L';
    molds_blas_int lda = n;
    dsymv(&uploA, &n, &alpha, a, &lda, vectorX, &incrementX, &beta, vectorY, &incrementY);
+#else
+   double* a = const_cast<double*>(&matrixA[0][0]);
+   double* x = const_cast<double*>(&vectorX[0]);
+   CBLAS_UPLO uploA=CblasUpper;
+   int lda = n;
+   cblas_dsymv(CblasRowMajor, uploA, n, alpha, a, lda, x, incrementX, beta, vectorY, incrementY);
+#endif
 }
 
 // matrixA = alpha*vectorX*vectorX^T + matrixA
@@ -203,9 +246,16 @@ void Blas::Dsyr(molds_blas_int n, double alpha,
           double const* vectorX, molds_blas_int incrementX,
           double ** matrixA)const{
    double* a = &matrixA[0][0];
+#ifdef __INTEL_COMPILER
    char uploA='L';
    molds_blas_int lda = n;
    dsyr(&uploA, &n, &alpha, vectorX, &incrementX, a, &lda);
+#else
+   double* x = const_cast<double*>(&vectorX[0]);
+   CBLAS_UPLO uploA=CblasUpper;
+   int lda = n;
+   cblas_dsyr(CblasRowMajor, uploA, n, alpha, x, incrementX, a, lda);
+#endif
 #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++){
@@ -241,6 +291,7 @@ void Blas::Dgemm(bool isColumnMajorMatrixA,
                  double const* const* matrixB,
                  double beta,
                  double** matrixC) const{
+#ifdef __INTEL_COMPILER
    double const* a = &matrixA[0][0];
    double const* b = &matrixB[0][0];
    double*       c = &matrixC[0][0];
@@ -266,9 +317,40 @@ void Blas::Dgemm(bool isColumnMajorMatrixA,
       transB = 'T'; //kb=k
       ldb = n;
    }
+#else
+   double* a = const_cast<double*>(&matrixA[0][0]);
+   double* b = const_cast<double*>(&matrixB[0][0]);
+   double*       c = &matrixC[0][0];
+
+   int lda;
+   CBLAS_TRANSPOSE transA;
+   if(isColumnMajorMatrixA){
+      transA = CblasNoTrans;
+      lda = m;
+   }
+   else{
+      transA = CblasTrans;
+      lda = k;
+   }
+
+   int ldb;
+   CBLAS_TRANSPOSE transB;
+   if(isColumnMajorMatrixB){
+      transB = CblasNoTrans;
+      ldb = k;
+   }
+   else{
+      transB = CblasTrans;
+      ldb = n;
+   }
+#endif
 
    double* tmpC;
+#ifdef __INTEL_COMPILER
    tmpC = (double*)mkl_malloc( sizeof(double)*m*n, 16 );
+#else
+   tmpC = (double*)malloc( sizeof(double)*m*n);
+#endif
    for(molds_blas_int i=0; i<m; i++){
       for(molds_blas_int j=0; j<n; j++){
          tmpC[i+j*m] = matrixC[i][j];
@@ -276,13 +358,21 @@ void Blas::Dgemm(bool isColumnMajorMatrixA,
    }
    molds_blas_int ldc = m;
    //call blas
+#ifdef __INTEL_COMPILER
    dgemm(&transA, &transB, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, tmpC, &ldc);
+#else
+   cblas_dgemm(CblasColMajor, transA, transB, m, n, k, alpha, a, lda, b, ldb, beta, tmpC, ldc);
+#endif
    for(molds_blas_int i=0; i<m; i++){
       for(molds_blas_int j=0; j<n; j++){
          matrixC[i][j] = tmpC[i+j*m];
       }
    }
+#ifdef __INTEL_COMPILER
    mkl_free(tmpC);
+#else
+   free(tmpC);
+#endif
 }
 
 }

--- a/src/wrappers/Lapack.cpp
+++ b/src/wrappers/Lapack.cpp
@@ -24,11 +24,30 @@
 #include<string>
 #include<stdexcept>
 #include<boost/format.hpp>
+#ifdef __INTEL_COMPILER
 #include"mkl.h"
+#else
+#if ( __WORDSIZE == 32 )
+#else
+#define HAVE_LAPACK_CONFIG_H
+#define LAPACK_ILP64
+#endif
+#include"lapacke.h"
+#endif
 #include"../base/PrintController.h"
 #include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
 #include"Lapack.h"
+
+#ifdef __INTEL_COMPILER
+#define MOLDS_LAPACK_malloc(a,b) mkl_malloc(a,b)
+#define MOLDS_LAPACK_free(a) mkl_free(a)
+#else
+#define MOLDS_LAPACK_malloc(a,b) malloc(a)
+#define MOLDS_LAPACK_free(a) free(a)
+#endif
+
+
 using namespace std;
 using namespace MolDS_base;
 
@@ -121,10 +140,10 @@ molds_lapack_int Lapack::Dsyevd(double** matrix, double* eigenValues, molds_lapa
    }
 
    // malloc
-   work = (double*)mkl_malloc( sizeof(double)*lwork, 16 );
-   iwork = (molds_lapack_int*)mkl_malloc( sizeof(molds_lapack_int)*liwork, 16 );
-   convertedMatrix = (double*)mkl_malloc( sizeof(double)*size*size, 16 );
-   tempEigenValues = (double*)mkl_malloc( sizeof(double)*size, 16 );
+   work = (double*)MOLDS_LAPACK_malloc( sizeof(double)*lwork, 16 );
+   iwork = (molds_lapack_int*)MOLDS_LAPACK_malloc( sizeof(molds_lapack_int)*liwork, 16 );
+   convertedMatrix = (double*)MOLDS_LAPACK_malloc( sizeof(double)*size*size, 16 );
+   tempEigenValues = (double*)MOLDS_LAPACK_malloc( sizeof(double)*size, 16 );
 
    for(molds_lapack_int i = 0; i < size; i++){
       for(molds_lapack_int j = i; j < size; j++){
@@ -133,7 +152,11 @@ molds_lapack_int Lapack::Dsyevd(double** matrix, double* eigenValues, molds_lapa
    }
 
    // call Lapack
+#ifdef __INTEL_COMPILER
    dsyevd(&job, &uplo, &size, convertedMatrix, &lda, tempEigenValues, work, &lwork, iwork, &liwork, &info);
+#else
+   info = LAPACKE_dsyevd_work(LAPACK_COL_MAJOR, job, uplo, size, convertedMatrix, lda, tempEigenValues, work, lwork, iwork, liwork);
+#endif
 
    for(molds_lapack_int i = 0; i < size; i++){
       for(molds_lapack_int j = 0; j < size; j++){
@@ -160,10 +183,10 @@ molds_lapack_int Lapack::Dsyevd(double** matrix, double* eigenValues, molds_lapa
    //this->OutputLog(boost::format("size=%d lwork=%d liwork=%d k=%d info=%d\n") % size % lwork % liwork % k % info);
 
    // free
-   mkl_free(work);
-   mkl_free(iwork);
-   mkl_free(convertedMatrix);
-   mkl_free(tempEigenValues);
+   MOLDS_LAPACK_free(work);
+   MOLDS_LAPACK_free(iwork);
+   MOLDS_LAPACK_free(convertedMatrix);
+   MOLDS_LAPACK_free(tempEigenValues);
   
    if(info != 0){
       stringstream ss;
@@ -199,9 +222,9 @@ molds_lapack_int Lapack::Dsysv(double const* const* matrix, double* b, molds_lap
    }
 
    // malloc
-   ipiv = (molds_lapack_int*)mkl_malloc( sizeof(molds_lapack_int)*2*size, 16 );
-   convertedMatrix = (double*)mkl_malloc( sizeof(double)*size*size, 16 );
-   tempB = (double*)mkl_malloc( sizeof(double)*size, 16 );
+   ipiv = (molds_lapack_int*)MOLDS_LAPACK_malloc( sizeof(molds_lapack_int)*2*size, 16 );
+   convertedMatrix = (double*)MOLDS_LAPACK_malloc( sizeof(double)*size*size, 16 );
+   tempB = (double*)MOLDS_LAPACK_malloc( sizeof(double)*size, 16 );
 
    for(molds_lapack_int i = 0; i < size; i++){
       for(molds_lapack_int j = i; j < size; j++){
@@ -218,24 +241,32 @@ molds_lapack_int Lapack::Dsysv(double const* const* matrix, double* b, molds_lap
    {
       lwork = -1;
       double tempWork[3]={0.0, 0.0, 0.0};
-      dsysv(&uplo, &size, &nrhs, convertedMatrix, &lda, ipiv, tempB, &ldb, tempWork, &lwork, &info);
+#ifdef __INTEL_COMPILER
+         dsysv(&uplo, &size, &nrhs, convertedMatrix, &lda, ipiv, tempB, &ldb, tempWork, &lwork, &info);
+#else
+         info = LAPACKE_dsysv_work(LAPACK_COL_MAJOR, uplo, size, nrhs, convertedMatrix, lda, ipiv, tempB, ldb, tempWork, lwork);
+#endif
       blockSize = tempWork[0]/size;
    }
    info = 0;
    lwork = blockSize*size;
-   work = (double*)mkl_malloc( sizeof(double)*lwork, 16 );
+   work = (double*)MOLDS_LAPACK_malloc( sizeof(double)*lwork, 16 );
 
    // call Lapack
+#ifdef __INTEL_COMPILER
    dsysv(&uplo, &size, &nrhs, convertedMatrix, &lda, ipiv, tempB, &ldb, work, &lwork, &info);
+#else
+   info = LAPACKE_dsysv_work(LAPACK_COL_MAJOR, uplo, size, nrhs, convertedMatrix, lda, ipiv, tempB, ldb, work, lwork);
+#endif
    for(molds_lapack_int i = 0; i < size; i++){
       b[i] = tempB[i];
    }
 
    // free
-   mkl_free(convertedMatrix);
-   mkl_free(ipiv);
-   mkl_free(work);
-   mkl_free(tempB);
+   MOLDS_LAPACK_free(convertedMatrix);
+   MOLDS_LAPACK_free(ipiv);
+   MOLDS_LAPACK_free(work);
+   MOLDS_LAPACK_free(tempB);
   
    if(info != 0){
       stringstream ss;
@@ -271,9 +302,9 @@ molds_lapack_int Lapack::Dgetrs(double const* const* matrix, double** b, molds_l
 
    try{
       // malloc
-      ipiv = (molds_lapack_int*)mkl_malloc( sizeof(molds_lapack_int)*2*size, 16 );
-      convertedMatrix = (double*)mkl_malloc( sizeof(double)*size*size, 16 );
-      convertedB = (double*)mkl_malloc( sizeof(double)*nrhs*size, 16 );
+      ipiv = (molds_lapack_int*)MOLDS_LAPACK_malloc( sizeof(molds_lapack_int)*2*size, 16 );
+      convertedMatrix = (double*)MOLDS_LAPACK_malloc( sizeof(double)*size*size, 16 );
+      convertedB = (double*)MOLDS_LAPACK_malloc( sizeof(double)*nrhs*size, 16 );
       for(molds_lapack_int i = 0; i < size; i++){
          for(molds_lapack_int j = 0; j < size; j++){
             convertedMatrix[i+j*size] = matrix[i][j];
@@ -285,7 +316,11 @@ molds_lapack_int Lapack::Dgetrs(double const* const* matrix, double** b, molds_l
          }
       }
       this->Dgetrf(convertedMatrix, ipiv, size, size);
+#ifdef __INTEL_COMPILER
       dgetrs(&trans, &size, &nrhs, convertedMatrix, &lda, ipiv, convertedB, &ldb, &info);
+#else
+      info = LAPACKE_dgetrs_work(LAPACK_COL_MAJOR, trans, size, nrhs, convertedMatrix, lda, ipiv, convertedB, ldb);
+#endif
       for(molds_lapack_int i = 0; i < nrhs; i++){
          for(molds_lapack_int j = 0; j < size; j++){
             b[i][j] = convertedB[j+i*size];
@@ -294,15 +329,15 @@ molds_lapack_int Lapack::Dgetrs(double const* const* matrix, double** b, molds_l
    }
    catch(MolDSException ex){
       // free
-      mkl_free(convertedMatrix);
-      mkl_free(convertedB);
-      mkl_free(ipiv);
+      MOLDS_LAPACK_free(convertedMatrix);
+      MOLDS_LAPACK_free(convertedB);
+      MOLDS_LAPACK_free(ipiv);
       throw ex;
    }
    // free
-   mkl_free(convertedMatrix);
-   mkl_free(convertedB);
-   mkl_free(ipiv);
+   MOLDS_LAPACK_free(convertedMatrix);
+   MOLDS_LAPACK_free(convertedB);
+   MOLDS_LAPACK_free(ipiv);
   
    if(info != 0){
       stringstream ss;
@@ -316,9 +351,9 @@ molds_lapack_int Lapack::Dgetrs(double const* const* matrix, double** b, molds_l
 // Argument "matrix" is sizeM*sizeN matrix.
 // Argument "matrix" will be LU-decomposed.
 molds_lapack_int Lapack::Dgetrf(double** matrix, molds_lapack_int sizeM, molds_lapack_int sizeN) const{
-   molds_lapack_int* ipiv = (molds_lapack_int*)mkl_malloc( sizeof(molds_lapack_int)*2*sizeM,        16 );
+   molds_lapack_int* ipiv = (molds_lapack_int*)MOLDS_LAPACK_malloc( sizeof(molds_lapack_int)*2*sizeM,        16 );
    this->Dgetrf(matrix, ipiv, sizeM, sizeN);
-   mkl_free(ipiv);
+   MOLDS_LAPACK_free(ipiv);
    molds_lapack_int info = 0;
    return info;
 }
@@ -326,7 +361,7 @@ molds_lapack_int Lapack::Dgetrf(double** matrix, molds_lapack_int sizeM, molds_l
 // Argument "matrix" is sizeM*sizeN matrix.
 // Argument "matrix" will be LU-decomposed.
 molds_lapack_int Lapack::Dgetrf(double** matrix, molds_lapack_int* ipiv, molds_lapack_int sizeM, molds_lapack_int sizeN) const{
-   double* convertedMatrix = (double*)mkl_malloc( sizeof(double)*sizeM*sizeN, 16 );
+   double* convertedMatrix = (double*)MOLDS_LAPACK_malloc( sizeof(double)*sizeM*sizeN, 16 );
    for(molds_lapack_int i=0; i<sizeM; i++){
       for(molds_lapack_int j=0; j<sizeN; j++){
          convertedMatrix[i+j*sizeM] = matrix[i][j];
@@ -338,7 +373,7 @@ molds_lapack_int Lapack::Dgetrf(double** matrix, molds_lapack_int* ipiv, molds_l
          matrix[i][j] = convertedMatrix[i+j*sizeM];
       }
    }
-   mkl_free(convertedMatrix);
+   MOLDS_LAPACK_free(convertedMatrix);
    molds_lapack_int info = 0;
    return info;
 }
@@ -348,7 +383,11 @@ molds_lapack_int Lapack::Dgetrf(double** matrix, molds_lapack_int* ipiv, molds_l
 molds_lapack_int Lapack::Dgetrf(double* matrix, molds_lapack_int* ipiv, molds_lapack_int sizeM, molds_lapack_int sizeN) const{
    molds_lapack_int info = 0;
    molds_lapack_int lda = sizeM;
+#ifdef __INTEL_COMPILER
    dgetrf(&sizeM, &sizeN, matrix, &lda, ipiv, &info);
+#else
+   info = LAPACKE_dgetrf_work(LAPACK_COL_MAJOR, sizeM, sizeN, matrix, lda, ipiv);
+#endif
    if(info != 0){
       stringstream ss;
       ss << errorMessageDgetrfInfo;