(Copied from an internal Mentor Graphics issue.)

Introduction
------------

C-VSIPL defines the following equation to compute the values y_n of a
minimum support 1-D convolution:

y_n = Sum k from 0 to M-1 of ( h_k * x_(n*D + (M-1) - k) )

Where:

N is the input vector length,
M is the kernel length,
D is the deimation,
x_j is the input vector (of size N)
h_k is the kernel vector (of size M)

C-VSIPL defines the length of the output vector (i.e. values n for
which y_n is defined) as:

floor( (N-1) / D ) - floor( (M-1) / D) + 1

In some cases, this equation requires the implementation to access
elements of x outside the x proper.

Example Illustrating the Problem
--------------------------------

To see this, first consider the case (where the equation is correct):

N = 5, M = 4, D = 1.

The expected output length is:

floor((5-1)/1) - floor((4-1)/1) + 1
= 4 - 3 + 1
= 2

Consider a diagram showing how each output is computed:

x_0 x_1 x_2 x_3 x_4

y_0 = h_3 h_2 h_1 h_0

y_1 = h_3 h_2 h_1 h_0

I.e.
y_0 = h_0*x_3 + h_1*x_2 + h_2*x_1 + h_3*x_0
y_1 = h_0*x_4 + h_1*x_3 + h_2*x_2 + h_3*x_1

In computing y_0 and y_1, only values

are required.

Next, consider the case:

N = 5, M = 4, D = 2.

The expected output length is:

floor((5-1)/2) - floor((4-1)/2) + 1
= floor(4/2) - floor(3/2) + 1
= 2 - 1 + 1
= 2

To compute y_n, we need to multiply h_k by x_(n D + (M-1) - k).
For n=1, k=0, the index to x is

n D + (M-1) - k
= 1 2 + (4-1) - 0
= 2 + 3 - 0
= 5

However, x_5 is not a valid element of x.

Pictorally:

x_0 x_1 x_2 x_3 x_4

y_0 = h_3 h_2 h_1 h_0

y_1 = h_3 h_2 h_1 h_0

Finally, consider two more cases:

N = 5, M = 4, D = 3 -> output length = 1

N = 5, M = 4, D = 4 -> output length = 2

(Why should D=3 be any different from the other cases?)

This does not make sense for several reasons:

• First, when D=2, a value (x_5) outside the support of x is required
  to compute the result. This might run counter to user's expectations
  of how the minimal output region of support should be defined.

• Second, changing the decimation

A New Equation for the Output Length
------------------------------------

What should the output length be?

Starting from the assumption that the formula for y_n is correct, we
need to determine the range of indices

for which all indices to x
are valid.

In particular, we want to find the largest n such that the following
holds for all k in 0 <= k <= M-1:

0 <= n*D + (M-1) - k < N

The lower bound is true for all n >= 0 since n*D >= 0 and (M-1)-k >= 0.

The upper bound:

n*D + (M-1) - k < N

The largest LHS occurs when k = 0:

n*D + (M-1) < N

n*D < N - M + 1

n < ceil((N-M+1) / D)

The largest n that we can compute y_n for is

n = ceil((N-M+1) / D)-1

Hence the output length is:

ceil((N-M+1) / D)

In the case where D=1, this reduces to:

N-M+1

(which is what the current equation reduces too when D=1).

Comparing the new equation with the old one for the previous examples:

N=5, M=4, D=1 old=2 new=2
N=5, M=4, D=2 old=2 new=1
N=5, M=4, D=3 old=1 new=1
N=5, M=4, D=4 old=2 new=1

Equation for Full Support Output Length
---------------------------------------

Similarly, the lengths for full support case are also incorrect.

For full support, the assumption is that each output vector should
have at least one value from x used to compute it:

n*D - k < N

(minimal index occurs when k = M-1)

n*D - (M-1) < N

n*D < N + (M-1)

n < ceil((N+M-1)/D)

However, here the error is less of an issue because the equation to
compute y_n is based on < x_k > rather than x_k. <x_k> is defined to
be 0 if k < 0 of k >= N

The rules for when a copy constructor "View v(w);" will create a view v
that aliases the view w are somewhat esoteric and likely to lead to
confusion – the result depends on details of the type of w, which do
not appear on this line and (with templates) may be in an entirely
separate section of the code.

It would be handy to have versions of this statement which are
predictable with regards to copy or alias, and which make the intent
explicit. Unfortunately this is not possible with the alias case, but
we can at least resolve it for the copy case.

Thus, I propose we add a deep_copy() function (or a method on views)
which would allow us to write "View v(deep_copy(w));" – with the result
that v has the same size as w, and the same values, but is guaranteed to
be a copy rather than an alias.

(Note that the "deep_copy" name is just a stand-in; there may be better
names for this functionality.)

This issue is copied from an internally-filed issue at Mentor that dates
back to the early days of the specification. Since that time, some
usage patterns have emerged that appear useful – most notably, the
implicit use of a length_type when a Domain<1> is required. This
particular case should either be preserved or additional overloads of
the relevant functions should be added.

However, in general this complaint is still quite relevant – for
instance, we should not be implicitly converting a length_type value to
a LU_solver object.

Note that I have not edited this to account for any changes to the
specification (either section numbers or normative changes) since the
issue was initially submitted in 2005.

In most cases, the resolution – adding an "explicit" keyword to many of
these constructors – seems straightforward and obvious. Chances are
very high that any existing code which is broken thereby was erroneous
in the first place.

Many of the constructors specified can take a single argument but are
not declared "explicit". A complete list:

Domain<1>(length_type)
Index<1>(index_type)
Dense<D>(Domain<D> const&)
const_Vector(Block&)
const_Vector(Vector const&)

• const_Vector<>(Vector const&)
  Vector<>(Block&)
  template <typename T0, typename Block0>
  Vector(Vector<T0, Block0> const&)
  const_Matrix<>(Block&)
• const_Matrix<>(Matrix const&)
  Matrix(Block&)
  template <typename T0, typename Block0>
  Matrix(Matrix<T0, Block0> const&)
  const_Tensor<>(Block&)
• const_Tensor<>(Tensor const&)
  Tensor(Block&)
  template <typename T0, typename Block0>
  Tensor(Tensor<T0, Block0> const&)
  lud<>(length_type) or LU_solver<>(length_type), see #29
  Rand<>(index_type)

Of these, only the starred conversions are defensible, and even
they will require some design adjustments. All the
rest are a problem; most are fixable simply by adding "explicit".
For the case of Domain<1>, this will require redesign of other
components to restore convenient usage; perhaps similarly for some
of the conversions from Block&.

Also, in the requirements list 6.3, the line View(Block &) should be
annotated to indicate that the conversion should not be implicit.

To expand on why these implicit conversions are dangerous...

First, implicit conversions from native numeric types are especially
problematic because of the promiscuous conversions among those types,
inherited from C, that occur even before calling the constructor.
As a result, any numeric type, whether int, float, char, or bool,
matches the constructor argument, and thus is accepted freely in
place of the class type by any function that takes an argument of
that type, regardless of whether it makes sense. These mistakes
can be very hard to spot; sometimes it results from calling f(a,b)
when f(b,a) was meant. (Functions declared to take numeric types
can sometimes be made safer by overloading with another that takes
a different, but acceptable, numeric type.)

Second, implicit conversions of any kind provoke annoying compiler
error messages complaining of ambiguities when more than one
conversion path is possible. The more implicit conversions that
are possible, the more frequently such ambiguities occur. The
workaround for such an ambiguity is to add a cast to the exact
type, but casts are themselves a source of errors, and besides
obscuring the logic of the code, they eliminate the convenience
the implicit conversion was supposed to provide.

Third, in cases where one conversion path is favored by the compiler
over another, the compiler may silently choose what the programmer
would consider the wrong one.

Fourth, template argument matching ignores conversions, making those
that do work less conveniently useful than they may seem.

The vsip_symmetry enumerator is used in a number of places, and in
general the VSIP_SYM_EVEN_LEN_* values are defined as "(Even) Symmetric,

length".

However, in the text there are sporadic mentions that imply that
symmetric here really means conjugate-symmetric; for instance, see the
bottom of page 567 ("The filter kernel can be even (conjugate) symmetric
or non-symmetric.") and the definition of the "kernel" argument to
vsip_*fir_create on page 612. This needs to be much more clearly and
consistently stated.

Regardless of whether the kernel is defined as symmetric or
conjugate-symmetric, the method for constructing the complete kernel
from the input in the symmetric cases is not defined. Is the input the
first half of the kernel, or the second half? Is the first half
conjugated, or the second half? This needs to be clearly specified.

This issue is copied from an internally-filed issue at Mentor that dates
back to the early days of the specification, circa 2005. I expect it
will be marked as "out of scope" for the VSIPL++ 1.2 FTF, but it should
be considered for a future 2.0 revision of the spec.

The specification requiring Blocks to implement reference counting
has several problems. First, it overconstrains implementations;
they may well have a better way to implement storage management,
if permitted. Second, as an interface it is clumsy and fails to
take advantage of C++ language strengths (most particularly,
destructors) for better safety. Finally, it is far from clear
that the library should ever dispose of a user block. If the
user created it, isn't the user's responsibility to dispose of
it when the user chooses, by whatever means the user chooses?

If in fact the library should take responsibility for discarding
a user's block when it is no longer needed, it suffices to provide
a way for the user to ask that ownership be transferred to the
library, without going into detail about how the library will
manage it afterward.

...

Note that this is also an exception-safety issue. In particular,
with the existing interface block storage could easily leak if an
exception were thrown before whoever was responsible to decrement
the count got to it.

If the spec were to say that all blocks derive from

template <bool isManaged, bool isConst,
typename valueT=..., typename mapT=...>
class vsip::Block_base;

then that base class and the views could handle all the refcounting
business invisibly. Besides being exception-safe and less messy,
it would be incrementally more convenient for users to define their
own block types. (Cf. std::iterator<>.)

Of course users could keep the library from deleting a managed block by
keeping their own view on it, too.

It is possible to implement the library so that blocks and programs
written to the old interface would still work as before.

The description of the ~Dense() destructor states "It should no longer
be used." The antecedent of "it" is grammatically unclear – although
users can guess from context, this should really say "The Dense object
should no longer be used."

It may perhaps be worth noting explicitly that any views bound to the
block should also not be used (although, in that case, what exactly does
"used" mean? Is it valid to access their size?).

The description of the ~Dense() destructor refers to the "reference
count". However, elsewhere the phrase "use count" is used. This should
be consistent

The description of how block reference-counting works with
view-allocated blocks is inconsistent.

In 6.1.1. [block.alloc], it is stated that "For a block allocated by a
view constructor, the effect of creation should include the effect of
invoking increment_count."

In 8.1 [view.view], however, the contrary requirement is given: "A
constructor that does allocate its block must not invoke its block's
increment_count member function."

These do not appear to be consistent – or, at the very least, are very
confusing. I suppose that the intended meaning is that, when a view
allocates a block, the act of creating that block should increment the
count, but that the view constructor should not otherwise increment the
count. However, that's not at all clear from the text.

Further, in 6.1.1, there is this note: "Note: For a block allocated by a
view constructor, the number of increment_count invocations should meet
or exceed the number of decrement_count invocations if the block is to
be used. When the number of decrement_count invocations exceeds, not
just equals, increment_count invocations, the block should be
deallocated." This seems at first glance to be inconsistent with the
definition of decrement_count, which states that the block will be
deallocated when the value reaches zero, not goes below zero. It is
only consistent if we assume that the phrase "the effect of invoking
increment_count" does not include the effect of incrementing the number
of times we consider increment_count to have been invoked!

Much of this would be clarified if we simply explicitly stated, in the
block requirements, that the "use count" should be set to one in the
block constructor, rather than the "effect of calling increment_count"
wording. Then we could omit the 6.1.1 paragraph 7 note entirely.

(That points out a related omission: We do not actually specify the
initial use-count value for the general case; for blocks that are not
Dense and not created by a view constructor, a conforming implementation
could initialize the use-count to infinity. Or three, just for
perversity. This is probably not desirable.)

In my previous issue, I stated with regards to clarifying the semantics
of a "View v(w)" statement: "Unfortunately this is not possible with the
alias case, but we can at least resolve it for the copy case.

This is not actually true. We could define an alias() function (or
method) on a View which returns a unspecified type, and then a View
constructor that accepts that type as input and either creates an alias
to the precursor view or throws a compile error. Thus, you could write

View v(alias(w));

and know that, if it compiles, you are guaranteed to get a true alias.

The Fir<> class is required to include both "static symmetry_type const
symmetry" and "symmetry_type symmetry()". This duplication of names is
illegal, and the duplication of functionality is pointless. (Thus, the
Mentor Graphics implementation provides only the const member, not the
method.)

Less illegally, the continuous_filter member and continuous_filtering()
method duplicate each other. Similarly, in [signal.convol] we have a
misspelled symmtry member and a symmetry() method, and in both
[signal.convol] and [signal.correl] we have a supprt member and
support() method. These duplicates and misspellings are confusing and
redundant, and should be removed.

The Fir<> copy constructor is declared "VSIP_NOTHROW". This is
erroneous and should instead be "VSIP_THROW((std::bad_alloc))", because
the copied Fir<> object must not share the storage allocated by the
source instance (in case the source is destroyed).

The same error applies to Convolution<>, Correlation<>, and potentially
others.

VSIPL++ implementations are not required to provide any member
definitions for tuple<>. However, both the old reference implementation
and the Mentor Graphics implementation define static members with the
values of the dimension ordering. These are certainly useful within the
implementations, and are likely to be useful to applications as well.

Recommendation: Require definition of 'dim0', 'dim1', and 'dim2' tuple
members, which match the tuple's template arguments.

There are several places where the VSIPL++ passes Domain<> arguments
where only the length portions of the domain are used. The Dense()
constructors are the most egregious cases, but there are others such as
Fftm(). This is confusing to users – we've noticed this in training
classes – and leads to unnecessarily initializing and passing around
the unused fields in the Domain objects.

It would be better to add a smaller object – e.g., Length<>, equivalent
to Index<> – that contains only the length values, and to redefine
these interfaces in terms of it. Implementations can continue to
support the Domain<> versions of the interfaces for backwards compatibility.

In most of these cases where N is fixed, it would also be an even
clearer (and unambiguous) interface to simply accept N length values as
integer arguments, and I suggest that those interfaces should be added
as well. We already effectively have them for the N=1 case because of
the implicit Domain<1>(int) constructor.

The effects of the decrement_count() accessor for blocks are described
as "Decrease the object’s use count. If the count becomes zero, the
block deallocates itself."

This is not always possible. Consider, for instance, this simple code:

Dense<...> block(...);
block.decrement_count();

Since block is created on the stack, it cannot deallocate itself.

Similar issues will occur with blocks created using placement-new into
existing allocated memory.

At minimum, this should be addressed by a clause stating that the
object's use-count shall not be decremented to zero (or allowed to
become decremented to zero) except for blocks that have been created
using "new".

It is clear from the existence of the increment_count and
decrement_count functions in the block requirements that the expectation
is that the reference-counting should apply to all VSIPL++ blocks.

However, there is no text at all that states this, aside from those two
function signatures.

There is no definition of "use count" anywhere. (See also my previous
issue noting that we need to explicitly state that the block constructor
initializes it to one.)

The only definition of the increment_count and decrement_count is in the
section describing the Dense block, although those definitions ought to
apply to all blocks.

This needs to be rearranged.

The specification for the SVD API is parametrized for the value-types to support real / complex single-precision / double-precision values.

However, the "decompose" method incorrectly requires the vector of singular values to use single-precision value-types only.

The specification defines complex_i as "The implementation shall define
this type to be complex<short>, complex<int>, or complex<long int> such
that the choice of underlying type is the same as scalar_i."

This is problematic, because the C++ spec states that complex<T> is not
necessarily defined when T is not float or double. In practice, GCC
provides only incomplete definitions for the int/short/long cases.

Since cscalar_i is only used in a very few places, I suggest that it
should simply be removed from the spec.

(The changelog in the back of the document mentions a change related to
this, but the actual text does not appear to be changed.)

The mag() and magsq() elementwise functions are defined to return real
views from complex views, as one would expect.

However, the maxmg, maxmgsq, minmg, and minmgsq reduction functions –
which ought to be simply the maximum/minimum of the mag and magsq
results, are defined to return complex values. They should likewise
return real values.

(Note that the Mentor Graphics implementation already does return real
values from these functions. The change is unlikely to break user code
in any case; if the user assigns the result to a complex variable, that
will still work when the functions return real values.)

The cmplx() function is placed at the end of section 10.3.5, which is
otherwise about "Scalar functions and their element-wise extensions" –
but cmplx() is defined only as an element-wise function, not as a scalar
function. It does not belong in this section.

In addition, the return value is given as "The element-wise extension of
the complex(T1, T2) constructor to views." This is not the return value
– it returns a view, not a constructor! This should be clarified.

The gather() and scatter() functions use an "Index<View<T,
Block0>::dim>" template to describe the dimension of the Indices that
they accept. However, "View" is not defined anywhere; the actual
input/output type is given as "const_View".

Thus, this should be changed to "Index<const_View<T, Block0>::dim>".

Copied from an internal bug report:

PROBLEM: Convolution accessor 'kernel_size()' is defined to return the
domain having the same length for each dimension as 'filter_coeffs'.
However, when constructing a convolution with a symmetric kernel
(sym_even_len_odd or sym_even_len_even), 'filter_coeffs' only holds
a subset of coefficients. In those cases, the true kernel size is
larder (For example, in 1D convolutions the true size will either '2 *
filter_coeffs.size()' or '2 * filter_coeff.size() + 1').

In contrast, the C-VSIPL spec defines the kernel size as M, and
specifies the size of 'filter_coeffs' as either M if symmetry =
non_sym, or 'floor(M/2)' if 'symmetry == sym_even_len_

'.
Queries to return the kernel size of a convolution return M.

PROPOSED FIX: Change wording of 'kernel_size()' accessor (and
Convolution constructor effects) to indicate correct kernel size.

TODO: Review specification of Correlation, FIR, and IIR for similar
error.

The VSIPL++ (and VSIPL) standards should define version macros, so that application code can be conditionalized on what version of the library they are being compiled with.

This should be done with numeric macros (i.e., VSIP_VERSION=1 and VSIP_VERSION_MINOR=2) to enable numeric comparisons.

In addition, there should be separate macros for VSIPL and VSIPL+, and either there should be separate version macros for VSIPL+ and Parallel VSIPL+, or there should be a macro ("VSIP_PARALLEL"?) to indicate whether the Parallel VSIPL+ pieces are supported.

The gather() and scatter() functions, as defined, are capable of
gathering from or scattering to multidimensional arrays, by virtue of
supporting multidimensional Index variables.

However, they cannot gather to or scatter from multidimensional
arrays, because they explicitly specify only Vectors of Index variables.
The extension to Matrices and Tensors of Index variables, with
corresponding-type gather() output and scatter() input, is obvious and
should be included.

There are only two functions in the entire VSIPL++ API that are required
to support computations on cscalar_i values: sumval and cumsum. This
seems quite inconsistent, and those requirements should probably be
removed for consistency.

(See also previous issue about removing cscalar_i entirely. The only
other uses of it are in view constructors, gather, and scatter.)

p48 says "Note to Implementors: memory allocation shall not be deferred." Why not? This is clearly an oblique reference to some unspoken semantics of the interface that isn't documented, and perhaps should be.

The ~Vector() destructor function is described as "If this object is the
only one using its block, the block is deleted. Otherwise, its block's
use count is decremented by one."

This is a redundant description of the decrement_block functionality,
and further spreads the description of the block reference-counting
throughout the standard.

It would be much better to either (a) explicitly state that "This calls
the block's decrement_count accessor, which may deallocate the block."

Or, if we do not want to constrain the implementation to literally call
the block method, we can alternately add a paragraph to the block
requirements section describing increment-count and decrement-count
functionality in general (with the increment_count and decrement_count
functions being interfaces to it), and the ~Vector description should
refer to that.

There are several oblique references to hidden interface semantics in "Notes to implementors", particularly the expectation that a development mode should keep reference counts on various data structures. However, since providing a development mode is not mandatory, this cannot be a normative part of the specification. This could be made clearer.