Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[stdlib_linalg] Update eye function. #481

Merged
merged 6 commits into from
Aug 23, 2021
Merged

[stdlib_linalg] Update eye function. #481

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
ebc61cd
Update eye function.
zoziha Aug 1, 2021
8745725
Improve `eye` function help docs.
zoziha Aug 3, 2021
c2f0abf
Update `eye` func:
zoziha Aug 5, 2021
c20f69f
Improve `eye` func test and doc.
zoziha Aug 6, 2021
c326353
Update `eye` func: `integer(int8)`
zoziha Aug 21, 2021
2e4d681
Merge branch 'master' into update_eye
zoziha Aug 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 32 additions & 5 deletions doc/specs/stdlib_linalg.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,30 +101,57 @@ end program demo_diag5

Experimental

### Class

Pure function.

### Description

Construct the identity matrix
Construct the identity matrix.

### Syntax

`I = [[stdlib_linalg(module):eye(function)]](n)`
`I = [[stdlib_linalg(module):eye(function)]](dim1 [, dim2])`

### Arguments

`n`: Shall be a scalar of default type `integer`.
`dim1`: Shall be a scalar of default type `integer`.
This is an `intent(in)` argument.

`dim2`: Shall be a scalar of default type `integer`.
This is an `intent(in)` and `optional` argument.

### Return value

Returns the identity matrix, i.e. a square matrix with ones on the main diagonal and zeros elsewhere. The return value is of type `integer(int8)`.
Return the identity matrix, i.e. a matrix with ones on the main diagonal and zeros elsewhere. The return value is of type `integer`.
zoziha marked this conversation as resolved.
Show resolved Hide resolved

#### Warning

Since the result of `eye` is of `integer` type, one should be careful about using it in arithmetic expressions. For example:
```fortran
real :: A(:,:)
!> Be careful
A = eye(2,2)/2 !! A == 0.0
!> Recommend
A = eye(2,2)/2.0 !! A == diag([0.5, 0.5])
```

### Example

```fortran
program demo_eye1
use stdlib_linalg, only: eye
implicit none
integer :: i(2,2)
real :: a(3,3)
A = eye(3)
real :: b(2,3) !! Matrix is non-square.
complex :: c(2,2)
I = eye(2) !! [1,0; 0,1]
A = eye(3) !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
A = eye(3,3) !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
B = eye(2,3) !! [1.0,0.0,0.0; 0.0,1.0,0.0]
C = eye(2,2) !! [(1.0,0.0),(0.0,0.0); (0.0,0.0),(1.0,0.0)]
C = (1.0,1.0)*eye(2,2) !! [(1.0,1.0),(0.0,0.0); (0.0,0.0),(1.0,1.0)]
end program demo_eye1
```

Expand Down
3 changes: 2 additions & 1 deletion src/Makefile.manual
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,8 @@ stdlib_string_type.o: stdlib_ascii.o \
stdlib_strings.o: stdlib_ascii.o \
stdlib_string_type.o \
stdlib_optval.o
stdlib_math.o: stdlib_kinds.o
stdlib_math.o: stdlib_kinds.o \
stdlib_optval.o
stdlib_math_linspace.o: \
stdlib_math.o
stdlib_math_logspace.o: \
Expand Down
38 changes: 24 additions & 14 deletions src/stdlib_linalg.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ module stdlib_linalg
!! ([Specification](../page/specs/stdlib_linalg.html))
use stdlib_kinds, only: sp, dp, qp, &
int8, int16, int32, int64
use stdlib_optval, only: optval
implicit none
private

Expand Down Expand Up @@ -82,20 +83,28 @@ module stdlib_linalg

contains

function eye(n) result(res)
!! version: experimental
!!
!! Constructs the identity matrix
!! ([Specification](../page/specs/stdlib_linalg.html#description_1))
integer, intent(in) :: n
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(:, :)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

@zoziha zoziha Aug 18, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

@milancurcic milancurcic Aug 21, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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


integer :: dim2_
integer :: i

dim2_ = optval(dim2, dim1)
allocate(result(dim1, dim2_))

result = 0
do i = 1, min(dim1, dim2_)
result(i, i) = 1
end do

end function eye

#:for k1, t1 in RCI_KINDS_TYPES
function trace_${t1[0]}$${k1}$(A) result(res)
Expand All @@ -108,4 +117,5 @@ contains
end do
end function trace_${t1[0]}$${k1}$
#:endfor
end module

end module stdlib_linalg
1 change: 1 addition & 0 deletions src/tests/Makefile.manual
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,3 +11,4 @@ all test clean:
$(MAKE) -f Makefile.manual --directory=stats $@
$(MAKE) -f Makefile.manual --directory=string $@
$(MAKE) -f Makefile.manual --directory=math $@
$(MAKE) -f Makefile.manual --directory=linalg $@
4 changes: 4 additions & 0 deletions src/tests/linalg/Makefile.manual
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
PROGS_SRC = test_linalg.f90


include ../Makefile.manual.test.mk
11 changes: 9 additions & 2 deletions src/tests/linalg/test_linalg.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,12 +81,19 @@ subroutine test_eye
integer :: i
write(*,*) "test_eye"

call check(all(eye(3,3) == diag([(1,i=1,3)])), &
msg="all(eye(3,3) == diag([(1,i=1,3)])) failed.",warn=warn)

rye = eye(3,4)
call check(sum(abs(rye(:,1:3) - diag([(1.0_sp,i=1,3)]))) < sptol, &
msg="sum(abs(rye(:,1:3) - diag([(1.0_sp,i=1,3)]))) < sptol failed", warn=warn)

call check(all(eye(5) == diag([(1,i=1,5)])), &
msg="all(eye(5) == diag([(1,i=1,5)] failed.",warn=warn)

rye = eye(6)
call check(sum(rye - diag([(1.0_sp,i=1,6)])) < sptol, &
msg="sum(rye - diag([(1.0_sp,i=1,6)])) < sptol failed.",warn=warn)
call check(sum(abs(rye - diag([(1.0_sp,i=1,6)]))) < sptol, &
msg="sum(abs(rye - diag([(1.0_sp,i=1,6)]))) < sptol failed.",warn=warn)

cye = eye(7)
call check(abs(trace(cye) - cmplx(7.0_sp,0.0_sp,kind=sp)) < sptol, &
Expand Down