docs: public static void solveR( final int blockLength ,
                                 final DSubmatrixD1 R,
                                 final DSubmatrixD1 B ,
                                 boolean transR ) {
  
          int lengthR = B.row1 - B.row0;
          if( R.getCols() != lengthR ) {
              throw new IllegalArgumentException("Number of columns in R must be equal to the number of rows in B");
          } else if( R.getRows() != lengthR ) {
              throw new IllegalArgumentException("Number of rows in R must be equal to the number of rows in B");
          }
  
          DSubmatrixD1 Y = new DSubmatrixD1(B.original);
  
          DSubmatrixD1 Rinner = new DSubmatrixD1(R.original);
          DSubmatrixD1 Binner = new DSubmatrixD1(B.original);
  
          int startI,stepI;
  
          if( transR ) {
              startI = 0;
              stepI = blockLength;
          } else {
              startI = lengthR - lengthR % blockLength;
              if( startI == lengthR && lengthR >= blockLength )
                  startI -= blockLength;
  
              stepI = -blockLength;
          }
  
          for( int i = startI; ; i += stepI ) {
              if( transR ) {
                  if( i >= lengthR ) break;
              } else {
                  if( i < 0 ) break;
              }
  
              // width and height of the inner T(i,i) block
              int widthT = Math.min(blockLength, lengthR-i);
  
              Rinner.col0 = R.col0 + i;    Rinner.col1 = Rinner.col0 + widthT;
              Rinner.row0 = R.row0 + i;    Rinner.row1 = Rinner.row0 + widthT;
  
              Binner.col0 = B.col0;       Binner.col1 = B.col1;
              Binner.row0 = B.row0 + i;   Binner.row1 = Binner.row0 + widthT;
  
              // solve the top row block
              // B(i,:) = T(i,i)^-1 Y(i,:)
              solveBlock(blockLength,true, Rinner,Binner,transR,false);
  
              boolean updateY;
              if( transR ) {
                  updateY = Rinner.row1 < R.row1;
              } else {
                  updateY = Rinner.row0 > 0;
              }
              if( updateY ) {
                  // Y[i,:] = Y[i,:] - sum j=1:i-1 { T[i,j] B[j,i] }
                  // where i is the next block down
                  // The summation is a block inner product
                  if( transR ) {
                      Rinner.col0 = Rinner.col1;
                      Rinner.col1 = Math.min(Rinner.col0+blockLength, R.col1);
                      Rinner.row0 = R.row0;
                      //Rinner.row1 = Rinner.row1;
  
                      Binner.row0 = B.row0;
                      //Binner.row1 = Binner.row1;
  
                      Y.row0 = Binner.row1;
                      Y.row1 = Math.min(Y.row0+blockLength,B.row1);
                  } else {
                      Rinner.row1 = Rinner.row0;
                      Rinner.row0 = Rinner.row1 - blockLength;
                      Rinner.col1 = R.col1;
  
  //                    Binner.row0 = Binner.row0;
                      Binner.row1 = B.row1;
  
                      Y.row0 = Binner.row0-blockLength;
                      Y.row1 = Binner.row0;
                  }
  
                  // step through each block column
                  for( int k = B.col0; k < B.col1; k += blockLength ) {
  
                      Binner.col0 = k;
                      Binner.col1 = Math.min(k+blockLength,B.col1);
  
                      Y.col0 = Binner.col0;
                      Y.col1 = Binner.col1;
  
                      if( transR ) {
                          // Y = Y - T^T * B
                          MatrixMult_DDRB.multMinusTransA(blockLength, Rinner,Binner,Y);
                      } else {
                          // Y = Y - T * B
                          MatrixMult_DDRB.multMinus(blockLength, Rinner,Binner,Y);
                      }
                  }
              }
          }
      }