Figure 1: Decomposition of v into a diagonal and triangular matrix

for(j = 2 to ndf) // assumes 1-based indexing
{
   nj1 = startrow(j);
  
   for (i = nj1 + 1 to j - 1)
   {
      ni1 = max(nj1, startrow(i))
      v(i, j) = 
         v(i, j) - dotproduct(v(k, i), v(k, j), k = ni1 .. i-1)
   }

   temp(k)  = v(k, j) / diag(k)  // k loops over the column
   diag(j)  = diag(j) - dotproduct(v(k, j), temp(k))
  
   v(k, j)  = temp(k)  // k loops over the column
}

- End of Figure -