1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
|
## PR#4582 %*% with NAs
stopifnot(is.na(NA %*% 0), is.na(0 %*% NA))
## depended on the BLAS in use.
## found from fallback test in slam 0.1-15
## most likely indicates an inaedquate BLAS.
x <- matrix(c(1, 0, NA, 1), 2, 2)
y <- matrix(c(1, 0, 0, 2, 1, 0), 3, 2)
(z <- tcrossprod(x, y))
stopifnot(identical(z, x %*% t(y)))
stopifnot(is.nan(log(0) %*% 0))
## depended on the BLAS in use: some (including the reference BLAS)
## had z[1,3] == 0 and log(0) %*% 0 as as.matrix(0).
## matrix products
for(mopt in c("default","internal","default.simd")) {
# matprod="blas" is excluded because some tests fail due to issues
# in NaN/Inf propagation even in Rblas
options(matprod=mopt)
m <- matrix(c(1,2,3,4), ncol=2)
v <- c(11,12)
rv <- v; dim(rv) <- c(1,2)
cv <- v; dim(cv) <- c(2,1)
Cm <- m+0*1i; # cast to complex keeping dimensions
Cv <- v+0*1i;
Ccv <- cv+0*1i;
Crv <- rv+0*1i;
cprod <- function(rres, cres, expected) {
stopifnot(identical(rres, expected))
stopifnot(identical(Re(cres), expected))
}
cprod(m %*% m, Cm %*% Cm, matrix(c(7,10,15,22), 2, 2) )
cprod(m %*% cv, Cm %*% Ccv, matrix(c(47,70), 2, 1) )
cprod(m %*% v, Cm %*% Cv, matrix(c(47,70), 2, 1) )
cprod(rv %*% m, Crv %*% Cm, matrix(c(35,81), 1, 2) )
cprod(v %*% m, Cv %*% Cm, matrix(c(35,81), 1, 2) )
cprod(rv %*% cv, Crv %*% Ccv, matrix(265,1,1) )
cprod(cv %*% rv, Ccv %*% Crv, matrix(c(121,132,132,144), 2, 2) )
cprod(v %*% v, Cv %*% Cv, matrix(265,1,1) )
cprod(crossprod(m, m), crossprod(Cm, Cm), matrix(c(5,11,11,25), 2, 2) )
cprod(crossprod(m), crossprod(Cm), matrix(c(5,11,11,25), 2, 2) )
cprod(crossprod(m, cv), crossprod(Cm, Ccv), matrix(c(35,81), 2, 1) )
cprod(crossprod(m, v), crossprod(Cm, Cv), matrix(c(35,81), 2, 1) )
cprod(crossprod(cv, m), crossprod(Ccv, Cm), matrix(c(35,81), 1, 2) )
cprod(crossprod(v, m), crossprod(Cv, Cm), matrix(c(35,81), 1, 2) )
cprod(crossprod(cv, cv), crossprod(Ccv, Ccv), matrix(265,1,1) )
cprod(crossprod(v, v), crossprod(Cv, Cv), matrix(265,1,1) )
cprod(crossprod(rv, rv), crossprod(Crv, Crv), matrix(c(121,132,132,144), 2, 2) )
cprod(tcrossprod(m, m), tcrossprod(Cm, Cm), matrix(c(10,14,14,20), 2, 2) )
cprod(tcrossprod(m), tcrossprod(Cm), matrix(c(10,14,14,20), 2, 2) )
cprod(tcrossprod(m, rv), tcrossprod(Cm, Crv), matrix(c(47,70), 2, 1) )
cprod(tcrossprod(rv, m), tcrossprod(Crv, Cm), matrix(c(47,70), 1, 2) )
cprod(tcrossprod(v, m), tcrossprod(Cv, Cm), matrix(c(47,70), 1, 2) )
cprod(tcrossprod(rv, rv), tcrossprod(Crv, Crv), matrix(265,1,1) )
cprod(tcrossprod(cv, cv), tcrossprod(Ccv, Ccv), matrix(c(121,132,132,144), 2, 2) )
cprod(tcrossprod(v, v), tcrossprod(Cv, Cv), matrix(c(121,132,132,144), 2, 2) )
## non-square matrix, with Inf
m1 <- matrix(c(1,2,Inf,4,5,6), ncol=2)
m2 <- matrix(c(1,2,3,4), ncol=2)
v <- c(11,12)
rv <- v; dim(rv) <- c(1,2)
cv <- v; dim(cv) <- c(2,1)
v1 <- c(11,12,13)
rv1 <- v1; dim(rv1) <- c(1,3)
cv1 <- v1; dim(cv1) <- c(3,1)
Cm1 <- m1+0*1i
Cm2 <- m2+0*1i
Cv <- v+0*1i
Crv <- rv+0*1i
Ccv <- cv+0*1i
Cv1 <- v1+0*1i
Crv1 <- rv1+0*1i
Ccv1 <- cv1+0*1i
cprod(m1 %*% m2, Cm1 %*% Cm2, matrix(c(9,12,Inf,19,26,Inf), 3, 2) )
cprod(m1 %*% cv, Cm1 %*% Ccv, matrix(c(59,82,Inf), 3, 1) )
# the following 7 lines fail with Rblas and matprod = "blas"
cprod(rv1 %*% m1, Crv1 %*% Cm1, matrix(c(Inf,182), 1, 2) )
cprod(crossprod(m1, m1), crossprod(Cm1, Cm1), matrix(c(Inf,Inf,Inf,77), 2, 2) )
cprod(crossprod(m1, cv1), crossprod(Cm1, Ccv1), matrix(c(Inf,182), 2, 1) )
cprod(crossprod(cv1, m1), crossprod(Ccv1, Cm1), matrix(c(Inf,182), 1, 2) )
cprod(tcrossprod(m1, m1), tcrossprod(Cm1, Cm1), matrix(c(17,22,Inf,22,29,Inf,Inf,Inf,Inf), 3,3) )
cprod(tcrossprod(m2, m1), tcrossprod(Cm2, Cm1), matrix(c(13,18,17,24,Inf,Inf), 2, 3) )
cprod(tcrossprod(rv, m1), tcrossprod(Crv, Cm1), matrix(c(59,82,Inf), 1, 3) )
# the previous 7 lines fail with Rblas and matprod = "blas"
cprod(tcrossprod(m1, rv), tcrossprod(Cm1, Crv), matrix(c(59,82,Inf), 3, 1) )
## complex
m1 <- matrix(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), ncol=2)
m2 <- matrix(c(1+1i,2+2i,3+3i,4+4i), ncol=2)
stopifnot(identical(m1 %*% m2, matrix(c(18i,24i,30i,38i,52i,66i), 3, 2) ))
stopifnot(identical(crossprod(m1, m1), t(m1) %*% m1))
stopifnot(identical(tcrossprod(m1, m1), m1 %*% t(m1)))
}
## check that propagation of NaN/Inf values in multiplication of complex
## numbers is the same as in multiplication of complex matrices
for(mopt in c("default","internal","default.simd")) {
# matprod="blas" is excluded because some tests fail due to issues
# in NaN/Inf propagation even in Rblas
options(matprod=mopt)
vals <- c(0, 1, NaN, Inf)
for(ar in vals)
for(ai in vals)
for(br in vals)
for(bi in vals) {
a = ar + 1i * ai
b = br + 1i * bi
stopifnot(identical(a * b, as.complex(a %*% b)))
stopifnot(identical(a * b, as.complex(crossprod(a,b))))
stopifnot(identical(a * b, as.complex(tcrossprod(a,b))))
}
}
|
