General questions about linalgcpp

[chakshinglee]

linalgcpp/include/vectorview.hpp

Lines 21 to 31 in 6186c09

	/*! @brief Vector view of data and size
	@note Views are only modifiable if you own the view
	that you plan to change it to.
	If you want view A to be equal to view B,
	you most own both A and B. Otherwise,
	it is trivial to subvert const restrictions.
	I'm not sure if this is good way or not to
	deal w/ this.
	*/

How do to determine if a view is owned?

By the way, should "most" be "must" in the documentation?

[gelever]

Perhaps this documentation is poorly worded. You can only create a copy of a view if the original view is not const. This prevents modifying a const view by simply copying it into a new non-const view (but the data they both view is the same).

Yes, good catch on the typo.

[chakshinglee]

linalgcpp/include/vectorview.hpp

Lines 616 to 622 in 6186c09

	/*! @brief Check if two vectors are equal
	@param lhs left hand side vector
	@param rhs right hand side vector
	@retval true if vectors are close enough to equal
	*/
	template <typename T, typename U>
	bool operator==(const VectorView<T>& lhs, const VectorView<U>& rhs)

Would it better to have "exact equal" for this operator (that's what I think of when I see ==)? If someone just want to measure if two vectors are close enough, they can check something like norm(v1-v2) < tol. This way the user knows what they are doing and has more flexible control?

[gelever]

Yes, that would definitely be an option for vector operator equal. I think that would imply that lhs and rhs are pointing to the same data. Since, if the vectors are non-trivial, the chances of them having the exact floating point entries is pretty slim.

[chakshinglee]

(Apparently I cannot reply to a comment right below it?)

For the ownership of VectorView, I still don't understand. So VectorView by definition does not own the data, right? The ownership here means const?

[chakshinglee]

linalgcpp/include/vectorview.hpp

Lines 247 to 257 in 6186c09

	template <typename T>
	VectorView<T>& VectorView<T>::operator=(VectorView<T>&& vect) noexcept
	{
	assert(vect.size_ == size_);
	std::copy(std::begin(vect), std::end(vect), std::begin(*this));
	//data_ = vect.data_;
	//size_ = vect.size_;
	return *this;
	}

Why the move assignment needs to copy the data? It is not safe the move the data?

[gelever]

Correct, VectorView does not own any data.

Operators are defined on the underlying data, not the view itself. I originally had it the other way around. In practice, this is much more useful. Because of the change, ownership now is not really an issue.

[chakshinglee]

linalgcpp/include/vectorview.hpp

Lines 476 to 481 in 6186c09

	/*! @brief Entrywise multiplication x_i = x_i * y_i
	@param lhs left hand side vector x
	@param rhs right hand side vector y
	*/
	template <typename T>
	VectorView<T> operator*=(VectorView<T> lhs, const VectorView<T>& rhs)

I found this function counter-intuitive. I understand what the function is doing, just to me mathematically there are two natural multiplications for vectors, one is vector * scalar = vector, one is vector * vector = scalar (of course, this is just my own intuition). I totally agree this is a useful function, I just dont find vector1 *= vector2 easy to follow.

[chakshinglee]

If by definition VectorView does not own data, then what is the meaning of "if you own the view" in the documentation?

[gelever]

Ownership refer(ed) to the VectorView itself. I will update the docs to clarify the restrictions on VectorView

[chakshinglee]

linalgcpp/include/sparsematrix.hpp

Lines 344 to 347 in 6186c09

	/*! @brief Scale rows by diagonal matrix
	@param values scale per row
	*/
	void ScaleRows(const SparseMatrix<T>& values);

For a general linear algebra library, this function does not make much sense to me (although I see how it can be used in smoothG). And in general if values is a non-diagonal matrix (with nnz = number of rows), the function will not complain, so it can be hard to debug if a non-diagonal matrix is passed in accidentally.

The same comment applies to three other similar functions.

[chakshinglee]

linalgcpp/src/blockoperator.cpp

Line 91 in 6186c09

x_ = input;

linalgcpp/src/blockoperator.cpp

Line 114 in 6186c09

output = y_;

Do these assignments make deep copies? If so, is it possible to avoid it?

[gelever]

For scaling, I think those are useful enough to belong in the linear algebra library. I agree that scaling by a matrix may be dangerous and requires either more clear documentation or stricter debug checks.

Yes, those deep copy for now, until I create BlockVectorView. It is a temp solution until I get there.

[chakshinglee]

I agree the functionality of the function itself is useful. I just cannot see what does it mean by scaling rows of a matrix by another matrix (thinking from the point of reading some code that uses this function).

By the way, how do you think about making a special matrix class for diagonal matrices? This needs not be in linalgcpp (as the class is probably too restrictive in general to be used), but can be created in smoothG (derived from Operator or something).

[chakshinglee]

linalgcpp/include/sparsematrix.hpp

Line 342 in 6186c09

T MaxNorm() const;

MaxNorm is not describing what it actually computes? This wiki page https://en.wikipedia.org/wiki/Matrix_norm calls it "entrywise" L1 norm (which is different from the induced L1 norm). Of course, that's just a reference.

[gelever]

Ah ok, I thought you meant scaling matrices in general. Yes, for scaling by matrix, the function name is misleading.

For matrix norm, you are correct. It is incorrectly named. I needed something that would determine if a matrix is close enough to 0 and named it w/o looking up the proper name. I'll modify it.