View difference between Paste ID: si2PmVZg and KCNEcvMi
SHOW: | | - or go back to the newest paste.
1
# RAM超過6GB,才可以run!! 
2
# 否則請更改矩陣大小(mat_size)或是列表的長度(length_list)
3
length_list = 30
4
mat_size = c(3000, 3000)
5
6
mat = matrix(rnorm(prod(mat_size)), mat_size[1], mat_size[2])
7
mat.list = lapply(1:length_list, function(i) mat)
8
9
res1 = function(){
10
	mat_sum = matrix(rep(0, prod(mat_size)), mat_size[1], mat_size[2])
11
	for (i in 1:length_list){
12
		mat_sum = mat_sum + mat.list[[i]]
13
	}
14
	mat_sum
15
}
16
res2 = function() Reduce('+', mat.list)
17
res3 = function() matrix(rowSums(matrix(unlist(mat.list), nrow=prod(mat_size))), nrow=mat_size[1])
18
res4 = function() matrix(colSums(do.call(rbind, lapply(mat.list, as.vector))), nrow=mat_size[1])
19
res5 = function() matrix(rowSums(do.call(cbind, lapply(mat.list, as.vector))), nrow=mat_size[1])
20
library(Rcpp)
21
library(RcppArmadillo)
22
library(inline)
23
sourceCpp(code = '
24
// [[Rcpp::depends(RcppArmadillo)]]
25
#include <RcppArmadillo.h>
26
#include <stdio.h>
27
using namespace Rcpp;
28
using namespace arma;
29
// [[Rcpp::export]]
30
mat list_mat_sum_f(List data_list){
31
  int n = data_list.size();
32
  SEXP mat1 = data_list[0];
33
  NumericMatrix mat1_rmat(mat1);
34
  mat result_mat(mat1_rmat.begin(), mat1_rmat.nrow(), mat1_rmat.ncol());
35
  for(int i = 1; i < n; i++)
36
  {
37
    SEXP tmp_m = data_list[i];
38
    NumericMatrix data_m(tmp_m);
39
	mat working_m(data_m.begin(), data_m.nrow(), data_m.ncol(), false);
40
    result_mat += working_m;
41
  }
42
  return result_mat;
43
}')
44
res6 = function() list_mat_sum_f(mat.list)
45
46
# cmpfun
47
library(compiler)
48
res1_cmp = cmpfun(res1)
49
res2_cmp = cmpfun(res2)
50
res3_cmp = cmpfun(res3)
51
res4_cmp = cmpfun(res4)
52
res5_cmp = cmpfun(res5)
53
54
all.equal(tmp <- res1(), res2())
55
# TRUE
56
all.equal(tmp, res3())
57
# TRUE
58
all.equal(tmp, res4())
59
# TRUE
60
all.equal(tmp, res5())
61
# TRUE
62
all.equal(tmp, res6())
63
# TRUE
64
65
library(rbenchmark)
66
benchmark(res1(),res2(),res3(),res4(),res5(),res6(),res1_cmp(),res2_cmp(),res3_cmp(),res4_cmp(),res5_cmp(),replications = 20, columns=c('test', 'replications', 'elapsed','relative', 'user.self'),  order='relative')
67
68
#          test replications elapsed relative user.self
69
# 6      res6()           20    6.49    1.000      6.13
70
# 2      res2()           20   15.58    2.401      8.45
71
# 8  res2_cmp()           20   16.19    2.495      8.61
72
# 1      res1()           20   18.44    2.841     10.61
73
# 7  res1_cmp()           20   19.33    2.978     11.15
74
# 5      res5()           20   56.60    8.721     40.24
75
# 11 res5_cmp()           20   56.79    8.750     40.36
76
# 4      res4()           20   70.48   10.860     54.63
77
# 10 res4_cmp()           20   70.53   10.867     55.17
78
# 3      res3()           20   92.87   14.310     76.21
79
# 9  res3_cmp()           20   93.74   14.444     77.19