Я попробовал ваш код на двухъядерном ноутбуке под управлением Debian Linux и R 3.5.1, связанных с OpenBLAS.Кроме того, я также использовал RcppArrayFire для выполнения этих вычислений с использованием графического процессора (ничего особенного, просто встроенная графика).Вот результаты теста:
Unit: microseconds
expr min lq mean median uq max neval cld
R10 55.488 90.142 117.9538 133.617 138.009 161.604 10 a
R100 698.883 711.488 1409.8286 773.352 901.339 5647.804 10 a
R1000 198612.971 213531.557 225713.6705 214993.916 226526.513 306675.313 10 d
Cpp10 16.157 31.478 38.3473 40.529 51.122 52.351 10 a
Cpp100 519.527 544.728 573.0099 570.956 610.985 613.400 10 a
Cpp1000 174121.655 184865.003 196224.3825 193142.715 207066.037 223900.830 10 c
AF10 700.805 790.203 1469.4980 1347.639 1410.717 3824.905 10 a
AF100 1376.737 1562.675 1818.8094 1898.797 1949.364 2222.889 10 a
AF1000 106398.431 109130.482 118817.6704 118612.679 123193.649 139474.579 10 b
При наименьшем размере (R10 и Cpp10) ваша машина быстрее моей.Но уже на R100 и Cpp100, в частности на R1000 и Cpp1000, я получаю более быстрое выполнение.Как отметил Дирк в комментариях, вы должны изучить оптимизированные / параллельные реализации BLAS.В Debian Linux (и производных, таких как Ubuntu) это так же просто, как
sudo apt-get install libopenblas-base
См. Также здесь .Теперь о результатах с использованием графического процессора: это довольно типично.Для более мелких матриц это хуже, чем как у основания R, так и у броненосца.Использование графического процессора приводит к значительным накладным расходам при перемещении данных между основной памятью и памятью графического процессора.Но для самого большого размера параллельное выполнение на графическом процессоре перевешивает эти издержки, и выигрыш в производительности довольно приятный.
Вот мой код для справки.Я позволил себе обновить inline::cxxfunction
до Rcpp::cppFunction
:
#Create the basic R implementation of the equation
flowCalc <- function(A,C){
B <- t(A) %*% C %*% A
Imp <- solve(B)
PTDF <- C %*% A %*% Imp
Out <- list(Imp = Imp, PTDF= PTDF)
return(Out)
}
# Create Armadillo function
Rcpp::cppFunction(depends = "RcppArmadillo", code = '
Rcpp::List flowCalcCpp(const arma::mat &Am, const arma::mat &Cm) {
arma::mat B = inv(trans(Am) * Cm * Am);
arma::mat PTDF = Cm * Am * B;
return Rcpp::List::create( Rcpp::Named("Imp") = B ,
Rcpp::Named("PTDF") = PTDF );
}')
# Create ArrayFire function
Rcpp::cppFunction(depends = "RcppArrayFire", code = '
Rcpp::List flowCalcAF(const RcppArrayFire::typed_array<f32> &A,
const RcppArrayFire::typed_array<f32> &C) {
af::array B = af::inverse(af::matmul(af::matmulTN(A, C), A));
af::array PTDF = af::matmul(af::matmul(C, A), B);
return Rcpp::List::create( Rcpp::Named("Imp") = B ,
Rcpp::Named("PTDF") = PTDF );
}')
library(igraph)
#Create a function to generate dummy data
MakeData <- function(n, edgep){#make an erdos renyi graph of size x
g <- erdos.renyi.game(n, edgep)
#convert to adjacency matrix
Adjmat <- as_adjacency_matrix(g, sparse = F)
#create random graph and mask the elements with not edge
Cmat <- matrix(rnorm(n*n), ncol = n)*Adjmat
##Output a list of the two matrices
list(A = Adjmat, C = Cmat)
}
#generate dummy data
set.seed(133)
Data10 <- MakeData(10, 1/1)
Data100 <- MakeData(100, 1/10)
Data1000 <- MakeData(1000, 1/100)
#Compare results
microbenchmark::microbenchmark(
R10 = flowCalc(Data10$A, Data10$C),
R100 = flowCalc(Data100$A, Data100$C),
R1000 = flowCalc(Data1000$A, Data1000$C),
Cpp10 = flowCalcCpp(Data10$A, Data10$C),
Cpp100 = flowCalcCpp(Data100$A, Data100$C),
Cpp1000 = flowCalcCpp(Data1000$A, Data1000$C),
AF10 = flowCalcAF(Data10$A, Data10$C),
AF100 = flowCalcAF(Data100$A, Data100$C),
AF1000 = flowCalcAF(Data1000$A, Data1000$C),
times = 10)