From 54c9ca120db7697217d89e6e1484078db5c7d82b Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 2 Dec 2024 05:42:36 -0600 Subject: [PATCH 1/5] `eye`: introduce optional `mold` to set return type --- src/stdlib_linalg.fypp | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ad36ad677..b552ac97c 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -188,6 +188,16 @@ module stdlib_linalg #:endfor end interface + ! Identity matrix + interface eye + !! version: experimental + !! + !! Constructs the identity matrix + !! ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix)) + #:for k1, t1 in RCI_KINDS_TYPES + module procedure eye_${t1[0]}$${k1}$ + #:endfor + end interface eye ! Outer product (of two vectors) interface outer_product @@ -1327,11 +1337,13 @@ contains !> !> Constructs the identity matrix. !> ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix)) - pure function eye(dim1, dim2) result(result) + #:for k1, t1 in RCI_KINDS_TYPES + pure function eye_${t1[0]}$${k1}$(dim1, dim2, mold) result(result) integer, intent(in) :: dim1 integer, intent(in), optional :: dim2 - integer(int8), allocatable :: result(:, :) + ${t1}$, intent(in) #{if not 'int8' in t1}#, optional #{endif}#:: mold + ${t1}$, allocatable :: result(:, :) integer :: dim2_ integer :: i @@ -1339,12 +1351,13 @@ contains dim2_ = optval(dim2, dim1) allocate(result(dim1, dim2_)) - result = 0_int8 + result = 0 do i = 1, min(dim1, dim2_) - result(i, i) = 1_int8 + result(i, i) = 1 end do - end function eye + end function eye_${t1[0]}$${k1}$ + #:endfor #:for k1, t1 in RCI_KINDS_TYPES function trace_${t1[0]}$${k1}$(A) result(res) From f0260974ca7b58072b4b2db75b1d456b072b916d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 2 Dec 2024 05:55:33 -0600 Subject: [PATCH 2/5] `eye`: update documentation --- doc/specs/stdlib_linalg.md | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 2f03f0fad..43a89b4bc 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -239,33 +239,34 @@ Pure function. ### Description -Construct the identity matrix. +Constructs the identity matrix. ### Syntax -`I = ` [[stdlib_linalg(module):eye(function)]] `(dim1 [, dim2])` +`I = ` [[stdlib_linalg(module):eye(function)]] `(dim1 [, dim2] [, mold])` ### Arguments -`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. +- `dim1`: A scalar of type `integer`. This is an `intent(in)` argument and specifies the number of rows. +- `dim2`: A scalar of type `integer`. This is an optional `intent(in)` argument specifying the number of columns. If not provided, the matrix is square (`dim1 = dim2`). +- `mold`: A scalar of any supported `integer`, `real`, or `complex` type. This is an optional `intent(in)` argument. If provided, the returned identity matrix will have the same type and kind as `mold`. If not provided, the matrix will be of type `integer(int8)` by default. ### Return value -Return the identity matrix, i.e. a matrix with ones on the main diagonal and zeros elsewhere. The return value is of type `integer(int8)`. -The use of `int8` was suggested to save storage. +Returns the identity matrix, with ones on the main diagonal and zeros elsewhere. + +- By default, the return value is of type `integer(int8)`, which is recommended for storage efficiency. +- If the `mold` argument is provided, the return value will match the type and kind of `mold`, allowing for arbitrary `integer`, `real`, or `complex` return types. #### Warning -Since the result of `eye` is of `integer(int8)` type, one should be careful about using it in arithmetic expressions. For example: +When using the default `integer(int8)` type, be cautious when performing arithmetic operations, as integer division may occur. For example: + ```fortran -!> Be careful -A = eye(2,2)/2 !! A == 0.0 -!> Recommend -A = eye(2,2)/2.0 !! A == diag([0.5, 0.5]) +!> Caution: default type is `integer` +A = eye(2,2)/2 !! A == 0.0 due to integer division +!> Recommend using a non-integer type for division +A = eye(2,2, mold=1.0)/2 !! A == diag([0.5, 0.5]) ``` ### Example From 758ab3a1d3362d93162d1511198533dbdffcaaaa Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 2 Dec 2024 05:57:30 -0600 Subject: [PATCH 3/5] Update example_eye1.f90 --- example/linalg/example_eye1.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/example/linalg/example_eye1.f90 b/example/linalg/example_eye1.f90 index 7edda482b..c8ded5683 100644 --- a/example/linalg/example_eye1.f90 +++ b/example/linalg/example_eye1.f90 @@ -5,10 +5,10 @@ program example_eye1 real :: a(3, 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)] + 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 example_eye1 From 406636510c378dcb0f53b8ae3e9fedad34c60919 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 2 Dec 2024 07:04:58 -0600 Subject: [PATCH 4/5] fix `optional` interface --- src/stdlib_linalg.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index b552ac97c..2f14e512d 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1342,7 +1342,7 @@ contains integer, intent(in) :: dim1 integer, intent(in), optional :: dim2 - ${t1}$, intent(in) #{if not 'int8' in t1}#, optional #{endif}#:: mold + ${t1}$, intent(in) #{if 'int8' in t1}#, optional #{endif}#:: mold ${t1}$, allocatable :: result(:, :) integer :: dim2_ From 57e3f89d367174c908b78bb8dd744e9fc336de12 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 11 Dec 2024 09:14:28 +0100 Subject: [PATCH 5/5] make `eye` default type: `real(dp)` --- doc/specs/stdlib_linalg.md | 18 +++++++----------- src/stdlib_linalg.fypp | 2 +- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 43a89b4bc..fe5ac698c 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -249,28 +249,24 @@ Constructs the identity matrix. - `dim1`: A scalar of type `integer`. This is an `intent(in)` argument and specifies the number of rows. - `dim2`: A scalar of type `integer`. This is an optional `intent(in)` argument specifying the number of columns. If not provided, the matrix is square (`dim1 = dim2`). -- `mold`: A scalar of any supported `integer`, `real`, or `complex` type. This is an optional `intent(in)` argument. If provided, the returned identity matrix will have the same type and kind as `mold`. If not provided, the matrix will be of type `integer(int8)` by default. +- `mold`: A scalar of any supported `integer`, `real`, or `complex` type. This is an optional `intent(in)` argument. If provided, the returned identity matrix will have the same type and kind as `mold`. If not provided, the matrix will be of type `real(real64)` by default. ### Return value Returns the identity matrix, with ones on the main diagonal and zeros elsewhere. -- By default, the return value is of type `integer(int8)`, which is recommended for storage efficiency. +- By default, the return value is of type `real(real64)`, which is recommended for arithmetic safety. - If the `mold` argument is provided, the return value will match the type and kind of `mold`, allowing for arbitrary `integer`, `real`, or `complex` return types. -#### Warning - -When using the default `integer(int8)` type, be cautious when performing arithmetic operations, as integer division may occur. For example: +### Example ```fortran -!> Caution: default type is `integer` -A = eye(2,2)/2 !! A == 0.0 due to integer division -!> Recommend using a non-integer type for division -A = eye(2,2, mold=1.0)/2 !! A == diag([0.5, 0.5]) +!> Return default type (real64) +A = eye(2,2)/2 !! A == diag([0.5_dp, 0.5_dp]) +!> Return 32-bit complex +A = eye(2,2, mold=(0.0,0.0))/2 !! A == diag([(0.5,0.5), (0.5,0.5)]) ``` -### Example - ```fortran {!example/linalg/example_eye1.f90!} ``` diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 2f14e512d..31a96e10e 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1342,7 +1342,7 @@ contains integer, intent(in) :: dim1 integer, intent(in), optional :: dim2 - ${t1}$, intent(in) #{if 'int8' in t1}#, optional #{endif}#:: mold + ${t1}$, intent(in) #{if t1 == 'real(dp)'}#, optional #{endif}#:: mold ${t1}$, allocatable :: result(:, :) integer :: dim2_