LDU decomposition and pivots

Making a pivot

pivot<-function(x,i,j) {
p<-x[i,j]; r<-x[i,]/p; c<-x[,j]/p
return(x-p*outer(c,r))
}

set.seed(12345)
print (a <- matrix (rnorm (9), 3, 3))

##            [,1]       [,2]       [,3]
## [1,]  0.5855288 -0.4534972  0.6300986
## [2,]  0.7094660  0.6058875 -0.2761841
## [3,] -0.1093033 -1.8179560 -0.2841597

pivot(a, 2, 3)

##            [,1]       [,2]          [,3]
## [1,]  2.2041358  0.9288011 -1.110223e-16
## [2,]  0.0000000  0.0000000  0.000000e+00
## [3,] -0.8392573 -2.4413402  0.000000e+00

print (a1<-pivot (a, 1, 1))

##      [,1]      [,2]       [,3]
## [1,]    0  0.000000  0.0000000
## [2,]    0  1.155375 -1.0396538
## [3,]    0 -1.902612 -0.1665364

print (a2 <- pivot (a1, 2, 3))

##      [,1]      [,2]         [,3]
## [1,]    0  0.000000 0.000000e+00
## [2,]    0  0.000000 0.000000e+00
## [3,]    0 -2.087685 2.775558e-17

print (a3 <- pivot (a2, 3, 2))

##      [,1] [,2] [,3]
## [1,]    0    0    0
## [2,]    0    0    0
## [3,]    0    0    0

Pivot along diagonal

pivot_along_diagonal<-function(x) {
n<-nrow(x); d<-diag(n);l<-r<-matrix(0,n,n)
for (i in 1:n) {
print(x)
p<-x[i,i]; u<-x[i,]/p; v<-x[,i]/p
x<-x-p*outer(v,u)
l[,i]<-v; r[i,]<-u;d[i,i]<-p
}
return(list(l=l,r=r,d=d))
}

print (lst <- pivot_along_diagonal (a))

##            [,1]       [,2]       [,3]
## [1,]  0.5855288 -0.4534972  0.6300986
## [2,]  0.7094660  0.6058875 -0.2761841
## [3,] -0.1093033 -1.8179560 -0.2841597
##      [,1]      [,2]       [,3]
## [1,]    0  0.000000  0.0000000
## [2,]    0  1.155375 -1.0396538
## [3,]    0 -1.902612 -0.1665364
##      [,1] [,2]      [,3]
## [1,]    0    0  0.000000
## [2,]    0    0  0.000000
## [3,]    0    0 -1.878585
## $l
##            [,1]      [,2] [,3]
## [1,]  1.0000000  0.000000    0
## [2,]  1.2116671  1.000000    0
## [3,] -0.1866745 -1.646749    1
## 
## $r
##      [,1]       [,2]      [,3]
## [1,]    1 -0.7745087  1.076119
## [2,]    0  1.0000000 -0.899841
## [3,]    0  0.0000000  1.000000
## 
## $d
##           [,1]     [,2]      [,3]
## [1,] 0.5855288 0.000000  0.000000
## [2,] 0.0000000 1.155375  0.000000
## [3,] 0.0000000 0.000000 -1.878585

What shall I pivot ?

pivot<-function(x) {
n<-nrow(x); m<-ncol(x); ii<-rep(0,n); jj<-
rep(0,m)
l<-u<-d<-matrix(0,n,m); k<-1
repeat{
print(x)
if (substr(readline("Continue ? "),1,1) == "n")
break()
i<-as.integer(readline("Row Index ? "))
j<-as.integer(readline("Column Index ? "))
ii[k]<-i; jj[k]<-j
p<-x[i,j]; r<-x[i,]/p; c<-x[,j]/p
x<-x-p*outer(c,r)
l[,i]<-c; u[j,]<-r; d[i,i]<-p
k<-k+1
}
return(list(l=l,u=u,d=d))
}