Current version of the draft is located here:
https://github.com/zertovitch/mathpaqs/blob/master/lin_alg/universal_matrices.ads

Dear Gautier,

thank you for the draft. Indeed something like that would be great possibly without polynomials or rational numbers etc. You just need doubles and/or floats and complex numbers. Sparse matrices are another animal, and once adopted they are another dimension of complexity. This is because apart from the triplet format which is easily expandable, one may also need to convert the triplet to CSR and there multiple entries with same row and column index need first to be summed so we have a single entry (i, j, value) for each (i, j) before converting to CSR. So I mean, one would like to build a library then and I am not so sure how easy would be to hardcode this library so that operations on matrices/vectors are as efficient as in Fortran. For Dense matrices it may still be possible to support mathematical expressions without introducing temporary objects. But if you also bring sparse into the game then?
(alpha,beta,gamma scalars, while x,y,z,w vectors, A, B_sparse matrices)

does this have a chance to be executed without temporaries?
y := alpha*x + beta*A*x + gamma*B_sparse*(z-w)

this has a chance to be executed without temporaries Fortran already does it.
y := alpha*x + beta*A*x

I would create a SparseMatrix library and break the statement above in 2
y1 := B_sparse*(z-w)
y := alpha*x + beta*A*x + gamma*y1;

I could survive with that. The problem dear Gautier is the temporary objects in such libraries. Very correctly J.P.Rosen suggested BLAS. The problem is that for simple vector operations (BLAS1) there is no speedup at all and the operations of BLAS1 cannot be executed for the whole expression, just in an `axpy (a*x+y)` manner one by one. This would require transversal of the arrays many times in a long expression involving vectors as (y := x1 + x2*alpha + beta*(x3-gamma*x4) + alpha*beta*(x1-x2)). For BLAS2 operations (matrix times vector) you may see a performance gain using BLAS but for large matrices, very large. Where you definitely see such a performance gain is for BLAS3 and there we definitely need something like that.

This discussion started about techniques to avoid generation of temporaries.. Even when you manipulate complex numbers you introduce temporaries. I need to come back to you guys with some number to understand the problem better. I will be back with these numbers. But until then, please let me know, do you feel that there is a need of a performance uplifting and mathematics enabling ADA library which is not Real_Arrays or Containers, but a real Math library allowing Real and Complex math manipulation and integrating possibly some Lapack routines as Real_Arrays do at the moment?

> > ADA should provide native support of vector-math or Vector.Numerics, or AdvancedNumerics or Numerics.LinearAlgebra and/or Numerics.SparseLinearAlgebra. Honestly with all the mission critical applications ADA is used for, I would expect something like that to enable people working in scientific computing and in general wherever Linear Algebra kernels are essential tool of the core algorithms development.
> Did you consider the facilities from Annex G? Moreover, there are
> bindings to BLAS, LAPACK, and more. Try "Ada mathematical libraries" on
> Google...

Thank you, very nice and positive suggestion. Please read my answer below to Gautier. Thank you again for the nice feedback.By the way, is there any other forum/discussion for ADA in github or elsewhere where we can also post code with colours etc?

Le 25/01/2023 à 10:52, Jim Paloander a écrit :
>
>>> ADA should provide native support of vector-math or Vector.Numerics, or AdvancedNumerics or Numerics.LinearAlgebra and/or Numerics.SparseLinearAlgebra. Honestly with all the mission critical applications ADA is used for, I would expect something like that to enable people working in scientific computing and in general wherever Linear Algebra kernels are essential tool of the core algorithms development.
>> Did you consider the facilities from Annex G? Moreover, there are
>> bindings to BLAS, LAPACK, and more. Try "Ada mathematical libraries" on
>> Google...
>
> Thank you, very nice and positive suggestion. Please read my answer below to Gautier. Thank you again for the nice feedback.By the way, is there any other forum/discussion for ADA in github or elsewhere where we can also post code with colours etc?

Gitter: https://gitter.im/ada-lang
LinkedIn: https://www.linkedin.com/groups/114211
Reddit: https://www.reddit.com/r/ada
StackOverflow: https://www.stackoverflow.com/questions/tagged/ada
Twitter: https://twitter.com/search?q=adaprogramming

--
J-P. Rosen
Adalog
2 rue du Docteur Lombard, 92441 Issy-les-Moulineaux CEDEX
https://www.adalog.fr https://www.adacontrol.fr

....and: https://forum.ada-lang.io/

> ...and: https://forum.ada-lang.io/

Which one do you think is better for general discussions on both ADA and ADA libraries?
Better in the sense that it is preferred by experienced open-minded people with engineering background or
people who may not necessarily hold an engineering degree or PhDs and Masters but who
have worked for many years on real world engineering problems ADA is supposed to assist.

On Sunday, January 22, 2023 at 2:34:19 PM UTC-7, ...@gmail.com wrote:
> Dear ADA lovers,
> with stack allocation of Real_Vector ( 1 .. N ) when N >= 100,000 I get STACK_OVERFLOW ERROR while trying to check how fast operator overloading is working for an expression

Leo's answer is responsive to the OP but seems to have gotten buried in an amazingly long discussion which I haven't read. This answer appeared here years ago (was it you, Leo?) and I have used it ever since. If there needs to be a long discussion maybe it could be about why this workaround is necessary. I will pull out Leo's answer from his post and put my version here:

type Real_Vector_Access is access Real_Vector;
then
x_Ptr : Real_Vector_Access := new Real_Vector(0 .. N - 1);
x : Real_Vector renames x_Ptr.all;

That's all. All the overloaded operators work as though the vector was declared from the stack. If you have code which formerly used stack allocation (of x, in this example), it will work without modification with this trick.

Jerry

> > Dear ADA lovers,
> > with stack allocation of Real_Vector ( 1 .. N ) when N >= 100,000 I get STACK_OVERFLOW ERROR while trying to check how fast operator overloading is working for an expression
> Leo's answer is responsive to the OP but seems to have gotten buried in an amazingly long discussion which I haven't read. This answer appeared here years ago (was it you, Leo?) and I have used it ever since. If there needs to be a long discussion maybe it could be about why this workaround is necessary. I will pull out Leo's answer from his post and put my version here:
>
> type Real_Vector_Access is access Real_Vector;
> then
> x_Ptr : Real_Vector_Access := new Real_Vector(0 .. N - 1);
> x : Real_Vector renames x_Ptr.all;
>
> That's all. All the overloaded operators work as though the vector was declared from the stack. If you have code which formerly used stack allocation (of x, in this example), it will work without modification with this trick.
>
> Jerry

Thank you Jerry, that's true, this trick does the job indeed. Thanks. But to be honest among us, it seems to me that the Real_Arrays library has been implemented probably (note the probably) by people not targeting real world engineering / industrial applications. Examples are Finite Volume solutions of semiconductor device equations, or weather prediction problems. You cannot expect to solve these on the stack and of course ADA containers are impractical for such computations. Possibly a library as discussed above that provides Real_Arrays functionality but memory allocation is done internally on the heap or even on the stack if the user chooses to do so to take advantage of fast stack access. But then the mechanism for changing the stack allocation size transparently and easily should also be provided as simply as setting the memory pool in an access type.

Do you agree?

In the meantime, I am wondering whether or not a hybrid approach (stack/heap) would work provided such a mechanism to set the stack size is easy accessible as simply as a memory pool. The functionality would be during operator overloading we create fast stack allocated arrays to transfer the intermediate results stored in temporary objects. So all the temporaries described above would reside in the stack if this would translate to performance improvement until a smarter approach is introduced similar to C++ expression templates.

> Thank you Jerry, that's true, this trick does the job indeed. Thanks. But to be honest among us, it seems to me that the Real_Arrays library has been implemented probably (note the probably) by people not targeting real world engineering / industrial applications. Examples are Finite Volume solutions of semiconductor device equations, or weather prediction problems. You cannot expect to solve these on the stack and of course ADA containers are impractical for such computations. Possibly a library as discussed above that provides Real_Arrays functionality but memory allocation is done internally on the heap or even on the stack if the user chooses to do so to take advantage of fast stack access. But then the mechanism for changing the stack allocation size transparently and easily should also be provided as simply as setting the memory pool in an access type.
>
> Do you agree?
>
Yes, with the minor caveat that I have no idea what I'm talking about.

When I ran into this limitation a few years ago, I learned how to increase the stack size before running my program. This helped but only a little, as large arrays still were not accommodated. (Define "large" as anything too big for stack expanded to OS limits.) I have never compared speed with stack versus heap allocation because, for me, it would serve no purpose. I can't just use smaller arrays in my work because large ones are too slow, if you follow my drift here.

The applications you mention are great but the list must be nearly endless. My applications run towards audio signals and radar (simulation) images.

Jerry