Data Structures, Algorithms, & Applications in Java
Chapter 8, Exercise 31

The general forumla to compute the product c of two n x n matrices a and b is
c(i,j) = sum (k=1 to n) (a(i,k) * b(k,j))

Since a and b are lower triangular matrices, a(i,j) = 0 and b(i,j) = 0 for i < j. Using this knowledge in the formula for c(i,j), we see that when i < j, if k > i, a(i,k) = 0; and if k <= i, then k < j and so b(k,j) = 0. Therefore, c(i,j) = 0. Hence, the product of two lower-triangular matrices is also a lower-triangular matrix.

If we eliminate terms in the matrix product formula that involve an element in the upper triangle of a or b, the formula becomes
c(i,j) = sum (k=i to j) (a(i,k) * b(k,j))

The preceding development results in the code that is given below. The complexity of is O(rows3).
public class LowerTriangularMatrixWithMultiply extends LowerTriangularMatrix
{
   public LowerTriangularMatrixWithMultiply(int theRows, Object theZero)
      {super(theRows, theZero);}
   
   /** return the this * b */
   public LowerTriangularMatrixWithMultiply multiply
                         (LowerTriangularMatrixWithMultiply b)
   {
     if (rows != b.rows)
        throw new IllegalArgumentException
                  ("Matrices must have same dimensions");
  
     // create result array w
     LowerTriangularMatrixWithMultiply w =
                    new LowerTriangularMatrixWithMultiply(rows, zero);

     for (int i = 1; i <= rows; i++)
        for (int j = 1; j <= i; j++)
        {// compute w(i,j)
           // compute this(i,j) * b(j,j)
           Computable sum = (Computable)
                            ((Computable) element[i * (i - 1) / 2 + j - 1]).
                            multiply(b.element[j * (j - 1) / 2 + j - 1]);
           // add in remaining terms
           for (int k = j + 1; k <= i; k++)
              sum.increment(((Computable) element[i * (i - 1) / 2 + k - 1]).
                            multiply(b.element[k * (k - 1) / 2 + j - 1]));

           // save as w(i,j)
           w.element[i * (i - 1) / 2 + j - 1] = sum;
      }
  
     return w;
   }
}



We can improve the run time slightly by avoiding the recomputation of i*(i-1)/2 - 1 and j - 1 inside the innermost for loop.