code: public static DoubleMatrix2D subdiagonalMultiply(final DoubleMatrix2D A, final DoubleMatrix2D B){
  		final int r = A.rows();
  		final int rc = A.columns();
  		final int c = B.columns();
  
  		if(r != c){
  			throw new IllegalArgumentException("The result must be square");
  		}
  
  		boolean useSparsity = A instanceof SparseDoubleMatrix2D;
  		DoubleFactory2D F2 = (useSparsity)? DoubleFactory2D.sparse : DoubleFactory2D.dense; 
  		final DoubleMatrix2D ret = F2.make(r, c);
  
  		if(useSparsity){
  			IntIntDoubleFunction myFunct = new IntIntDoubleFunction() {
  				@Override
  				public double apply(int t, int s, double pts) {
  					int i = t;
  					for (int j = 0; j < i + 1; j++) {
  						ret.setQuick(i, j, ret.getQuick(i, j) + pts * B.getQuick(s, j));
  					}
  					return pts;
  				}
  			};
  
  			//view A row by row
  			A.forEachNonZero(myFunct);
  		}else{
  			for (int i = 0; i < r; i++) {
  				for (int j = 0; j < i + 1; j++) {
  					double s = 0;
  					for (int k = 0; k < rc; k++) {
  						s += A.getQuick(i, k) * B.getQuick(k, j);
  					}
  					ret.setQuick(i, j, s);
  				}
  			}
  		}
  
  		return ret;
  	}
docs: Return the sub-diagonal result of the multiplication.
  If A is sparse, returns a sparse matrix (even if, generally speaking,
  the multiplication of two sparse matrices is not sparse) because the result
  is at least 50% (aside the diagonal elements) sparse.