>  # mxmult.r         common pitfall in matrix multiplication
>  #                  September 2008, revised September 2009
>  #
>  # slr design matrix X
>   X <- cbind( rep(1,6), (1:6) )
>   X
     [,1] [,2]
[1,]    1    1
[2,]    1    2
[3,]    1    3
[4,]    1    4
[5,]    1    5
[6,]    1    6
>  # weight vector w
>   w <- sqrt( (1:6) )          # just so they're different
>   w
[1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490
>  # weight matrix
>   W <- diag(w)
>  # now (X'inv(W) X) four ways:
>  #   first          correct, but slow, W takes lots of space
>   W <- diag(w)
>   xwx1 <- t(X) %*% solve(W) %*% X      
>   xwx1
          [,1]     [,2]
[1,]  3.639919 10.83182
[2,] 10.831822 42.90186
>  #   second         correct, less slow, takes lots of space 
>   xwx2 <- t(X) %*% diag(1/w) %*% X    
>   xwx2
          [,1]     [,2]
[1,]  3.639919 10.83182
[2,] 10.831822 42.90186
>  #   third          wrong, but faster -- recycles w wrong
>   xwx3 <- ( t(X) / w ) %*% X           # *****wrong*****
>   xwx3
          [,1]     [,2]
[1,]  4.049128 13.06637
[2,] 10.709769 44.89199
>  #   fourth         correct, but looks wrong
>   xwx4 <- t(X) %*% (X/w)               # uses recycling correctly, fast
>   xwx4
          [,1]     [,2]
[1,]  3.639919 10.83182
[2,] 10.831822 42.90186
>  #                  looks different but same execution
>   crossprod(X,X/w)                     # correct, fast
          [,1]     [,2]
[1,]  3.639919 10.83182
[2,] 10.831822 42.90186
>  # done
>   rm( list=ls() )
>   q()