Real_Arrays on heap with overloaded operators and clean syntax

[Jim Paloander]

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

X := A + B + C + C + A + B, where
A,B,C,X are all Real_Vector ( 1 .. N ).

So my only option was to allocate on the heap using new. But then I lost the clean syntax

X := A + B + C + C + A + B

and I had to write instead:

X.all := A.all + B.all + C.all + C.all + A.all + B.all.

This is really ugly and annoying because when you are using Real_Arrays for implementing some linear algebra method who relies heavilly on matrix vector products and vector updates, you do need to allocate on the heap (sizes are determined in runtime) and you do need a clean syntax. So, is there any way to simplify my life without using the .all or even without declaring A,B,C,X as access Real_Vector?
Thanks for your time!

[Joakim Strandberg]

Easiest solution is probably to declare a new task and specify the stack size using the Storage_Size aspect. Allocate as much stack space as you need to be able to do the calculations and do all the allocations on the declared task, not on the environment task. You will avoid the unnecessary heap allocations and have nice clean syntax.

Best regards,
Joakim

[Jim Paloander]

Thank you for your reply,
since I am a newbie I was under the impression that tasks are used only when you want to write a parallel code that takes advantage of multicore architectures. You suggest I have a single task and single thread something like this? I see, but there should be a way to do this also for the main program. But thanks anyway. Are you aware of any libraries similar to Real_Arrays, but who allocated memory internally using heap? This is the natural way to do such things. Similarly to the Containers.Vector. But Vector has such an awful syntax. There should be something like an indexer [i] similarly to the C++ std::vector to make things simpler and overloaded operators similarly to Real_Arrays. It is a no brainer. Most programs need to allocate on the heap, why did they restrict Real_Arrays on the stack?

[Dmitry A. Kazakov]

You can define "+" on the access type, which should probably be an arena
pointer for performance reasons:

Arena : Mark_And_Release_Pool;
type Real_Vector_Ptr is access Real_Vector;
for Real_Vector_Ptr'Storage_Pool use Arena;

function "+" (Left, Right : Real_Vector_Ptr)
return Real_Vector_Ptr is
begin
if Left'Length /= Right'Length then
raise Constraint_Error;
end if;
return Result : Real_Vector_Ptr := new Real_Vector (Left'Range) do
for I in Result'Range loop
Result (I) :=
Left (I) + Right (I - Left'First + Right'First);
end loop;
end return;
end "+";

You can overload that with

function "+" (Left : Real_Vector_Ptr; Right : Real_Vector)
return Real_Vector_Ptr is
begin
if Left'Length /= Right'Length then
raise Constraint_Error;
end if;
return Result : Real_Vector_Ptr := new Real_Vector (Left'Range) do
for I in Result'Range loop
Result (I) :=
Left (I) + Right (I - Left'First + Right'First);
end loop;
end return;
end "+";

and with

function "+" (Left : Real_Vector; Right : Real_Vector_Ptr)
return Real_Vector_Ptr;

Then you will be able to write:

X := A + B + C + C + A + B;
-- Use X
Free (X); -- Pop all arena garbage

But of course, the optimal way to work large linear algebra problems is
by using in-place operations, e.g.

procedure Add (Left : in out Real_Vector; Right : Real_Vector);

etc.

--
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

[Jim Paloander]

Privet Dmitry,
thank you for your answer, I was thinking about that, but I was not sure whether or not it can be avoided with Implicit_Dereference,

type Accessor (Data: not null access Element) is limited private
with Implicit_Dereference => Data;

somehow... Otherwise what you described for operator+ one has to do for every operator overloaded inside Real_Arrays package. The optimal way to work large linear algebra problem is what you describe because unfortunately ADA does not allow what Fortran does since 30 years ago or more. But in C++ you can reproduce the same functionality as Fortran using Expression Templates and Template Metaprogramming. Perhaps ADA should allow something like that. Because for maintainability reasons the best would be to write the mathematical expressions as close as possible to the mathematical formulas.

Spacebo Bolshoe,
Jim.

[Joakim Strandberg]

It my impression that in the Ada community the preferred way of working is in general stack only. Heap allocations are avoided for a number of reasons for example performance, the application needs to ask the operating system for memory which one doesn't know how much time that will take nor if it always will succeed. In addition, applications that use the heap may be susceptible to heap memory fragmentation. In Ada, it is easy to specify stack sizes when declaring tasks. It is not part of the Ada standard to specify stack size of the environment task. The Ada way is to declare new tasks and do work on them. Prefer the bounded containers over the unbounded containers. If you really need to allocate objects I recommend using storage pools with pre-allocated memory also known as arena pools.

[Jim Paloander]

With great depression I realized that the preferred way is of stack only. This is very restrictive excluding all scientific modelling involving solvers for partial differential equations, linear algebra kernels, etc. It is insane. Completely insane. 3D simulations of physical phenomena may involve billions of grid-cells and at each grid-cell several unknowns are defined (velocity, pressure, temperature, energy, density, etc). That is why they are using Fortran or C++, but ADA has really cool stuff for so many things, why not vectors and matrices and heap allocation? Would you please give me an example, I googled and I cannot find a single example demonstrating how to use a task with the declaration of stack size. Why is there so little information online about so important things such as allocation?

[Joakim Strandberg]

For what you described above being able to write
X := A + B + C + C + A + B;
working with a bigger stack is the easiest solution. Don't understand why you think it would restrict all scientific modelling? Maybe I haven't described it well enough.

From the top of my head:
task T with Storage_Size => 10_000_000 is

end T;

task body T is
begin
null;
end T;

Tasks need to be defined inside Ada packages, they cannot be stand-alone.

An optional type is for example in Ada:
type Optional_Integer (Exists : Boolean := False) is record
case Exists is
when True => Value : Integer;
when False => null;
end Exists;
end record;

[Gautier write-only address]

Note that Real_Arrays does not specify where things are allocated (heap or stack).
Only when you define "x : Real_Vector (1 .. n)", it is on stack. You can always write something like the snippet below.
Anyway, after a certain size, you may have to find compromises, like avoiding operators (they do too many allocations & deallocations in the background, even assuming elegant heap-allocated objects) and also give up plain matrices, against sparse matrices or band-stored matrices, typically for solving Partial Differential Equations.

with Ada.Numerics.Generic_Real_Arrays;

procedure Test_Large is

type Float_15 is digits 15;

package F15_R_A is new Ada.Numerics.Generic_Real_Arrays (Float_15);

use F15_R_A;

procedure Solve_it
(x : in Real_Vector;
y : out Real_Vector;
A : in Real_Matrix) is
begin
null; -- Here, the big number-crunching
end;

n : constant := 10_000;

type Vector_Access is access Real_Vector;
type Matrix_Access is access Real_Matrix;

x, y : Vector_Access := new Real_Vector (1 .. n);
A : Matrix_Access := new Real_Matrix (1 .. n, 1 .. n);

begin
Solve_it (x.all, y.all, A.all);
-- !! Deallocation here
end;

[Rod Kay]

If you are on linux, then you could set the stack size with

$ ulimit -s unlimited
$ launch_my_app



Regards.

[Jim Paloander]

On Windows 10 with mingw64?

[Rod Kay]

Not sure. I don't have a windows machine.

What happens when try ?

$ ulimit -a

[Joakim Strandberg]

Something came up and I had to send my previous reply/e-mail as is. I wanted to find the video where Jean Pierre Rosen talks about how memory is handled in the Ada language from FOSDEM perhaps 2018-2019. Unfortunately I have been unable to find it.

Secondly, example of Ada code where Storage pool usage is demonstrated: https://github.com/joakim-strandberg/advent_of_code

Example of task with rendez-vous mechanism:

task T with Storage_Size => 10_000_000 is

entry Run (I : Integer);

end T;

task body T is
begin

loop
begin
select
accept Run (I : Integer) do
null; -- Here save the value of the integer I for later processing in the task
end Run;
null; -- Whenever another task calls Run on this task, the work to be done should be put here.
-- Put the mathematical calcuations here.
or
terminate; -- To end the select-statement with "or terminate;" means the task will terminate when the environment task has finished execution of the "Main" procedure of the application. No need to tell the task T that now it is time to shutdown.
end select;
exception
when Error : others =>
-- Print error to standard out here using subprograms from Ada.Exceptions and Ada.Text_IO.
end;
end loop;
end T;

Best regards,
Joakim

[Leo Brewin]

Here is a slight variation on the solution suggested by Gautier. It uses
Aad's "rename" syntax so that you can avoid all the .all stuff. I use
this construction extensively in my large scale scientific computations.

with Ada.Numerics.Generic_Real_Arrays;
with Ada.Unchecked_Deallocation;

procedure Test_Large is

type Float_15 is digits 15;

package F15_R_A is new Ada.Numerics.Generic_Real_Arrays (Float_15);

use F15_R_A;

procedure Solve_it
(x : in Real_Vector;
y : out Real_Vector;
A : in Real_Matrix) is
begin
null; -- Here, the big number-crunching
end;

n : constant := 10_000;

type Vector_Access is access Real_Vector;
type Matrix_Access is access Real_Matrix;

x_ptr, y_ptr : Vector_Access := new Real_Vector (1 .. n);
A_ptr : Matrix_Access := new Real_Matrix (1 .. n, 1 .. n);

x : Real_Vector renames x_ptr.all;
y : Real_Vector renames y_ptr.all;
A : Real_Matrix renames A_ptr.all;

procedure FreeVector is new
Ada.Unchecked_Deallocation (Real_Vector,Vector_Access);
procedure FreeMatrix is new
Ada.Unchecked_Deallocation (Real_Matrix,Matrix_Access);

begin
Solve_it (x, y, A);
-- Deallocation here
FreeVector (x_ptr);
FreeVector (y_ptr);
FreeMatrix (A_ptr);
end;

[Jim Paloander]

Thank you very much, would a for Real_Vector_Access'Storage_Pool use localPool; save you from the need to free the vectors and matrix yourself?

On the other hand, is there any way of avoiding temporaries? Possibly a modern version of the Real_Array using expression generic syntax or something similar? Since you are using scientific computationas extensively, you must be aware of Fortran. Have you compared Fortran's complex numbers with ADA's for inner products or similar computations to see who is faster? You see, I like a lot of things about ADA, but the syntax is really difficult to follow. Sometimes it gives me the impression that it is more difficult than really needed to be. For example there should be a way for Real_Arrays to allocate memory internally and not on the stack directly like containers. And for containers to provide an indexer Vector(i) and overloaded operators similarly to Real_Vectors. But the fact that they do not give me the impression that this Language although being designed by the army for mission critical applications never realized that modern mission critical need to simplify mathematical calculations providing an easy syntax. I am surprised that after so many years and so many updates to the Standard the designers of the Language did not realized that such mathematical syntax should be simplified to attract more people from scientific computing, who are tired with Fortran 10000 ways of declaring something a variable.

[Rod Kay]

ulimit is available on cygwin.

It is not available on mingw64 then ?

[Jim Paloander]

It is, but I am not sure if it works since -D400m -d400m together with ulimit were not able to solve the problem. So I thought that under windows 10 and mingw things related to stack size settings are not equivalent to Linux. I checked your repository. Looks good. We need some elegant solution, because sometimes it is not easy to predict the stack size that will be needed. Imagine you write a commercial application that is running from a GUI and the user loads a particular case from a file. The memory that will be needed may even exceed the available memory installed on the computer. There must be some elegant way to set the stack size globally. On the other hand storage pool seem that does not need deallocation. Isnt there any storage pool for stack based allocation so that you enjoy best of both worlds?

[Rod Kay]

>>>>> If you are on linux, then you could set the stack size with
>>>>>
>>>>> $ ulimit -s unlimited
>>>>> $ launch_my_app
>>>>>
>>>>>
>>>>>
>>>>> Regards.
>>>> On Windows 10 with mingw64?
>>>
>>> Not sure. I don't have a windows machine.
>>>
>>> What happens when try ?
>>>
>>> $ ulimit -a
>> ulimit is available on cygwin.
>>
>> It is not available on mingw64 then ?
> It is, but I am not sure if it works since -D400m -d400m together with ulimit were not able to solve the problem. So I thought that under windows 10 and mingw things related to stack size settings are not equivalent to Linux. I checked your repository. Looks good. We need some elegant solution, because sometimes it is not easy to predict the stack size that will be needed. Imagine you write a commercial application that is running from a GUI and the user loads a particular case from a file. The memory that will be needed may even exceed the available memory installed on the computer. There must be some elegant way to set the stack size globally. On the other hand storage pool seem that does not need deallocation. Isnt there any storage pool for stack based allocation so that you enjoy best of both worlds?

Afraid I'm not familiar with '-D400m -d400m'. I did see a few examples
on the net which seemed to set stack size as a linker parameter.

My repository ? Which one ?


Regards.

[Egil H H]

On Monday, January 23, 2023 at 12:53:15 AM UTC+1, joak...@kth.se wrote:
> Something came up and I had to send my previous reply/e-mail as is. I wanted to find the video where Jean Pierre Rosen talks about how memory is handled in the Ada language from FOSDEM perhaps 2018-2019. Unfortunately I have been unable to find it.
>

It was in 2016:

https://archive.fosdem.org/2016/schedule/event/ada_memory/

[Dmitry A. Kazakov]

On 2023-01-22 23:36, Jim Paloander wrote:

> I was not sure whether or not it can be avoided with Implicit_Dereference,
>
> type Accessor (Data: not null access Element) is limited private
> with Implicit_Dereference => Data;

If you create a new wrapper type, anyway, then it is easier to define
operations directly on that new type.

> Otherwise what you described for operator+ one has to do for every operator overloaded inside Real_Arrays package.

You should not use the standard library anyway. It is not intended for
large problems, which require specific approaches and methods, like
sparse matrices, concurrent processing and so on.

> The optimal way to work large linear algebra problem is what you describe because unfortunately ADA does not allow what Fortran does since 30 years ago or more.

I am not sure what you mean. It is quite possible to design a wrapper
datatype allocating vectors/matrices in the pool. E.g. Ada's
Unbounded_String is such a thing. Real_Arrays were not designed this way
because see above.

> But in C++ you can reproduce the same functionality as Fortran using Expression Templates and Template Metaprogramming.

Nothing prevents you from wrapping Real_Array in a generic way:

generic
with package Real_Arrays is new Numerics.Generic_Real_Arrays (<>);
package Generic_Pool_Real_Arrays is
...
end Generic_Pool_Real_Arrays;

> Perhaps ADA should allow something like that. Because for maintainability reasons the best would be to write the mathematical expressions as close as possible to the mathematical formulas.

There is no problem with that as you can define operations on pointers.

[G.B.]

On 22.01.23 23:07, Jim Paloander 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
>>>
>>> X := A + B + C + C + A + B, where
>>> A,B,C,X are all Real_Vector ( 1 .. N ).
>>>
>>> So my only option was to allocate on the heap using new. But then I lost the clean syntax
>>>
>>> X := A + B + C + C + A + B
>>>
>>> and I had to write instead:
>>>
>>> X.all := A.all + B.all + C.all + C.all + A.all + B.all.
>>>
>>> This is really ugly and annoying because when you are using Real_Arrays for implementing some linear algebra method who relies heavilly on matrix vector products and vector updates, you do need to allocate on the heap (sizes are determined in runtime) and you do need a clean syntax. So, is there any way to simplify my life without using the .all or even without declaring A,B,C,X as access Real_Vector?
>>> Thanks for your time!
>> Easiest solution is probably to declare a new task and specify the stack size using the Storage_Size aspect. Allocate as much stack space as you need to be able to do the calculations and do all the allocations on the declared task, not on the environment task. You will avoid the unnecessary heap allocations and have nice clean syntax.
>>
>> Best regards,
>> Joakim
>
> Thank you for your reply,
> since I am a newbie I was under the impression that tasks are used only when you want to write a parallel code that takes advantage of multicore architectures.

Tasks---their name is an indication that they just
perform something that is specified---are not restricted to
parallelism, or to multiple cores. The can do the job on
a single processor, or move to some other, they are more
abstract than contemporary hardware, or hardware for that
matter.

> You suggest I have a single task and single thread something like this? I see, but there should be a way to do this also for the main program.

It is unfortunate that GCC had to introduce a tiny scale
stack(frame) model some time ago.

> But thanks anyway. Are you aware of any libraries similar to Real_Arrays, but who allocated memory internally using heap? This is the natural way to do such things.

The most natural way to work with an array of FPT numbers is
for the programmer to declare an array indexed by some index type.
Done. If GNAT gets in the way there, it might be worth a note sent to
its maintainers. Whenever a programmer is tasked with considering
memory allocation, then depending on one's propensity towards working
on memory allocation it is inconvenient and distracting.
Math programs don't make you do this, I think.

Also, std::vector an its relatives shield the programmer from
the absurdly clever, yet unreadable memory allocation that
needs to be stuffed behind the scenes.
More importantly, though, C++ introduced std::move semantics
after a few decades of its existence, to address copying when
using chains of +. It might be interesting to see Ada's in-situ
construction of return values in comparison.

> Similarly to the Containers.Vector. But Vector has such an awful syntax. There should be something like an indexer [i] similarly to the C++ std::vector to make things simpler

.at() does some of what Ada does.
Is

v.at(k) = 4;

less awful than

v(k) := 4;

?

Another thing: Mathematical notation has ellipsis, thus

A + B + ... + Y + Z;

Most general purpose languages don't have ellipsis for this kind
of expression. However, even mathematical formulas use what
programmers can usually achieve, too. The usual

\sum_k A_k.

No "+" at all, and an array of vectors, not single ones.
Going further, some like to write

reduce("+", A);

In Ada, you could have a generic function for this, or use a
function pointer.

The .all thing vanishes automatically whenever you want to refer
to a particular component of the pointed-at object, as opposed
to all of them. So, A.all(K) is the same as A(K).
Likewise, .all can be dropped if want to invoke the pointed-at
subprogram if it has parameters.

[J-P. Rosen]

Thanks Egil, you were faster than me...
I have also a full tutorial at several Ada-Europe conferences. No video,
but I can send the slides to those interested.

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

[Jim Paloander]

Dmitry,
you mentioned in a previous email

> You should not use the standard library anyway. It is not intended for
> large problems, which require specific approaches and methods, like
> sparse matrices, concurrent processing and so on.
> > The optimal way to work large linear algebra problem is what you describe because unfortunately ADA does not allow what Fortran does since 30 years ago or more.
> I am not sure what you mean. It is quite possible to design a wrapper

that the optimal way is to work with specific approaches and methods. Nevertheless, at least for vector operations, and/or dense matrices operations, 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.

So it seems that as you correctly implied even dense linear vector math are not enough. But it is at least what I hoped for. Sparse matrix times vector could be stored in a temporary vector. I would not expect that an expression such as

y := alpha*a + beta*(z-w*gamma) + alpha*beta*SparseMat*x;

is computed in a way similar to Fortran, which is

for I in 1 .. N loop
y[I] := alpha*a[I] + beta*(z[I] - w[I]*gamma) + alpha*beta*Sum_j ( SparseMat[I,j]*x[j])
end

where Sum_j is implemented in another loop. I would instead compute v := SparseMat*x, and then replace SparseMat*x with v in the above expression. So I would be happy if ADA would give me a library that allows Vector Math with performance similar to Fortran without temporaries for both Complex and Real numbers. If they could provide that also for Sparse Linear Algebra even better even for CSR format, because there are plenty Sparse Matrix formats for Sparse linear Algebra.

Thanks.

[J-P. Rosen]

Le 24/01/2023 à 02:04, 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...

[Gautier write-only address]

There are plenty of legitimate complaints and good ideas in this thread.

- need for a package for large matrices (issue with stack allocation)
- should support approximations (at least, floating-point) real and complex numbers (and perhaps others)
- should support various matrix storages: sparse storages, band storages, ...

Fortunately Ada (not "ADA" = "American Dentist Association" and some others...) is powerful enough to have such components written in Ada itself (with GNAT, it is done exactly like this).
So _we_ don't need to wait that _they_ do anything about it :-)

So, if I summarize the ideas discussed and combine with stuff grabbed from a toolbox of mine ( https://github.com/zertovitch/mathpaqs
), I obtain the following specification.
The choice between different kind of matrix storages would be trivial.
Comments are welcome!

---8<--------8<--------8<--------8<--------8<--------8<-----
-- Draft of a specification for an universal matrix package.

-- The elements of the vectors and matrices can be of any algebraic ring
-- and their implementation, eventually approximate:
--
-- - integers (stored as Integer, Integer_n or Big_Integer)
-- - modular integers
-- - real numbers (approximated as floating-point, fixed-point or Big_Real)
-- - complex numbers
-- - rational numbers
-- - polynomials (of real, complex, rational, ...)
-- - other matrices
-- - ...

with Ada.Containers.Vectors;

generic

type Ring_Element is private; -- Element of an algebraic ring

zero, one : Ring_Element; -- 0 and 1 elements

with function "-" (a : Ring_Element) return Ring_Element; -- Unary operators

with function "+" (a, b : Ring_Element) return Ring_Element; -- Binary operators
with function "-" (a, b : Ring_Element) return Ring_Element;
with function "*" (a, b : Ring_Element) return Ring_Element;

package Universal_Matrices is

package UM_Vectors is new Ada.Containers.Vectors (Positive, Ring_Element);

subtype Vector is UM_Vectors.Vector;

-------------------------
-- Vector operations --
-------------------------

function "*" (v : Vector; factor : Ring_Element) return Vector;
function "*" (factor : Ring_Element; v : Vector) return Vector;

-- v := factor * v :
procedure Scale_Left (v : in out Vector; factor : Ring_Element);

function "-" (v : Vector) return Vector;

function "+" (v, w : Vector) return Vector;
-- v := v + w :
procedure Add (v : in out Vector; w : Vector);
-- v := v + factor * w :
procedure Add_Left_Scaled (v : in out Vector; factor : Ring_Element; w : Vector);

function "-" (v, w : Vector) return Vector;
-- v := v - w :
procedure Subtract (v : in out Vector; w : Vector);

function "*" (v, w : Vector) return Ring_Element;

-- Euclidean norm and distance:
function Square_L2_Norm (v : Vector) return Ring_Element;
function Square_L2_Distance (v, w : Vector) return Ring_Element;

-------------------------------------------------------------------
-- Root matrix type. --
-- Possible derivations: dense, sparse, band storage matrices. --
-------------------------------------------------------------------

type Matrix is interface;

-------------------------
-- Matrix operations --
-------------------------

function Transpose (A : Matrix) return Matrix is abstract;
function Identity (order : Positive) return Matrix is abstract;
function "*" (factor : Ring_Element; A : Matrix) return Matrix is abstract;
function "*" (A : Matrix; factor : Ring_Element) return Matrix is abstract;
function "*" (A, B : Matrix) return Matrix is abstract;
function "+" (A, B : Matrix) return Matrix is abstract;
function "-" (A, B : Matrix) return Matrix is abstract;

-- Matrix-Vector operations

function "*" (A : Matrix; x : Vector) return Vector is abstract;

end Universal_Matrices;

[Gautier write-only address]

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

[Jim Paloander]

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?

[Jim Paloander]

> > 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?

[J-P. Rosen]

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

[Gautier write-only address]

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

[Jim Paloander]

> ...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.

[Jerry]

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

[Jim Paloander]

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.

[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?
>

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