def lstsq(a, b, driver:'lsd', rcond:-1)
    a = NArray.asarray(a)
    b = NArray.asarray(b)
    b_orig = nil
    if b.shape.size==1
      b_orig = b
      b = b_orig[true,:new]
    end
    m = a.shape[-2]
    n = a.shape[-1]
    #nrhs = b.shape[-1]
    if m != b.shape[-2]
      raise NArray::ShapeError, "size mismatch: A-row and B-row"
    end
    if m < n   # need to extend b matrix
      shp = b.shape
      shp[-2] = n
      b2 = b.class.zeros(*shp)
      b2[false,0...m,true] = b
      b = b2
    end
    case driver.to_s
    when /^(ge)?lsd$/i
      # x, s, rank, info
      x, s, rank, = Lapack.call(:gelsd, a, b, rcond:rcond)
    when /^(ge)?lss$/i
      # v, x, s, rank, info
      _, x, s, rank, = Lapack.call(:gelss, a, b, rcond:rcond)
    when /^(ge)?lsy$/i
      jpvt = Int32.zeros(*a[false,0,true].shape)
      # v, x, jpvt, rank, info
      _, x, _, rank, = Lapack.call(:gelsy, a, b, jpvt, rcond:rcond)
      s = nil
    else
      raise ArgumentError, "invalid driver: #{driver}"
    end
    resids = nil
    if m > n
      if /ls(d|s)$/i =~ driver
        case rank
        when n
          resids = (x[n..-1,true].abs**2).sum(axis:0)
        when NArray
          resids = (x[false,n..-1,true].abs**2).sum(axis:-2)
          ## NArray does not suppurt this yet.
          # resids = x[false,0,true].new_zeros
          # mask = rank.eq(n)
          # resids[mask,true] = (x[mask,n..-1,true].abs**2).sum(axis:-2)
        end
      end
      x = x[false,0...n,true]
    end
    if b_orig && b_orig.shape.size==1
      x = x[true,0]
      resids &&= resids[false,0]
    end
    [x, resids, rank, s]
  end
def lstsq(a, b, driver:'lsd', rcond:-1)
    a = NArray.asarray(a)
    b = NArray.asarray(b)
    b_orig = nil
    if b.shape.size==1
      b_orig = b
      b = b_orig[true,:new]
    end
    m = a.shape[-2]
    n = a.shape[-1]
    #nrhs = b.shape[-1]
    if m != b.shape[-2]
      raise NArray::ShapeError, "size mismatch: A-row and B-row"
    end
    if m < n   # need to extend b matrix
      shp = b.shape
      shp[-2] = n
      b2 = b.class.zeros(*shp)
      b2[false,0...m,true] = b
      b = b2
    end
    case driver.to_s
    when /^(ge)?lsd$/i
      # x, s, rank, info
      x, s, rank, = Lapack.call(:gelsd, a, b, rcond:rcond)
    when /^(ge)?lss$/i
      # v, x, s, rank, info
      _, x, s, rank, = Lapack.call(:gelss, a, b, rcond:rcond)
    when /^(ge)?lsy$/i
      jpvt = Int32.zeros(*a[false,0,true].shape)
      # v, x, jpvt, rank, info
      _, x, _, rank, = Lapack.call(:gelsy, a, b, jpvt, rcond:rcond)
      s = nil
    else
      raise ArgumentError, "invalid driver: #{driver}"
    end
    resids = nil
    if m > n
      if /ls(d|s)$/i =~ driver
        case rank
        when n
          resids = (x[n..-1,true].abs**2).sum(axis:0)
        when NArray
          resids = (x[false,n..-1,true].abs**2).sum(axis:-2)
          ## NArray does not suppurt this yet.
          # resids = x[false,0,true].new_zeros
          # mask = rank.eq(n)
          # resids[mask,true] = (x[mask,n..-1,true].abs**2).sum(axis:-2)
        end
      end
      x = x[false,0...n,true]
    end
    if b_orig && b_orig.shape.size==1
      x = x[true,0]
      resids &&= resids[false,0]
    end
    [x, resids, rank, s]
  end
Computes the minimum-norm solution to a linear least squares
 problem:

         minimize 2-norm(| b - A*x |)

 using the singular value decomposition (SVD) of A.
 A is an M-by-N matrix which may be rank-deficient.
 @param a [Numo::NArray] m-by-n matrix A (>= 2-dimensinal NArray)
 @param b [Numo::NArray] m-by-nrhs right-hand-side matrix b
   (>= 1-dimensinal NArray)
 @param driver [String or Symbol] choose LAPACK driver from
   'lsd','lss','lsy' (optional, default='lsd')
 @param rcond [Float] (optional, default=-1)
   RCOND is used to determine the effective rank of A.
   Singular values `S(i) <= RCOND*S(1)` are treated as zero.
   If RCOND < 0, machine precision is used instead.
 @return [[x, resids, rank, s]]
   - **x** -- The solution matrix/vector X.
   - **resids** -- Sums of residues, squared 2-norm for each column in
     `b - a x`. If matrix_rank(a) < N or > M, or 'gelsy' is used,
     this is an empty array.
   - **rank** --  The effective rank of A, i.e.,
     the number of singular values which are greater than RCOND*S(1).
   - **s** --  The singular values of A in decreasing order.
     Returns nil if 'gelsy' is used.