matprod.R 4.22 KB
 Radford Neal committed Nov 07, 2013 1 2 3 ``````# Test matrix multiplication, with %*%, crossprod, and tcrossprod, by # BLAS and 'matprod' routines. # `````` Radford Neal committed Sep 25, 2018 4 ``````# Added for pqR, 2013, 2018 Radford M. Neal. `````` Radford Neal committed Nov 07, 2013 5 6 7 8 9 10 11 12 13 `````` # Check matrix multiplication with various sizes of matrices, setting # matrix elements to random values that are integer multiples of 1/8 so # that floating point arithmetic will be exact (but accidental conversions # to integer will be detected). The crossprod and tcrossprod routines # are also checked when give one argument (producing a symmetric result). # # Returns the last (largest) result matrix. `````` Radford Neal committed Sep 25, 2018 14 15 16 17 18 19 20 ``````check_matprod <- function (print=TRUE) { if (print) { cat("\n") print(options()[c("mat_mult_with_BLAS","helpers_disable")]) } `````` Radford Neal committed Nov 07, 2013 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 `````` # Matrix multiply the hard way, to check results. matmult <- function (A,B) { n <- nrow(A) m <- ncol(B) k <- ncol(A) stopifnot(nrow(B)==k) C <- matrix(0,n,m) for (i in seq_len(n)) { for (j in seq_len(m)) { C[i,j] <- sum (A[i,] * B[,j]) } } C } `````` Radford Neal committed Sep 25, 2018 37 38 39 40 41 42 `````` # Do checks with given matrix sizes. check <- function (n, m, k) { A <- matrix (rgeom(n*k,0.1)/8, n, k) B <- matrix (rgeom(k*m,0.1)/8, k, m) C <- matmult(A,B) `````` Radford Neal committed Sep 25, 2018 43 44 `````` C1 <- matmult(A,B+1) C2 <- matmult(A,B)+2 `````` Radford Neal committed Sep 25, 2018 45 46 47 48 49 `````` At <- t(A) Bt <- t(B) stopifnot(identical(A%*%B,C)) stopifnot(identical(crossprod(At,B),C)) stopifnot(identical(tcrossprod(A,Bt),C)) `````` Radford Neal committed Sep 25, 2018 50 51 52 53 54 55 `````` stopifnot(identical(A%*%(B+1),C1)) stopifnot(identical(crossprod(At,B+1),C1)) stopifnot(identical(tcrossprod(A,Bt+1),C1)) stopifnot(identical(A%*%B+2,C2)) stopifnot(identical(crossprod(At,B)+2,C2)) stopifnot(identical(tcrossprod(A,Bt)+2,C2)) `````` Radford Neal committed Sep 25, 2018 56 57 58 59 60 61 62 63 `````` if (n==m) { D <- matmult(At,A) stopifnot(identical(crossprod(A),D)) stopifnot(identical(tcrossprod(At),D)) } C } `````` Radford Neal committed Nov 10, 2018 64 65 66 67 `````` # Set seed to get consistent results. set.seed(1) `````` Radford Neal committed Sep 25, 2018 68 `````` # Try various sizes systematically. `````` Radford Neal committed Nov 07, 2013 69 `````` `````` Radford Neal committed Nov 10, 2018 70 71 `````` cat("---- (0..11) (0..11) (0..11)\n") `````` 72 73 74 `````` for (n in 0..11) { for (m in 0..11) { for (k in 0..11) `````` Radford Neal committed Sep 25, 2018 75 `````` { C <- check(n,m,k) `````` Radford Neal committed Nov 07, 2013 76 77 78 79 `````` } } } `````` Radford Neal committed Sep 25, 2018 80 81 82 83 84 `````` if (print) { cat("\n") print(C) } `````` Radford Neal committed Nov 10, 2018 85 86 87 88 89 `````` cat("---- n (1 2 3 4 5) m\n") for (n in c(2,3,4,5,8,100,1000,2000,4000)) { for (m in c(2,3,4,5,8,100,1000,2000,4000)) { if (n*m <= 1000000) `````` 90 `````` { if (print) print(c(n,m)) `````` Radford Neal committed Nov 10, 2018 91 `````` for (k in c(1,2,3,4,5)) check(n,m,k) `````` 92 93 94 95 `````` } } } `````` Radford Neal committed Nov 10, 2018 96 97 98 99 100 101 102 103 104 105 106 107 108 `````` cat("---- (1 2 3 4 5) k m\n") for (k in c(2,3,4,5,8,100,1000,2000)) { for (m in c(2,3,4,5,8,100,1000,2000)) { if (k*m <= 1000000) { if (print) print(c(k,m)) for (n in c(1,2,3,4,5)) check(n,m,k) } } } cat("---- random\n") `````` Radford Neal committed Sep 25, 2018 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 `````` # Try various sizes randomly. Can increase loop count for more testing # by setting the environment variable R_MATPROD_TEST_COUNT (minimum 200). count <- as.integer(Sys.getenv("R_MATPROD_TEST_COUNT")) if (is.na(count) || count < 200) count <- 200 for (i in 1..count) { n <- if (runif(1)<0.6) 15*rpois(1,10) + rpois(1,10) else if (runif(1)>0.93) 50*rpois(1,12) + rpois(1,30) else rpois(1,5) m <- if (runif(1)<0.6) 15*rpois(1,10) + rpois(1,10) else if (runif(1)>0.93) 50*rpois(1,12) + rpois(1,30) else rpois(1,5) k <- if (runif(1)<0.6) 15*rpois(1,10) + rpois(1,10) else if (runif(1)>0.93) 50*rpois(1,12) + rpois(1,30) else rpois(1,5) if (n > 1200) n <- 1200 if (m > 1200) m <- 1200 if (k > 1200) k <- 1200 # Print n,m,k, but only for the first 200, so increasing the loop # repetitions won't invalidate saved output. if (print && i <= 200) print(c(n,m,k)) check(n,m,k) } `````` Radford Neal committed Nov 07, 2013 133 134 ``````} `````` Radford Neal committed Sep 25, 2018 135 136 ``````# Check matrix products using BLAS and using 'matprod' routines, with # or without helper threads enabled. `````` Radford Neal committed Nov 07, 2013 137 `````` `````` Radford Neal committed Sep 25, 2018 138 ``````sv <- options()[c("mat_mult_with_BLAS","helpers_disable")] `````` Radford Neal committed Nov 07, 2013 139 `````` `````` Radford Neal committed Sep 25, 2018 140 ``````options(mat_mult_with_BLAS=FALSE,helpers_disable=FALSE) `````` Radford Neal committed Nov 10, 2018 141 ``````cat("\nNot BLAS, Helpers not disabled\n\n") `````` Radford Neal committed Sep 25, 2018 142 ``````check_matprod() `````` Radford Neal committed Nov 07, 2013 143 `````` `````` Radford Neal committed Sep 25, 2018 144 ``````options(mat_mult_with_BLAS=FALSE,helpers_disable=TRUE) `````` Radford Neal committed Nov 10, 2018 145 ``````cat("\nNot BLAS, Helpers disabled\n\n") `````` Radford Neal committed Sep 25, 2018 146 147 148 149 ``````check_matprod() if (identical(Sys.getenv("R_MATPROD_TEST_BLAS"),"TRUE")) { `````` Radford Neal committed Nov 10, 2018 150 `````` cat("\nBLAS, Helpers not disabled\n\n") `````` Radford Neal committed Sep 25, 2018 151 152 153 `````` options(mat_mult_with_BLAS=TRUE,helpers_disable=FALSE) check_matprod(print=FALSE) `````` Radford Neal committed Nov 10, 2018 154 `````` cat("\nBLAS, Helpers disabled\n\n") `````` Radford Neal committed Sep 25, 2018 155 156 157 `````` options(mat_mult_with_BLAS=TRUE,helpers_disable=TRUE) check_matprod(print=FALSE) } `````` Radford Neal committed Nov 07, 2013 158 `````` `````` Radford Neal committed Sep 25, 2018 159 ``options(sv)``