Is my version of gfortran just too old?

[James Van Buskirk]

I encountered the following bug:

D:\gfortran\james\cheby2>type bug1.f90
function insert(lower,mi)
implicit none
integer, intent(in) :: lower(:)
integer, intent(in) :: mi
integer(size(lower)+1) insert
integer loc

loc = maxloc(lower,mask = mi > lower,dim = 1)
insert = [lower(1:loc),mi,lower(loc+1:)]
end function insert

D:\gfortran\james\cheby2>gfortran -c bug1.f90
bug1.f90:5:22:

integer(size(lower)+1) insert
1
Error: Deferred array 'lower' at (1) is not permitted in an initialization
expre
ssion
bug1.f90:1:6:

function insert(lower,mi)
1
Error: Function 'insert' at (1) has no IMPLICIT type

gfortran -v says
gcc version 8.1.0 (x86_64-posix-seh-rev0, Built by MinGW-W64 project)

Is there a newer version that fixes this?

[Thomas Koenig]

James Van Buskirk <not_...@comcast.net> schrieb:

> I encountered the following bug:
>
> D:\gfortran\james\cheby2>type bug1.f90
> function insert(lower,mi)
> implicit none
> integer, intent(in) :: lower(:)
> integer, intent(in) :: mi
> integer(size(lower)+1) insert
> integer loc
>
> loc = maxloc(lower,mask = mi > lower,dim = 1)
> insert = [lower(1:loc),mi,lower(loc+1:)]
> end function insert
>
> D:\gfortran\james\cheby2>gfortran -c bug1.f90
> bug1.f90:5:22:
>
> integer(size(lower)+1) insert
> 1
> Error: Deferred array 'lower' at (1) is not permitted in an initialization
> expre
> ssion

Do you want to select the kind of the integer dynamically? That
is not possible, the KIND has to be a compile-time constant.

If you meant to write

integer :: insert(size(lower)+1)

that should work (without having tested it).

[James Van Buskirk]

"Thomas Koenig" wrote in message
news:tjp54r$2il75$1...@newsreader4.netcologne.de...

> integer :: insert(size(lower)+1)

Thanks. I was translating code from Matlab to Fortran and sort
of forgot the syntax rules. This was indeed my mistake.

[gah4]

And, reading along, I knew what it was supposed to do, and didn't notice.

Or, maybe, it is time for lunch.

I do hope you don't use this too much. If I understand it, it
allocates a whole new array each time. I am not sure what
happens to the old one.

[Thomas Koenig]

gah4 <ga...@u.washington.edu> schrieb:

> On Monday, October 31, 2022 at 12:40:13 PM UTC-7, James Van Buskirk wrote:
>> "Thomas Koenig" wrote in message
>> news:tjp54r$2il75$1...@newsreader4.netcologne.de...
>
>> > integer :: insert(size(lower)+1)
>
>> Thanks. I was translating code from Matlab to Fortran

A step in the right direction :-)

>and sort
>> of forgot the syntax rules. This was indeed my mistake.
>
> And, reading along, I knew what it was supposed to do, and didn't notice.
>
> Or, maybe, it is time for lunch.
>
> I do hope you don't use this too much. If I understand it, it
> allocates a whole new array each time. I am not sure what
> happens to the old one.

When it goes out of scope, it no longer uses memory. (I was
going to say it is deallocated, but that has a precise meaning
in Fortran, so I didn't).

There are ways to create memory leaks in Fortran, but they
all involve pointers.

[gah4]

On Monday, October 31, 2022 at 1:20:35 PM UTC-7, Thomas Koenig wrote:

(snip, I wrote)

> > I do hope you don't use this too much. If I understand it, it
> > allocates a whole new array each time. I am not sure what
> > happens to the old one.

> When it goes out of scope, it no longer uses memory. (I was
> going to say it is deallocated, but that has a precise meaning
> in Fortran, so I didn't).

> There are ways to create memory leaks in Fortran, but they
> all involve pointers.

It would help to show how the function is called.

And the calling routine might even use pointers.

In any case, there is at least one new array allocated for each call.
If you aren't careful, there might be two or three.

OK, say you do:

integer, allocatable :: y(:)

y = insert(x,3)

As well as I know it, insert allocates an array to return the result, a new
copy of y is allocated (allocate on assignment) and the return array goes away.

Even more interesting:

do i=1,10
y = insert(y, i)
end do

Now the original y gets deallocated, though two allocations each
time through the loop.

integer, pointer :: z(:)

z => insert(x,3)

I don't know if that is allowed., but I think

z => insert(z, 3)

does leak memory for z.

You could also have a function return an allocatable array, and I think:

integer, allocatable :: w(:)

! (allocate and use w)

move_alloc(insert(w,3), w)

So, only one allocation each time.

[James Van Buskirk]

"gah4" wrote in message
news:d604914b-d8a2-4857...@googlegroups.com...

> y = insert(x,3)

This was close to my usage, which actually was:

lower_new = insert(lower,mi)

It replaced the Matlab code:

lower_new = sort([mi lower])

Both methods are highly efficient in that they reliably insert
the new element mi into the sorted array lower, while still
occupying perhaps a couple of ppm of overall CPU time.
Of course, I admit to having gotten a little cross-eyed in
composing the Fortran replacement :)
Replacing code like

err = jac\f;

or

plot(dplot',xplot,'b-');

is much more problematic.

[gah4]

On Monday, October 31, 2022 at 5:45:20 PM UTC-7, James Van Buskirk wrote:
> "gah4" wrote in message
> news:d604914b-d8a2-4857...@googlegroups.com...
>
> > y = insert(x,3)
>
> This was close to my usage, which actually was:
>
> lower_new = insert(lower,mi)
>
> It replaced the Matlab code:
>
> lower_new = sort([mi lower])
>
> Both methods are highly efficient in that they reliably insert
> the new element mi into the sorted array lower, while still
> occupying perhaps a couple of ppm of overall CPU time.
> Of course, I admit to having gotten a little cross-eyed in
> composing the Fortran replacement :)
> Replacing code like

So I presume you are not doing this for millions of billions of elements.

Some years ago, I was working on a C program that reads a file into memory.
It had a loop that reads a line, and uses strcat() to add it to the buffer.

And this was done for maybe 1,000,000 elements. (From 60 character lines.)

That is O(n**2), as is I believe yours.

But if n isn't so big, it should be fine.

And as previously, as far as I know it allocates two arrays each
time through the loop.

[gah4]

(snip)

> And this was done for maybe 1,000,000 elements. (From 60 character lines.)

> That is O(n**2), as is I believe yours.

OK, I am running this one:

integer, allocatable :: x(:)
interface
function insert(y, mi)
integer, intent(in) :: y(:), mi
integer insert(size(y)+1)
end function
end interface

allocate(x(1))
x(1)=1
do i=1000000,1,-1
x=insert(x,i)
enddo
print *,x(1), x(size(x)), size(x)
end

function insert(lower,mi)
implicit none
integer, intent(in) :: lower(:)
integer, intent(in) :: mi

integer insert(size(lower)+1)

integer loc

loc = maxloc(lower,mask = mi > lower,dim = 1)
insert = [lower(1:loc),mi,lower(loc+1:)]
end function insert

Any suggestions on how long it should take?

I will post when it is done.

[gah4]

On Tuesday, November 1, 2022 at 11:10:12 PM UTC-7, gah4 wrote:
> (snip)

> do i=1000000,1,-1
> x=insert(x,i)
> enddo

(snip)

> Any suggestions on how long it should take?

> I will post when it is done.

And the answer is 82 minutes.

[Ron Shepard]

You should consider a different underlying data structure than a linear
array. I suggest looking at a binary search tree (BST). You can code a
simple one in fortran in about the same number of lines as your posted
example.

Your insertion step requires O(n) effort for each insertion where n is
the current number of data points. A BST insertion requires O(log(n))
effort for each insertion. There are some other data structures that
might be useful too, but this is simple and should get you started in
the right direction.

$.02 -Ron Shepard

[gah4]

On Wednesday, November 2, 2022 at 7:44:12 AM UTC-7, Ron Shepard wrote:

(snip, I wrote)

> > And the answer is 82 minutes.

> You should consider a different underlying data structure than a linear
> array. I suggest looking at a binary search tree (BST). You can code a
> simple one in fortran in about the same number of lines as your posted
> example.

> Your insertion step requires O(n) effort for each insertion where n is
> the current number of data points. A BST insertion requires O(log(n))
> effort for each insertion. There are some other data structures that
> might be useful too, but this is simple and should get you started in
> the right direction.

Just to be sure, it is James data structure and not mine.
I was trying to see how to call insert, and this was my first time
with an INTERFACE for a function returning a variable number of
array elements.

There is a popular binary search tree using pointers, but there is
also the heap, and especially binary heap:

https://en.wikipedia.org/wiki/Heap_(data_structure)

The binary heap is a binary tree without the pointers.

But that only helps if you don't reallocate it (once, or twice,
or more) each time.

But if N is 10, or 100, or even 1000, then it isn't so bad.

My calculation says for N=1,000,000,000 it is 156 years.

[James Van Buskirk]

"gah4" wrote in message
news:b4284a2d-ba24-48fb...@googlegroups.com...

> On Wednesday, November 2, 2022 at 7:44:12 AM UTC-7, Ron Shepard wrote:

> (snip, I wrote)

> > > And the answer is 82 minutes.

> > You should consider a different underlying data structure than a linear
> > array. I suggest looking at a binary search tree (BST). You can code a
> > simple one in fortran in about the same number of lines as your posted
> > example.

> > Your insertion step requires O(n) effort for each insertion where n is
> > the current number of data points. A BST insertion requires O(log(n))
> > effort for each insertion. There are some other data structures that
> > might be useful too, but this is simple and should get you started in
> > the right direction.

> Just to be sure, it is James data structure and not mine.
> I was trying to see how to call insert, and this was my first time
> with an INTERFACE for a function returning a variable number of
> array elements.

> There is a popular binary search tree using pointers, but there is
> also the heap, and especially binary heap:

For this particular data structure there was no impact on performance,
but there was another one where the results were accumulated in
arrays. That really did slow things down tremendously but I hadn't
noticed before because the rest of my program was too slow for
the arrays to get that large. Since the results were to be plotted,
their order of generation was what mattered so I dropped them into
a linked list and the longest-running cases sped up by an order of
magnitude.

In the original example the sorted arrays had to be swept through
in order several times after each update so only arrays made
sense.

[gah4]

On Thursday, November 10, 2022 at 10:39:17 AM UTC-8, James Van Buskirk wrote:

(snip, I wrote)

> > There is a popular binary search tree using pointers, but there is
> > also the heap, and especially binary heap:

> For this particular data structure there was no impact on performance,
> but there was another one where the results were accumulated in
> arrays. That really did slow things down tremendously but I hadn't
> noticed before because the rest of my program was too slow for
> the arrays to get that large. Since the results were to be plotted,
> their order of generation was what mattered so I dropped them into
> a linked list and the longest-running cases sped up by an order of
> magnitude.

> In the original example the sorted arrays had to be swept through
> in order several times after each update so only arrays made
> sense.

The binary heap is an interesting data structure, and one well
designed for languages were arrays start at 1. (It is usual to
ignore the 0th element in C.)

If you consider a balanced binary tree, and read off the elements
row by row, and then store them in an array:

the root goes in element 1.
The next two in elements 10 and 11. (Binary addressing.)
The next four in elements 100, 101, 110, 111.

Because of the way the addressing works, you can go up the
tree, just shifting the address right one bit at a time.

There is a convenient algorithm to go through the elements
in sorted order, based on binary arithmetic. Maybe not quite
as fast as a sorted array, but pretty fast.

And it is O(log N) to add elements to the heap.

[James Van Buskirk]

"gah4" wrote in message
news:1868633b-01d9-4c22...@googlegroups.com...

> The binary heap is an interesting data structure, and one well
> designed for languages were arrays start at 1. (It is usual to
> ignore the 0th element in C.)

> If you consider a balanced binary tree, and read off the elements
> row by row, and then store them in an array:

> the root goes in element 1.
> The next two in elements 10 and 11. (Binary addressing.)
> The next four in elements 100, 101, 110, 111.

> Because of the way the addressing works, you can go up the
> tree, just shifting the address right one bit at a time.

> There is a convenient algorithm to go through the elements
> in sorted order, based on binary arithmetic. Maybe not quite
> as fast as a sorted array, but pretty fast.

> And it is O(log N) to add elements to the heap.

I think that I have previously post AVL tree code to this forum.

It is unfortunate that you have not fleshed your ideas out by
posting code like:

integer n
real(REAL64) harvest
.
.
.
call RANDOM_NUMBER(harvest)
n = HUGE(n)*harvest
! Now insert n into your data structure...

So as to compare performance between the example you posted
and your proposed data structure solution...

[gah4]

On Thursday, November 10, 2022 at 3:19:35 PM UTC-8, James Van Buskirk wrote:

(snip)

> I think that I have previously post AVL tree code to this forum.

> It is unfortunate that you have not fleshed your ideas out by
> posting code like:

> integer n
> real(REAL64) harvest

> call RANDOM_NUMBER(harvest)
> n = HUGE(n)*harvest
> ! Now insert n into your data structure...

> So as to compare performance between the example you posted
> and your proposed data structure solution...

Many years ago for a CS class assignment I did a heapsort program.
(In PDP-10 ALGOL!) That was based on the explanation in
Knuth's volume 3, which I think I bought to do that assignment.

I have not actually tried to write one since then.

I think, though, it depends somewhat on how, and how often, you
go through the sorted data.


As well as I know, your original program goes through the array
four times for each call. One to create the mask. One to extract
the location from the mask. One to copy into the return array.
One to copy the return array to the called program array.

Though that only means a 4 in front of the N**2.

OK, now another Fortran tree story.

The IBM Fortran H compiler stores the symbol table in six
balanced trees, instead of a hash table. In the manual, they
suggest that for faster compilation, you distribute your names
equally between 1 and 6 characters!

No suggestion that you choose names for readability, though.