Proposed resolution:

This wording is relative to N5014.

[Drafting note: The wording below implements option 1 of the issue discussion]

1. Modify [linalg.syn], header <linalg> synopsis, as indicated:

   namespace std::linalg {
     […]
     // [linalg.algs.blas1.ssq], scaled sum of squares of a vector's elements
     template<class Scalar>
       struct sum_of_squares_result {
         Scalar scaling_factor;
       };
     template<in-vector InVec, class Scalar>
       sum_of_squares_result<Scalar>
         vector_sum_of_squares(InVec v, sum_of_squares_result<Scalar> init);
     template<class ExecutionPolicy, in-vector InVec, class Scalar>
       sum_of_squares_result<Scalar>
         vector_sum_of_squares(ExecutionPolicy&& exec,
                               InVec v, sum_of_squares_result<Scalar> init); 
     […]
   }
   

2. Delete the entire [linalg.algs.blas1.ssq] as indicated:

   29.9.13.8 Scaled sum of squares of a vector's elements [linalg.algs.blas1.ssq]

   template<in-vector InVec, class Scalar>
     sum_of_squares_result<Scalar> vector_sum_of_squares(InVec v, sum_of_squares_result<Scalar> init);
   template<class ExecutionPolicy, in-vector InVec, class Scalar>
     sum_of_squares_result<Scalar> vector_sum_of_squares(ExecutionPolicy&& exec,
                                                               InVec v, sum_of_squares_result<Scalar> init);
   

   -1- [Note 1: These functions correspond to the LAPACK function xLASSQ[20]. — end note]

   -2- Mandates: decltype(abs-if-needed(declval<typename InVec::value_type>())) is convertible to `Scalar`.

   -3- Effects: Returns a value `result` such that

   1. (3.1) — `result.scaling_factor` is the maximum of `init.scaling_factor` and abs-if-needed(x[i]) for all `i` in the domain of `v`; and

   2. (3.2) — let `s2init` be

      init.scaling_factor * init.scaling_factor * init.scaled_sum_of_squares
      

      then `result.scaling_factor * result.scaling_factor * result.scaled_sum_of_squares` equals the sum of `s2init` and the squares of abs-if-needed(x[i]) for all `i` in the domain of `v`.

   -4- Remarks: If `InVec::value_type`, and `Scalar` are all floating-point types or specializations of `complex`, and if `Scalar` has higher precision than `InVec::value_type`, then intermediate terms in the sum use `Scalar`'s precision or greater.

3. Modify [linalg.algs.blas1.nrm2] as indicated:

   template<in-vector InVec, class Scalar>
     Scalar vector_two_norm(InVec v, Scalar init);
   template<class ExecutionPolicy, in-vector InVec, class Scalar>
     Scalar vector_two_norm(ExecutionPolicy&& exec, InVec v, Scalar init);
   

   -1- [Note 1: […] ]

   -2- Mandates: […]

   -3- Returns: […]

   [Note 2: […] ]

   -4- Remarks: […]

   [Note 3: An implementation of this function for floating-point types `T` can use the `scaled_sum_of_squares` result from `vector_sum_of_squares(x, {.scaling_factor=1.0, .scaled_sum_of_squares=init})`. — end note]

The current wording for `vector_sum_of_squares` [linalg.algs.blas1.ssq] has three problems with its specification of the value of `result.scaling_factor`.

1. The function permits `InVec::value_type` and `Scalar` to be any linear algebra value types. However, computing `result.scaling_factor` that satisfies both (3.1) and (3.2) requires more operations, such as division. Even if those operations are defined, they might not make `result.scaling_factor` satisfy the required properties. For example, integers have division, but integer division won't help here.

2. LAPACK's xLASSQ (the algorithm to which Note 1 in [linalg.algs.blas1.ssq] refers) changed its algorithm recently (see Reference-LAPACK/lapack/pull/#494) so that the scaling factor is no longer necessarily the maximum of the input scaling factor and the absolute value of all the input elements. It's a better algorithm and we would like to be able to use it.

3. Both members of sum_of_squares_result<Scalar> have the same type, `Scalar`. If the input `mdspan`'s `value_type` represents a quantity with units, this would not be correct. For example, if `value_type` has units of distance (say [m]), the sum of squares should have units of area ([m²]), while the scaling factor should have units of distance ([m]).

Problem (1) means that the current wording is broken. I suggest two different ways to fix this.

1. Remove `vector_sum_of_squares` entirely (both overloads from [linalg.syn], and the entire [linalg.algs.blas1.ssq]). That way, we won't be baking an old, less good algorithm into the Standard. Remove Note 3 from [linalg.algs.blas1.nrm2], which is the only other reference to `vector_sum_of_squares` in the Standard.

2. Fix [linalg.algs.blas1.ssq] by adding to the Mandates element (para 2) that `InVec::value_type` and `Scalar` are both floating-point types (so that we could fix this later if we want), and remove [linalg.algs.blas1.ssq] 3.1. Optionally add Recommended Practice, though Note 1 already suggests the intent.

I prefer just removing `vector_sum_of_squares`. Implementers who care about QoI of `vector_two_norm` should already know what to do. If somebody cares sufficiently, they can propose it back for C++29 and think about how to make it work for generic number types.