[stdlib_linalg] Update eye function.

[zoziha]

• Update eye function, add support that matrix is not square.
• Based on storage considerations, integer(int8) will used as the return value type.
• allocation or automatic: allocation because of optional dim2!(see [stdlib_linalg] Update eye function. #481 (comment))

Description

That eye creates a non-square matrix would be a nice update.

Notes

More routines to do, not this PR
see #476.

[Jim-215-Fisher]

I will review this one.

[awvwgk]

Thanks for sharing, @zoziha. Thanks for reviewing @Jim-215-Fisher. Looks good to me, except for one issue with merge.

[nshaffer]

Made a few suggestions, but mostly seems good to me. Let me know if you have any questions about my review.

[jvdp1]

Thanks for this PR. I have a few suggestions/comments/questions.

[jvdp1]

Could you avoid the allocatable status, like for other functions in stdlib?

[zoziha]

(see #478 (comment), #477 (comment), and https://community.intel.com/t5/Intel-Fortran-Compiler/how-to-set-the-stack-size/td-p/933530)
It seems that ifort places the automatic array in the heap by default, which leads to a stack overflow when the dimension of the automatic array is large. Or is there any disadvantage of the allocatable array, need to avoid it?

[jvdp1]

My comment was based on @certik and @milancurcic 's comments (this comment and following ones).

Furthermore, the rules for stack vs heap can be different between compilers/OS, and can be also changed within a compiler through a compiler option (e.g., using Gfortran with -Ofast will imply -fstack-arrays, that will put all arrays of unknown size and temporary arrays onto stack memory, which can create issues for very large arrays).

Therefore, I am in favor to always use automatic arrays when possible.

@certik @milancurcic @awvwgk @ivan-pi Any comments on that topic?

[zoziha]

This may indeed be a big controversy, it seems to be: @awvwgk agrees with allocation array, @jvdp1 agrees with automatic array.
By default, ifort's automatic array solution faces large-scale arrays, it is easier to report an error: stack overflow. (see #478 (comment), #477 (comment), and https://community.intel.com/t5/Intel-Fortran-Compiler/how-to-set-the-stack-size/td-p/933530)

[milancurcic]

Yes, I suggest using automatic array wherever possible, and allocatable otherwise. In this case, it seems like you need to make it allocatable because the second dimension depends on the presence of the optional dummy argument.

However, if you re-write this to be two specific procedures under one generic name, instead of using the optional argument, then I think you don't need allocatable for the result.

[awvwgk]

Yes, I'm not a fan of automatic arrays, since usage with ifort will lead to stack overflows. The option to allocate automatic arrays on the heap rather than on the stack is there, for ifort it happens to be buggy (heap allocation is not always freed). gfortran is doing a little better here, but it is an issue as well.

So far automatic array usage has resulted in a significant amount of support request from users for me. Setting the system stack with ulimit -s unlimited works, still this is something the user has to remember.

[milancurcic]

@awvwgk I don't think that's a good reason to avoid a language feature. If a feature (and a very old one) doesn't work with a compiler, let's report it. Avoiding it only helps permeate the bug with the compiler. I've been using this with ifort since 2010 and haven't seen an issue. It's possible that I didn't notice it.

The downside to the allocatable approach is that it's always on the heap. The user doesn't have the choice and always gets the less-performant implementation. Also, in the optional argument approach, it's possible (but I'm not sure) that the compiler would evaluate the branch inside the function at run-time rather than compile-time. If at run-time, that would be another performance penalty. It's possible that performance is less important for eye, I don't know. And the advantage to the allocatable approach is, of course, that it's always on the heap. :)

@zoziha When resolving conversations as a GitHub feature, I suggest to leave a comment about what you decided, otherwise discussions get hidden but aren't really resolved. Maybe we should document this in the Code Review guide.

[jvdp1]

I agree with @milancurcic: issues with some compilers/options/OS/... are not good reasons to avoid a feature of the language.
Futhermore, because it is a small function that would be most often (I guess) used in expressions (e.g., arr1 = eye(5) * scalar), I am also in favor to apply @milancurcic suggestion (i.e., to have 2 specific procedures under one generic name)

[milancurcic]

@jvdp1 okay if we do it in a separate PR and merge this one first? We could evaluate the performance difference too.

[zoziha]

@milancurcic Ok, got it.

1. I changed the return value type of the eye function to integer(int8);
2. allocatable is used because of the existence of optional dim2 argument.

At present, the questions raised by @jvdp1 , I think they are answered and resolved.

[milancurcic]

Okay, thanks. If needed we can revisit the allocatable result later. If there are no objections let's merge this tomorrow.

[jvdp1]

Ok for me. Le dim. 22 août 2021 à 17:22, Milan Curcic ***@***.***> a écrit :

…

***@***.**** commented on this pull request. ------------------------------ In src/stdlib_linalg.fypp <#481 (comment)>: > - integer(int8) :: res(n, n) - integer :: i - res = 0 - do i = 1, n - res(i, i) = 1 - end do - end function eye + !> Version: experimental + !> + !> Constructs the identity matrix. + !> ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix)) + pure function eye(dim1, dim2) result(result) + + integer, intent(in) :: dim1 + integer, intent(in), optional :: dim2 + integer, allocatable :: result(:, :) @jvdp1 <https://github.com/jvdp1> okay if we do it in a separate PR and merge this one first? We could evaluate the performance difference too. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#481 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AD5RO7BKDXEGU3C6QX4DY3DT6EI3HANCNFSM5BLFGHUA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email> .