I am interested in solving the linear system of equations Ax = b where A is a lower-triangular matrix (n × n) and b is a (n × 1) vector where n ≈ 600k.

I coded up backsubstitution in R and it works fast for matrices of size upto 1000, but is really slow for larger n (≈ 600k). i know that the naive backsubstitution is O (n ^ 2).

My R function is below; does anyone know of a more efficient (vectorial, parallel, etc.) way of doing it that scales to large n?

Backsubstitution

backsub=function(X,y)
{ 
 l=dim(X)
 n=l[1]  
 p=l[2]
 y=as.matrix(y)

   for (j in seq(p,1,-1))
    {  
      y[j,1]=y[j,1]/X[j,j]
      if((j-1)>0)
         y[1:(j-1),1]=y[1:(j-1),1]-(y[j,1]*X[1:(j-1),j])
    }  
    return(y)
}