docs: public static void computeZ(final int blockLength , final DSubmatrixD1 Y , final DSubmatrixD1 W,
                                  final int col , final double []temp , final double beta )
      {
          final int width = Y.col1-Y.col0;
  
          final double dataW[] = W.original.data;
          final double dataY[] = Y.original.data;
  
          final int colsW = W.original.numCols;
  
          final double beta_neg = -beta;
  
          for( int i = Y.row0; i < Y.row1; i += blockLength ) {
              int heightW = Math.min( blockLength , Y.row1 - i );
  
              int indexW = i*colsW + heightW*W.col0;
              int indexZ = i*colsW + heightW*W.col0 + col;
              int indexV = i*Y.original.numCols + heightW*Y.col0 + col;
  
              if( i == Y.row0 ) {
                  // handle the triangular portion with the leading zeros and the one
                  for( int k = 0; k < heightW; k++ , indexZ += width, indexW += width , indexV += width ) {
                      // compute the rows of W * h
                      double total = 0;
  
                      for( int j = 0; j < col; j++ ) {
                          total += dataW[indexW+j] * temp[j];
                      }
  
                      // add the two vectors together and multiply by -beta
                      if( k < col ) {  // zeros
                          dataW[indexZ] = -beta*total;
                      } else if( k == col ) { // one
                          dataW[indexZ] = beta_neg*(1.0 + total);
                      } else { // normal data
                          dataW[indexZ] = beta_neg*(dataY[indexV] + total);
                      }
                  }
              } else {
                  int endZ = indexZ + width*heightW;
  //                for( int k = 0; k < heightW; k++ ,
                  while( indexZ != endZ ) {
                      // compute the rows of W * h
                      double total = 0;
  
                      for( int j = 0; j < col; j++ ) {
                          total += dataW[indexW+j] * temp[j];
                      }
  
                      // add the two vectors together and multiply by -beta
                      dataW[indexZ] = beta_neg*(dataY[indexV] + total);
  
                      indexZ += width; indexW += width; indexV += width;
                  }
              }
          }
      }