Addition of diag, eye, and trace #170

ivan-pi · 2020-04-10T11:52:23Z

See proposal in #169 . This is a draft pull request.

This is my first time using the fypp preprocessor, so I decided to only focus on rank-2 arrays to keep things simple and submit a pull request before my focus shifts away from stdlib. Hopefully the API will not need to change if we decide to generalize these to higher rank-matrices, otherwise the generalizations could be placed into a stdlib_tensor module in the future.

I am still missing the tests, but thought it would be beneficial to get some feedback on the functions and documentation first.

@jvdp1 Did you use fypp to generate the tests for the statistical routines and then committed only the generated .f90 files?

certik

This looks great!

jvdp1

Thanks @ivan-pi for these functions. It looks nice.
Some comments:

It would maybe good to move the functions in submodules. It avoids the compilation of alll functions when modifying only some of them.
I think it could be easily extended to ranks >2 with the same API quite easily using fypp. However, EDIT: I am NOT the target user for such functions.
I wrote the tests by hand.

src/stdlib_experimental_linalg.fypp

src/stdlib_experimental_linalg.md

jvdp1 · 2020-04-10T17:42:18Z

src/tests/linalg/test_linalg.f90

+
+  print *, eye(5)
+
+  print *, all(eye(4) == diag([1,1,1,1]))


Could the check function be used here?

I will add more tests with calls to check. This was just temporary for development purposes.

Nice. For internal stdlib tests, I suggest providing custom message in each check call so we can identify which failed, and also with warn=.true. to allow multiple tests to fail in a single test program.

and also with warn=.true. to allow multiple tests to fail in a single test program.

@milancurcic With warn = .true. multiple tests can fail, but ctest will report that all tests passed. It is ok when we develop it. But before merging into master, we should probably have warn=.false..

You're right, I forgot about that. We shouldn't do it then. It's important that ctest correctly reports failure.

ivan-pi · 2020-04-10T18:20:30Z

It would maybe good to move the functions in submodules. It avoids the compilation of alll functions when modifying only some of them.

I agree. I expect that in the future we will likely have several linear algebra submodules (array creation routines, solvers, matrix properties...).

I think it could be easily extended to ranks >2 with the same API quite easily using fypp. However, I am the target user for such functions.

I don't think a tensor version of eye is trivial, at least not under a common interface like eye(n [,rank]). We could have multiple versions for different ranks, e.g. eye (equal to eye2), eye3, eye4, eye5, etc.?

I can try to make trace rank-generic for the main diagonal, and also have a version similar to numpy

numpy.trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None)

which calculates the trace of each two-dimensional slice returning an array. But since the name trace is typically linked with linear algebra, for the tensor "trace" it might be a better idea to have a general tensor_contraction function (https://en.wikipedia.org/wiki/Tensor_contraction).

With diag comes again the question, how to specify the rank of the output? Also the diagonal specifier is not immediately clear. Is the aim to operate on the individual 2d slices, or somehow set/extract a diagonal/slice in >=3 dimensions?
(image taken from here)

I wrote the tests by hand.

I applaud your effort. 👏

jvdp1 · 2020-04-10T18:32:22Z

Is the aim to operate on the individual 2d slices, or somehow set/extract a diagonal/slice in >=3 dimensions?

Good point. I didn't think so far.
Anyway, this is a nice proposal that could help many people as is.

milancurcic

Looks great, IMO good to go pending any requested changes.

jvdp1

Pending the tests, it looks good to me. Thank you.

ivan-pi · 2020-04-10T22:20:56Z

In the routine to extract the k-th diagonal from a matrix, I have the following interface:

     function diag_${t1[0]}$${k1}$_mat_k(A,k) result(res)
        ${t1}$, intent(in) :: A(:,:)
        integer, intent(in) :: k
        ${t1}$ :: res(minval(shape(A))-abs(k))
    end function

Interestingly, if abs(k) > minval(shape(A)), the size of res is negative, however this gets correctly translated into a zero-size array. I've double checked this with the following program:

implicit none
real :: res(-5)
print *, size(res)
end

which prints "0" for both the ifort and gfortan compilers without any warnings. So it looks like Fortran had my back covered 😎 .

Edit: a related issue appeared in the Intel Fortran Forum: https://software.intel.com/en-us/forums/intel-fortran-compiler/topic/297666. Here is a quote from the Fortran standard pertaining to this issue in the context of allocatable arrays:

If the lower bound is omitted, the default value is 1. If the upper bound is less than the lower bound, the extent in that dimension is zero and the array has zero size.

Explicit-shape arrays also follow the same behavior. (e.g. https://groups.google.com/forum/#!topic/comp.lang.fortran/xGwU2diHm60)

ivan-pi · 2020-04-11T00:57:33Z

Is there any interest to extend diag also to logical arrays?

jvdp1 · 2020-04-11T20:23:41Z

Is there any interest to extend diag also to logical arrays?

Not from my side. For what would you use such a matrix?

ivan-pi · 2020-04-13T14:21:12Z

It could be used as a mask to unpack values into a certain diagonal or more generally as a mask in other functions. The numpy diag also allows it, in MATLAB I think logical masks are binary arrays, so it probably also just works.

jvdp1 · 2020-04-13T17:12:14Z

It could be used as a mask to unpack values into a certain diagonal or more generally as a mask in other functions. The numpy diag also allows it, in MATLAB I think logical masks are binary arrays, so it probably also just works.

If other languages support it, it is probably good to support it too. Up to you to decide. I don't think I would need such a logical function.

ivan-pi · 2020-05-03T20:05:02Z

I was planning to have a go at the logical kind version, but would like to clarify first if it would suffice to support only the default logical kind or all those supported by the compiler (for some odd reason?):

ipribec@ipribec-T530:~/fortran$ cat test_logical_kinds.f90 
use iso_fortran_env
print *, logical_kinds
end
ipribec@ipribec-T530:~/fortran$ gfortran test_logical_kinds.f90 && ./a.out
           1           2           4           8          16
ipribec@ipribec-T530:~/fortran$ ifort test_logical_kinds.f90 && ./a.out
           1           2           4           8
ipribec@ipribec-T530:~/fortran$

Any thoughts?

If not, pending a review of the tests, I would merge this.

certik · 2020-05-03T20:12:41Z

I think we can merge this now. We are in experimental, so any such changes / modifications can always be done later.

jvdp1

LGTM. I only have a couple of minor comments.

jvdp1 · 2020-05-05T19:17:08Z

src/tests/linalg/test_linalg.f90

+          msg="all(diag(a) == pack(a,mask))", warn=warn)
+    call check(all(diag(diag(a)) == merge(a,0_int16,mask)), &
+          msg="all(diag(diag(a)) == merge(a,0_int16,mask)) failed.", warn=warn)
+    a = unpack(int([1,2,3,4],int16),eye(n)==1,a)


I believe that this line could be removed.

jvdp1 · 2020-05-05T19:22:34Z

src/tests/linalg/test_linalg.f90

+    a = reshape([(i**2,i=1,n*(n-1))],[n,n-1])
+    ans = sum([1._dp,36._dp,121._dp])
+
+    call check(abs(trace(a) - ans) < epsilon(1.0_dp), &


When submitting tests for mean and var, it was mentioned that epsilon() is too strict. Therefore, you should multiply it by e.g., 1000 (see here for the current implementation) to avoid any issues.

Co-Authored-By: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

ivan-pi · 2020-05-05T21:00:53Z

Thanks @jvdp1, @milancurcic, and @certik for the review and approvals. I've squashed the commits and will merge after the tests are completed.

ivan-pi requested review from milancurcic and jvdp1 April 10, 2020 11:52

ivan-pi marked this pull request as draft April 10, 2020 11:52

ivan-pi mentioned this pull request Apr 10, 2020

diag, eye, and trace #169

Closed

certik approved these changes Apr 10, 2020

View reviewed changes

jvdp1 reviewed Apr 10, 2020

View reviewed changes

milancurcic approved these changes Apr 10, 2020

View reviewed changes

jvdp1 approved these changes Apr 10, 2020

View reviewed changes

ivan-pi marked this pull request as ready for review May 3, 2020 20:05

jvdp1 reviewed May 5, 2020

View reviewed changes

Implementation and tests for diag, eye, and trace.

d01f657

Co-Authored-By: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

ivan-pi force-pushed the ivan-pi branch from ca439e6 to d01f657 Compare May 5, 2020 20:58

ivan-pi merged commit bbecb52 into fortran-lang:master May 5, 2020

ivan-pi deleted the ivan-pi branch May 5, 2020 21:07

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Addition of diag, eye, and trace #170

Addition of diag, eye, and trace #170

ivan-pi commented Apr 10, 2020

certik left a comment

jvdp1 left a comment •

edited

Loading

jvdp1 Apr 10, 2020

ivan-pi Apr 10, 2020

milancurcic Apr 10, 2020

jvdp1 Apr 10, 2020

milancurcic Apr 10, 2020

ivan-pi commented Apr 10, 2020 •

edited

Loading

jvdp1 commented Apr 10, 2020

milancurcic left a comment

jvdp1 left a comment

ivan-pi commented Apr 10, 2020 •

edited

Loading

ivan-pi commented Apr 11, 2020

jvdp1 commented Apr 11, 2020

ivan-pi commented Apr 13, 2020

jvdp1 commented Apr 13, 2020

ivan-pi commented May 3, 2020

certik commented May 3, 2020

jvdp1 left a comment

jvdp1 May 5, 2020

jvdp1 May 5, 2020

ivan-pi commented May 5, 2020

Addition of diag, eye, and trace #170

Addition of diag, eye, and trace #170

Conversation

ivan-pi commented Apr 10, 2020

certik left a comment

Choose a reason for hiding this comment

jvdp1 left a comment • edited Loading

Choose a reason for hiding this comment

jvdp1 Apr 10, 2020

Choose a reason for hiding this comment

ivan-pi Apr 10, 2020

Choose a reason for hiding this comment

milancurcic Apr 10, 2020

Choose a reason for hiding this comment

jvdp1 Apr 10, 2020

Choose a reason for hiding this comment

milancurcic Apr 10, 2020

Choose a reason for hiding this comment

ivan-pi commented Apr 10, 2020 • edited Loading

jvdp1 commented Apr 10, 2020

milancurcic left a comment

Choose a reason for hiding this comment

jvdp1 left a comment

Choose a reason for hiding this comment

ivan-pi commented Apr 10, 2020 • edited Loading

ivan-pi commented Apr 11, 2020

jvdp1 commented Apr 11, 2020

ivan-pi commented Apr 13, 2020

jvdp1 commented Apr 13, 2020

ivan-pi commented May 3, 2020

certik commented May 3, 2020

jvdp1 left a comment

Choose a reason for hiding this comment

jvdp1 May 5, 2020

Choose a reason for hiding this comment

jvdp1 May 5, 2020

Choose a reason for hiding this comment

ivan-pi commented May 5, 2020

jvdp1 left a comment •

edited

Loading

ivan-pi commented Apr 10, 2020 •

edited

Loading

ivan-pi commented Apr 10, 2020 •

edited

Loading