code: def tsqr(data, compute_svd=False, _max_vchunk_size=None):
      """Direct Tall-and-Skinny QR algorithm
  
      As presented in:
  
          A. Benson, D. Gleich, and J. Demmel.
          Direct QR factorizations for tall-and-skinny matrices in
          MapReduce architectures.
          IEEE International Conference on Big Data, 2013.
          https://arxiv.org/abs/1301.1071
  
      This algorithm is used to compute both the QR decomposition and the
      Singular Value Decomposition.  It requires that the input array have a
      single column of blocks, each of which fit in memory.
  
      Parameters
      ----------
      data: Array
      compute_svd: bool
          Whether to compute the SVD rather than the QR decomposition
      _max_vchunk_size: Integer
          Used internally in recursion to set the maximum row dimension
          of chunks in subsequent recursive calls.
  
      Notes
      -----
      With ``k`` blocks of size ``(m, n)``, this algorithm has memory use that
      scales as ``k * n * n``.
  
      The implementation here is the recursive variant due to the ultimate
      need for one "single core" QR decomposition. In the non-recursive version
      of the algorithm, given ``k`` blocks, after ``k`` ``m * n`` QR
      decompositions, there will be a "single core" QR decomposition that will
      have to work with a ``(k * n, n)`` matrix.
  
      Here, recursion is applied as necessary to ensure that ``k * n`` is not
      larger than ``m`` (if ``m / n >= 2``). In particular, this is done
      to ensure that single core computations do not have to work on blocks
      larger than ``(m, n)``.
  
      Where blocks are irregular, the above logic is applied with the "height" of
      the "tallest" block used in place of ``m``.
  
      Consider use of the ``rechunk`` method to control this behavior.
      Taller blocks will reduce overall memory use (assuming that many of them
      still fit in memory at once).
  
      See Also
      --------
      dask.array.linalg.qr
          Powered by this algorithm
      dask.array.linalg.svd
          Powered by this algorithm
      dask.array.linalg.sfqr
          Variant for short-and-fat arrays
      """
      nr, nc = len(data.chunks[0]), len(data.chunks[1])
      cr_max, cc = max(data.chunks[0]), data.chunks[1][0]
  
      if not (
          data.ndim == 2 and nc == 1
      ):  # Is a matrix  # Only one column block
          raise ValueError(
              "Input must have the following properties:\n"
              "  1. Have two dimensions\n"
              "  2. Have only one column of blocks\n\n"
              "Note: This function (tsqr) supports QR decomposition in the case of\n"
              "tall-and-skinny matrices (single column chunk/block; see qr)\n"
              "Current shape: {},\nCurrent chunksize: {}".format(
                  data.shape, data.chunksize
              )
          )
  
      token = "-" + tokenize(data, compute_svd)
  
      m, n = data.shape
      numblocks = (nr, 1)
  
      qq, rr = np.linalg.qr(np.ones(shape=(1, 1), dtype=data.dtype))
  
      layers = data.__dask_graph__().layers.copy()
      dependencies = data.__dask_graph__().dependencies.copy()
  
      # Block qr
      name_qr_st1 = "qr" + token
      dsk_qr_st1 = blockwise(
          _wrapped_qr,
          name_qr_st1,
          "ij",
          data.name,
          "ij",
          numblocks={data.name: numblocks},
      )
      layers[name_qr_st1] = dsk_qr_st1
      dependencies[name_qr_st1] = set(data.__dask_layers__())
  
      # Block qr[0]
      name_q_st1 = "getitem" + token + "-q1"
      dsk_q_st1 = {
          (name_q_st1, i, 0): (operator.getitem, (name_qr_st1, i, 0), 0)
          for i in range(numblocks[0])
      }
      layers[name_q_st1] = dsk_q_st1
      dependencies[name_q_st1] = {name_qr_st1}
  
      # Block qr[1]
      name_r_st1 = "getitem" + token + "-r1"
      dsk_r_st1 = {
          (name_r_st1, i, 0): (operator.getitem, (name_qr_st1, i, 0), 1)
          for i in range(numblocks[0])
      }
      layers[name_r_st1] = dsk_r_st1
      dependencies[name_r_st1] = {name_qr_st1}
  
      # Next step is to obtain a QR decomposition for the stacked R factors, so either:
      # - gather R factors into a single core and do a QR decomposition
      # - recurse with tsqr (if single core computation too large and a-priori "meaningful
      #   reduction" possible, meaning that chunks have to be well defined)
  
      single_core_compute_m = nr * cc
      chunks_well_defined = not any(
          np.isnan(c) for cs in data.chunks for c in cs
      )
      prospective_blocks = np.ceil(single_core_compute_m / cr_max)
      meaningful_reduction_possible = (
          cr_max if _max_vchunk_size is None else _max_vchunk_size
      ) >= 2 * cc
      can_distribute = (
          chunks_well_defined and int(prospective_blocks) > 1
      )
  
      if (
          chunks_well_defined
          and meaningful_reduction_possible
          and can_distribute
      ):
          # stack chunks into blocks and recurse using tsqr
  
          # Prepare to stack chunks into blocks (from block qr[1])
          all_blocks = []
          curr_block = []
          curr_block_sz = 0
          for idx, a_m in enumerate(data.chunks[0]):
              m_q = a_m
              n_q = min(m_q, cc)
              m_r = n_q
              # n_r = cc
              if curr_block_sz + m_r > cr_max:
                  all_blocks.append(curr_block)
                  curr_block = []
                  curr_block_sz = 0
              curr_block.append((idx, m_r))
              curr_block_sz += m_r
          if len(curr_block) > 0:
              all_blocks.append(curr_block)
  
          # R_stacked
          name_r_stacked = "stack" + token + "-r1"
          dsk_r_stacked = {
              (name_r_stacked, i, 0): (
                  np.vstack,
                  (
                      tuple,
                      [
                          (name_r_st1, idx, 0)
                          for idx, _ in sub_block_info
                      ],
                  ),
              )
              for i, sub_block_info in enumerate(all_blocks)
          }
          layers[name_r_stacked] = dsk_r_stacked
          dependencies[name_r_stacked] = {name_r_st1}
  
          # retrieve R_stacked for recursion with tsqr
          vchunks_rstacked = tuple(
              sum(map(lambda x: x[1], sub_block_info))
              for sub_block_info in all_blocks
          )
          graph = HighLevelGraph(layers, dependencies)
          # dsk.dependencies[name_r_stacked] = {data.name}
          r_stacked_meta = meta_from_array(
              data, len((sum(vchunks_rstacked), n)), dtype=rr.dtype
          )
          r_stacked = Array(
              graph,
              name_r_stacked,
              shape=(sum(vchunks_rstacked), n),
              chunks=(vchunks_rstacked, n),
              meta=r_stacked_meta,
          )
  
          # recurse
          q_inner, r_inner = tsqr(r_stacked, _max_vchunk_size=cr_max)
          layers = toolz.merge(q_inner.dask.layers, r_inner.dask.layers)
          dependencies = toolz.merge(
              q_inner.dask.dependencies, r_inner.dask.dependencies
          )
  
          # Q_inner: "unstack"
          name_q_st2 = "getitem" + token + "-q2"
          dsk_q_st2 = {
              (name_q_st2, j, 0): (
                  operator.getitem,
                  (q_inner.name, i, 0),
                  ((slice(e[0], e[1])), (slice(0, n))),
              )
              for i, sub_block_info in enumerate(all_blocks)
              for j, e in zip(
                  [x[0] for x in sub_block_info],
                  _cumsum_blocks([x[1] for x in sub_block_info]),
              )
          }
          layers[name_q_st2] = dsk_q_st2
          dependencies[name_q_st2] = set(q_inner.__dask_layers__())
  
          # R: R_inner
          name_r_st2 = "r-inner" + token
          dsk_r_st2 = {(name_r_st2, 0, 0): (r_inner.name, 0, 0)}
          layers[name_r_st2] = dsk_r_st2
          dependencies[name_r_st2] = set(r_inner.__dask_layers__())
  
          # Q: Block qr[0] (*) Q_inner
          name_q_st3 = "dot" + token + "-q3"
          dsk_q_st3 = blockwise(
              np.dot,
              name_q_st3,
              "ij",
              name_q_st1,
              "ij",
              name_q_st2,
              "ij",
              numblocks={name_q_st1: numblocks, name_q_st2: numblocks},
          )
          layers[name_q_st3] = dsk_q_st3
          dependencies[name_q_st3] = {name_q_st1, name_q_st2}
      else:
          # Do single core computation
  
          # Stacking for in-core QR computation
          to_stack = [(name_r_st1, i, 0) for i in range(numblocks[0])]
          name_r_st1_stacked = "stack" + token + "-r1"
          dsk_r_st1_stacked = {
              (name_r_st1_stacked, 0, 0): (np.vstack, (tuple, to_stack))
          }
          layers[name_r_st1_stacked] = dsk_r_st1_stacked
          dependencies[name_r_st1_stacked] = {name_r_st1}
  
          # In-core QR computation
          name_qr_st2 = "qr" + token + "-qr2"
          dsk_qr_st2 = blockwise(
              np.linalg.qr,
              name_qr_st2,
              "ij",
              name_r_st1_stacked,
              "ij",
              numblocks={name_r_st1_stacked: (1, 1)},
          )
          layers[name_qr_st2] = dsk_qr_st2
          dependencies[name_qr_st2] = {name_r_st1_stacked}
  
          # In-core qr[0]
          name_q_st2_aux = "getitem" + token + "-q2-aux"
          dsk_q_st2_aux = {
              (name_q_st2_aux, 0, 0): (
                  operator.getitem,
                  (name_qr_st2, 0, 0),
                  0,
              )
          }
          layers[name_q_st2_aux] = dsk_q_st2_aux
          dependencies[name_q_st2_aux] = {name_qr_st2}
  
          chucks_are_all_known = not any(
              np.isnan(c) for cs in data.chunks for c in cs
          )
          if chucks_are_all_known:
              # when chunks are all known...
              # obtain slices on q from in-core compute (e.g.: (slice(10, 20), slice(0, 5)))
              q2_block_sizes = [min(e, n) for e in data.chunks[0]]
              block_slices = [
                  (slice(e[0], e[1]), slice(0, n))
                  for e in _cumsum_blocks(q2_block_sizes)
              ]
              dsk_q_blockslices = {}
              deps = set()
          else:
              # when chunks are not already known...
  
              # request shape information: vertical chunk sizes & column dimension (n)
              name_q2bs = "shape" + token + "-q2"
              dsk_q2_shapes = {
                  (name_q2bs, i): (
                      min,
                      (getattr, (data.name, i, 0), "shape"),
                  )
                  for i in range(numblocks[0])
              }
              name_n = "getitem" + token + "-n"
              dsk_n = {
                  name_n: (
                      operator.getitem,
                      (getattr, (data.name, 0, 0), "shape"),
                      1,
                  )
              }
  
              # cumulative sums (start, end)
              name_q2cs = "cumsum" + token + "-q2"
              dsk_q2_cumsum = {(name_q2cs, 0): [0, (name_q2bs, 0)]}
  
              for i in range(1, numblocks[0]):
                  dsk_q2_cumsum[(name_q2cs, i)] = (
                      _cumsum_part,
                      (name_q2cs, i - 1),
                      (name_q2bs, i),
                  )
  
              # obtain slices on q from in-core compute (e.g.: (slice(10, 20), slice(0, 5)))
              name_blockslice = "slice" + token + "-q"
              dsk_block_slices = {
                  (name_blockslice, i): (
                      tuple,
                      [
                          (apply, slice, (name_q2cs, i)),
                          (slice, 0, name_n),
                      ],
                  )
                  for i in range(numblocks[0])
              }
  
              dsk_q_blockslices = toolz.merge(
                  dsk_n, dsk_q2_shapes, dsk_q2_cumsum, dsk_block_slices
              )
  
              deps = {data.name}
              block_slices = [
                  (name_blockslice, i) for i in range(numblocks[0])
              ]
  
          layers["q-blocksizes" + token] = dsk_q_blockslices
          dependencies["q-blocksizes" + token] = deps
  
          # In-core qr[0] unstacking
          name_q_st2 = "getitem" + token + "-q2"
          dsk_q_st2 = {
              (name_q_st2, i, 0): (
                  operator.getitem,
                  (name_q_st2_aux, 0, 0),
                  b,
              )
              for i, b in enumerate(block_slices)
          }
          layers[name_q_st2] = dsk_q_st2
          if chucks_are_all_known:
              dependencies[name_q_st2] = {name_q_st2_aux}
          else:
              dependencies[name_q_st2] = {
                  name_q_st2_aux,
                  "q-blocksizes" + token,
              }
  
          # Q: Block qr[0] (*) In-core qr[0]
          name_q_st3 = "dot" + token + "-q3"
          dsk_q_st3 = blockwise(
              np.dot,
              name_q_st3,
              "ij",
              name_q_st1,
              "ij",
              name_q_st2,
              "ij",
              numblocks={name_q_st1: numblocks, name_q_st2: numblocks},
          )
          layers[name_q_st3] = dsk_q_st3
          dependencies[name_q_st3] = {name_q_st1, name_q_st2}
  
          # R: In-core qr[1]
          name_r_st2 = "getitem" + token + "-r2"
          dsk_r_st2 = {
              (name_r_st2, 0, 0): (
                  operator.getitem,
                  (name_qr_st2, 0, 0),
                  1,
              )
          }
          layers[name_r_st2] = dsk_r_st2
          dependencies[name_r_st2] = {name_qr_st2}
  
      if not compute_svd:
          is_unknown_m = np.isnan(data.shape[0]) or any(
              np.isnan(c) for c in data.chunks[0]
          )
          is_unknown_n = np.isnan(data.shape[1]) or any(
              np.isnan(c) for c in data.chunks[1]
          )
  
          if is_unknown_m and is_unknown_n:
              # assumption: m >= n
              q_shape = data.shape
              q_chunks = (data.chunks[0], (np.nan,))
              r_shape = (np.nan, np.nan)
              r_chunks = ((np.nan,), (np.nan,))
          elif is_unknown_m and not is_unknown_n:
              # assumption: m >= n
              q_shape = data.shape
              q_chunks = (data.chunks[0], (n,))
              r_shape = (n, n)
              r_chunks = (n, n)
          elif not is_unknown_m and is_unknown_n:
              # assumption: m >= n
              q_shape = data.shape
              q_chunks = (data.chunks[0], (np.nan,))
              r_shape = (np.nan, np.nan)
              r_chunks = ((np.nan,), (np.nan,))
          else:
              q_shape = (
                  data.shape
                  if data.shape[0] >= data.shape[1]
                  else (data.shape[0], data.shape[0])
              )
              q_chunks = (
                  data.chunks
                  if data.shape[0] >= data.shape[1]
                  else (data.chunks[0], data.chunks[0])
              )
              r_shape = (
                  (n, n)
                  if data.shape[0] >= data.shape[1]
                  else data.shape
              )
              r_chunks = r_shape
  
          # dsk.dependencies[name_q_st3] = {data.name}
          # dsk.dependencies[name_r_st2] = {data.name}
          graph = HighLevelGraph(layers, dependencies)
          q_meta = meta_from_array(data, len(q_shape), qq.dtype)
          r_meta = meta_from_array(data, len(r_shape), rr.dtype)
          q = Array(
              graph,
              name_q_st3,
              shape=q_shape,
              chunks=q_chunks,
              meta=q_meta,
          )
          r = Array(
              graph,
              name_r_st2,
              shape=r_shape,
              chunks=r_chunks,
              meta=r_meta,
          )
          return q, r
      else:
          # In-core SVD computation
          name_svd_st2 = "svd" + token + "-2"
          dsk_svd_st2 = blockwise(
              np.linalg.svd,
              name_svd_st2,
              "ij",
              name_r_st2,
              "ij",
              numblocks={name_r_st2: (1, 1)},
          )
          # svd[0]
          name_u_st2 = "getitem" + token + "-u2"
          dsk_u_st2 = {
              (name_u_st2, 0, 0): (
                  operator.getitem,
                  (name_svd_st2, 0, 0),
                  0,
              )
          }
          # svd[1]
          name_s_st2 = "getitem" + token + "-s2"
          dsk_s_st2 = {
              (name_s_st2, 0): (
                  operator.getitem,
                  (name_svd_st2, 0, 0),
                  1,
              )
          }
          # svd[2]
          name_v_st2 = "getitem" + token + "-v2"
          dsk_v_st2 = {
              (name_v_st2, 0, 0): (
                  operator.getitem,
                  (name_svd_st2, 0, 0),
                  2,
              )
          }
          # Q * U
          name_u_st4 = "getitem" + token + "-u4"
          dsk_u_st4 = blockwise(
              dotmany,
              name_u_st4,
              "ij",
              name_q_st3,
              "ik",
              name_u_st2,
              "kj",
              numblocks={name_q_st3: numblocks, name_u_st2: (1, 1)},
          )
  
          layers[name_svd_st2] = dsk_svd_st2
          dependencies[name_svd_st2] = {name_r_st2}
          layers[name_u_st2] = dsk_u_st2
          dependencies[name_u_st2] = {name_svd_st2}
          layers[name_u_st4] = dsk_u_st4
          dependencies[name_u_st4] = {name_q_st3, name_u_st2}
          layers[name_s_st2] = dsk_s_st2
          dependencies[name_s_st2] = {name_svd_st2}
          layers[name_v_st2] = dsk_v_st2
          dependencies[name_v_st2] = {name_svd_st2}
  
          uu, ss, vvh = np.linalg.svd(
              np.ones(shape=(1, 1), dtype=data.dtype)
          )
  
          k = _nanmin(
              m, n
          )  # avoid RuntimeWarning with np.nanmin([m, n])
  
          m_u = m
          n_u = int(k) if not np.isnan(k) else k
          n_s = n_u
          m_vh = n_u
          n_vh = n
          d_vh = max(
              m_vh, n_vh
          )  # full matrix returned: but basically n
          graph = HighLevelGraph(layers, dependencies)
          u_meta = meta_from_array(data, len((m_u, n_u)), uu.dtype)
          s_meta = meta_from_array(data, len((n_s,)), ss.dtype)
          vh_meta = meta_from_array(data, len((d_vh, d_vh)), vvh.dtype)
          u = Array(
              graph,
              name_u_st4,
              shape=(m_u, n_u),
              chunks=(data.chunks[0], (n_u,)),
              meta=u_meta,
          )
          s = Array(
              graph,
              name_s_st2,
              shape=(n_s,),
              chunks=((n_s,),),
              meta=s_meta,
          )
          vh = Array(
              graph,
              name_v_st2,
              shape=(d_vh, d_vh),
              chunks=((n,), (n,)),
              meta=vh_meta,
          )
          return u, s, vh
  
imports: from __future__ import annotations
  import operator
  import warnings
  from functools import partial
  from numbers import Number
  import numpy
  import tlz
  from dask.array.core import Array, concatenate, dotmany, from_delayed
  from dask.array.creation import eye
  from dask.array.random import RandomState, default_rng
  from dask.array.utils import array_safe, meta_from_array, solve_triangular_safe, svd_flip
  from dask.base import tokenize, wait
  from dask.blockwise import blockwise
  from dask.delayed import delayed
  from dask.highlevelgraph import HighLevelGraph
  from dask.utils import apply, derived_from
label_desc: Use image node_id 5 for calling a global function with example usage: tsqr(data, compute_svd, _max_vchunk_size) and returns: q,  r, u,  s,  vh