     k -= kb;
   }
 
-  for (int ll = 0; ll < k; ll++) {
-    saxpy(ll + 1, B[0], a[ll], x);
-    B += ldb;
-  }
+  if (k >  0) { INNER_RIGHT_LOOP_INC( 1); B += ldb;
+  if (k >  1) { INNER_RIGHT_LOOP_INC( 2); B += ldb;
+  if (k >  2) { INNER_RIGHT_LOOP_INC( 3); B += ldb;
+  if (k >  3) { INNER_RIGHT_LOOP_INC( 4); B += ldb;
+  if (k >  4) { INNER_RIGHT_LOOP_INC( 5); B += ldb;
+  if (k >  5) { INNER_RIGHT_LOOP_INC( 6); B += ldb;
+  if (k >  6) { INNER_RIGHT_LOOP_INC( 7); B += ldb;
+  if (k >  7) { INNER_RIGHT_LOOP_INC( 8); B += ldb;
+  if (k >  8) { INNER_RIGHT_LOOP_INC( 9); B += ldb;
+  if (k >  9) { INNER_RIGHT_LOOP_INC(10); B += ldb;
+  if (k > 10) { INNER_RIGHT_LOOP_INC(11); B += ldb;
+  if (k > 11) { INNER_RIGHT_LOOP_INC(12); B += ldb;
+  if (k > 12) { INNER_RIGHT_LOOP_INC(13); B += ldb;
+  if (k > 13) { INNER_RIGHT_LOOP_INC(14); B += ldb;
+  if (k > 14) { INNER_RIGHT_LOOP_INC(15); B += ldb;
+  if (k > 15) { INNER_RIGHT_LOOP_INC(16); B += ldb; }}}}}}}}}}}}}}}}
 
   // Process non-diagonal blocks as for SGEMM
   k = n - bj - nb;


    if (k < kb) break;

    if (diag != CBlasNonUnit)
      x[0] += B[0];
    else
      saxpy( 1, B[0], a[ 0], x);
    B += ldb;
    INNER_RIGHT_LOOP_INC( 2); B += ldb;
    INNER_RIGHT_LOOP_INC( 3); B += ldb;
    INNER_RIGHT_LOOP_INC( 4); B += ldb;

    k -= kb;
  }

  if (k >  0) { INNER_RIGHT_LOOP_INC( 1); B += ldb;
  if (k >  1) { INNER_RIGHT_LOOP_INC( 2); B += ldb;
  if (k >  2) { INNER_RIGHT_LOOP_INC( 3); B += ldb;
  if (k >  3) { INNER_RIGHT_LOOP_INC( 4); B += ldb;
  if (k >  4) { INNER_RIGHT_LOOP_INC( 5); B += ldb;
  if (k >  5) { INNER_RIGHT_LOOP_INC( 6); B += ldb;
  if (k >  6) { INNER_RIGHT_LOOP_INC( 7); B += ldb;
  if (k >  7) { INNER_RIGHT_LOOP_INC( 8); B += ldb;
  if (k >  8) { INNER_RIGHT_LOOP_INC( 9); B += ldb;
  if (k >  9) { INNER_RIGHT_LOOP_INC(10); B += ldb;
  if (k > 10) { INNER_RIGHT_LOOP_INC(11); B += ldb;
  if (k > 11) { INNER_RIGHT_LOOP_INC(12); B += ldb;
  if (k > 12) { INNER_RIGHT_LOOP_INC(13); B += ldb;
  if (k > 13) { INNER_RIGHT_LOOP_INC(14); B += ldb;
  if (k > 14) { INNER_RIGHT_LOOP_INC(15); B += ldb;
  if (k > 15) { INNER_RIGHT_LOOP_INC(16); B += ldb; }}}}}}}}}}}}}}}}

  // Process non-diagonal blocks as for SGEMM
  k = n - bj - nb;
