38
39 # multiply one ut matrix by another, returning another ut instance;
40 # implement as a binary operation
41 "%mut%" <- function(utmat1,utmat2) {
42 n <- length(utmat1$ix) # numbers of rows and cols of matrix
43 utprod <- ut(matrix(0,nrow=n,ncol=n))
44 for (i in 1:n) { # compute col i of product
45 # let a[j] and bj denote columns j of utmat1 and utmat2, respectively,
46 # so that, e.g. b2[1] means element 1 of column 2 of utmat2
47 # then column i of product is equal to
48 # bi[1]*a[1] + ... + bi[i]*a[i]
49 # find index of start of column i in utmat2
50 startbi <- utmat2$ix[i]
51 # initialize vector that will become bi[1]*a[1] + ... + bi[i]*a[i]
52 prodcoli <- rep(0,i)
53 for (j in 1:i) { # find bi[j]*a[j], add to prodcoli
54 startaj <- utmat1$ix[j]
55 bielement <- utmat2$mat[startbi+j-1]
56 prodcoli[1:j] <- prodcoli[1:j] +
57 bielement*utmat1$mat[startaj:(startaj+j-1)]
58 }
59 # now need to tack on the lower 0s
60 startprodcoli <- sum1toi(i-1)+1
61 utprod$mat[startbi:(startbi+i-1)] <- prodcoli
62 }
63 return(utprod)
64 }Let’s test it.> test
function() {
utm1 <- ut(rbind(1:2,c(0,2)))
utm2 <- ut(rbind(3:2,c(0,1)))
utp <- utm1 %mut% utm2
print(utm1)
print(utm2)
print(utp)
utm1 <- ut(rbind(1:3,0:2,c(0,0,5)))
utm2 <- ut(rbind(4:2,0:2,c(0,0,1)))
utp <- utm1 %mut% utm2
print(utm1)
print(utm2)
print(utp)
}216 Chapter 9