matprod.R 4.22 KB
Newer Older
Radford Neal's avatar
Radford Neal committed
1 2 3
# Test matrix multiplication, with %*%, crossprod, and tcrossprod, by
# BLAS and 'matprod' routines.
#
4
# Added for pqR, 2013, 2018 Radford M. Neal.
Radford Neal's avatar
Radford Neal committed
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.

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's avatar
Radford Neal committed
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
  }

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)
43 44
    C1 <- matmult(A,B+1)
    C2 <- matmult(A,B)+2
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))
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))
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's avatar
Radford Neal committed
64 65 66 67
  # Set seed to get consistent results.

  set.seed(1)

68
  # Try various sizes systematically.
Radford Neal's avatar
Radford Neal committed
69

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)
75
      { C <- check(n,m,k)
Radford Neal's avatar
Radford Neal committed
76 77 78 79
      }
    }
  }

80 81 82 83 84
  if (print)
  { cat("\n")
    print(C)
  }

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))
91
        for (k in c(1,2,3,4,5)) check(n,m,k)
92 93 94 95
      }
    }
  }

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")

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's avatar
Radford Neal committed
133 134
}

135 136
# Check matrix products using BLAS and using 'matprod' routines, with
# or without helper threads enabled.
Radford Neal's avatar
Radford Neal committed
137

138
sv <- options()[c("mat_mult_with_BLAS","helpers_disable")]
Radford Neal's avatar
Radford Neal committed
139

140
options(mat_mult_with_BLAS=FALSE,helpers_disable=FALSE)
141
cat("\nNot BLAS, Helpers not disabled\n\n")
142
check_matprod()
Radford Neal's avatar
Radford Neal committed
143

144
options(mat_mult_with_BLAS=FALSE,helpers_disable=TRUE)
145
cat("\nNot BLAS, Helpers disabled\n\n")
146 147 148 149
check_matprod()

if (identical(Sys.getenv("R_MATPROD_TEST_BLAS"),"TRUE")) {

150
    cat("\nBLAS, Helpers not disabled\n\n")
151 152 153
    options(mat_mult_with_BLAS=TRUE,helpers_disable=FALSE)
    check_matprod(print=FALSE)

154
    cat("\nBLAS, Helpers disabled\n\n")
155 156 157
    options(mat_mult_with_BLAS=TRUE,helpers_disable=TRUE)
    check_matprod(print=FALSE)
}
Radford Neal's avatar
Radford Neal committed
158

159
options(sv)