I’m implementing a linear map. This operator basically applies three matrices one after the other, to a vector. Example : you have an SPD matrix A, and you have L its incomplete Cholesky factor, and I want to have an entity PrecondA that is multipliable, such that :
PrecondA * v = L \ ( A * ( L’ \ v ) )
So what I did was defining a structure PrecondMatrix, with member fields A and L, and I overloaded the * operator.
Now let’s say I want to implement a CG method. I create a function CGSolve(A::Matrix, b…). Let’s say now that I want to have a CG method that uses the preconditioned system instead of the original system. Then, at the moment, I’m forced to create a new function CGSolve(A::PrecondMatrix, b …) which has the exact same implementation, since I overloaded the * operator for my PrecondMatrix structure !
Is there a way of avoiding this ? If I declare my PrecondMatrix as an AbstractArray, I have to actually give its coefficients, which requires computation, I don’t wan’t to do this. I would prefer an abstract type that basically says : “this thing is a vector multiplier, its application to a vector is provided. It is not known by its coefficients”. And both Array{Float, 2} and PrecondMatrix whould implement this interface.
I hope I was clear enough. Thanks in advance for your help!
Are you referring to some specific implementation available in some package here?
In principle, if the operations inside the function are compatible with the input, you don’t need to have the new type as a subtype of anything, just drop the Matrix from A::Matrix (that is, leave it untyped), or put Union{Matrix,PrecondMatrix} instead.
Hi, thanks for your quick reply!
Sorry if I wasn’t clear, I was not referring to specific implementation, just any naïve implementation of CG method. It was just an example of situation. I have to admit that your answer covers the troubles I experienced so far. But as to the general question : there isn’t linear algebra interfaces ?
There are no general “interfaces” in Julia (lots of discussions about them). Specific functions, even within the umbrella of “Linear Algebra” may required different properties from the input. This is sometimes specified on dispatch (like different methods for diagonal or sparse arrays), sometimes if the input does not satisfy the expected requirements you will get an error down the line (like if the * operation was not defined for the elements of the matrix, the error will appear the first time that is tried).
If you are implementing your CG algorithm, you can well take care of it being generic for the types you are willing to use. Or, if you think the type you are creating is important enough, you may even test the CG algorithm implemented in some well maintained package, and see if you can contribute to it by making it generic to your type as well (if it isn’t already).