# Linear algebra in R

#   Create four column vectors with three components each

c1 = c(2, 1, 0)      # column 1
c2 = c(1, 2, 1)      # column 2
c3 = c(0, 1, 2)      # column 3
c4 = c(0, 0, 1)      # column 4

# use matrix and cbind to create a 3 X 3 matrix with the first 3 columns

A = matrix( cbind( c1, c2, c3), nrow=3, ncol=3)

cat("\n Print some entries of A\n\n")

cat("A[1,1] is ", A[1,1], "\n")
cat("A[1,2] is ", A[1,2], "\n")
cat("A[3,1] is ", A[3,1], "\n")

cat("\n Print A, row by row\n\n")

cat("Row 1 of A is ", A[1,1:3], "\n")
cat("Row 2 of A is ", A[2,1:3], "\n")
cat("Row 3 of A is ", A[3,1:3], "\n")

B = A %*% A    # Matrix multiplication is %*%, not just *

cat("\n Print B = A*A, row by row\n\n")

cat("Row 1 of B is ", B[1,1:3], "\n")
cat("Row 2 of B is ", B[2,1:3], "\n")
cat("Row 3 of B is ", B[3,1:3], "\n")

B = A %*% A + 2*A   # But matrix addition is + and scalar multiplication is *

cat("\n Print B = A*A + 2*A, row by row\n\n")

cat("Row 1 of B is ", B[1,1:3], "\n")
cat("Row 2 of B is ", B[2,1:3], "\n")
cat("Row 3 of B is ", B[3,1:3], "\n")

x = c(-1, 3, 2)   #  create a column vector in 3D

cat("\n Print the entries of x\n\n")

cat("x[1] is ", x[1], "\n")
cat("x[2] is ", x[2], "\n")
cat("x[3] is ", x[3], "\n")

cat("\n Printing x as a whole looks like a row vector: ", x, "\n")

y = A %*% x      # Matrix vector multiplication is also %*%

cat("Print all the entries of y = A*x: ", y , "\n")

C = matrix( cbind( c1, c2, c3, c4), nrow=3, ncol=4)  # Now include the 4-th column, a non-square matrix

cat("\n Print the rows of C, which has 3 rows and 4 columns \n\n")

cat("Row 1 of C is ", C[1,1:4], "\n")
cat("Row 2 of C is ", C[2,1:4], "\n")
cat("Row 3 of C is ", C[3,1:4], "\n")

cat("\n C*C is not compatible for multiplication, uncomment and see the error message \n\n")

# D = C %*% C    # uncomment this to see an error message

Ct = t(C)  # Ct is transpose(C)

cat("\n Print the rows of E = C * (transpose(C)), which has 3 rows and 3 columns \n\n")

E = C %*% Ct  # Ct is transpose(C)

cat("Row 1 of E is ", E[1,1:3], "\n")
cat("Row 2 of E is ", E[2,1:3], "\n")
cat("Row 3 of E is ", E[3,1:3], "\n")

cat("\n Print the rows of Ai = inverse(A) \n\n")

Ai = solve(A)     # Ai is inverse of A, computed by the solve() function in R

cat("Row 1 of Ai is ", Ai[1,1:3], "\n")
cat("Row 2 of Ai is ", Ai[2,1:3], "\n")
cat("Row 3 of Ai is ", Ai[3,1:3], "\n")

Id = Ai %*% A     # A * inverse(A) should be the identity matrix

cat("\n Print the rows of Id = identity matrix, if the inverse was right \n\n")

cat("Row 1 of Id is ", Id[1,1:3], "\n")
cat("Row 2 of Id is ", Id[2,1:3], "\n")
cat("Row 3 of Id is ", Id[3,1:3], "\n")