Skip to content

Commit bbecb52

Browse files
ivan-piivan-pi
authored
Merge pull request #170 from ivan-pi/ivan-pi
Addition of diag, eye, and trace
2 parents c95f7a9 + d01f657 commit bbecb52

8 files changed

+771
-0
lines changed

‎src/CMakeLists.txt

+2
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,8 @@
33
# Create a list of the files to be preprocessed
44
set(fppFiles
55
stdlib_experimental_io.fypp
6+
stdlib_experimental_linalg.fypp
7+
stdlib_experimental_linalg_diag.fypp
68
stdlib_experimental_optval.fypp
79
stdlib_experimental_stats.fypp
810
stdlib_experimental_stats_cov.fypp

‎src/Makefile.manual

+5
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,8 @@ SRC = f18estop.f90 \
22
stdlib_experimental_ascii.f90 \
33
stdlib_experimental_error.f90 \
44
stdlib_experimental_io.f90 \
5+
stdlib_experimental_linalg.f90 \
6+
stdlib_experimental_linalg_diag.f90 \
57
stdlib_experimental_kinds.f90 \
68
stdlib_experimental_optval.f90 \
79
stdlib_experimental_quadrature.f90 \
@@ -42,6 +44,7 @@ stdlib_experimental_io.o: \
4244
stdlib_experimental_error.o \
4345
stdlib_experimental_optval.o \
4446
stdlib_experimental_kinds.o
47+
stdlib_experimental_linalg_diag.o: stdlib_experimental_kinds.o
4548
stdlib_experimental_optval.o: stdlib_experimental_kinds.o
4649
stdlib_experimental_quadrature.o: stdlib_experimental_kinds.o
4750
stdlib_experimental_stats_mean.o: \
@@ -59,6 +62,8 @@ stdlib_experimental_stats_var.o: \
5962

6063
# Fortran sources that are built from fypp templates
6164
stdlib_experimental_io.f90: stdlib_experimental_io.fypp
65+
stdlib_experimental_linalg.f90: stdlib_experimental_linalg.fypp
66+
stdlib_experimental_linalg_diag.f90: stdlib_experimental_linalg_diag.fypp
6267
stdlib_experimental_quadrature.f90: stdlib_experimental_quadrature.fypp
6368
stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp
6469
stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp

‎src/stdlib_experimental_linalg.fypp

+80
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
#:include "common.fypp"
2+
#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
3+
module stdlib_experimental_linalg
4+
use stdlib_experimental_kinds, only: sp, dp, qp, &
5+
int8, int16, int32, int64
6+
implicit none
7+
private
8+
9+
public :: diag
10+
public :: eye
11+
public :: trace
12+
13+
interface diag
14+
!
15+
! Vector to matrix
16+
!
17+
#:for k1, t1 in RCI_KINDS_TYPES
18+
module function diag_${t1[0]}$${k1}$(v) result(res)
19+
${t1}$, intent(in) :: v(:)
20+
${t1}$ :: res(size(v),size(v))
21+
end function diag_${t1[0]}$${k1}$
22+
#:endfor
23+
#:for k1, t1 in RCI_KINDS_TYPES
24+
module function diag_${t1[0]}$${k1}$_k(v,k) result(res)
25+
${t1}$, intent(in) :: v(:)
26+
integer, intent(in) :: k
27+
${t1}$ :: res(size(v)+abs(k),size(v)+abs(k))
28+
end function diag_${t1[0]}$${k1}$_k
29+
#:endfor
30+
31+
!
32+
! Matrix to vector
33+
!
34+
#:for k1, t1 in RCI_KINDS_TYPES
35+
module function diag_${t1[0]}$${k1}$_mat(A) result(res)
36+
${t1}$, intent(in) :: A(:,:)
37+
${t1}$ :: res(minval(shape(A)))
38+
end function diag_${t1[0]}$${k1}$_mat
39+
#:endfor
40+
#:for k1, t1 in RCI_KINDS_TYPES
41+
module function diag_${t1[0]}$${k1}$_mat_k(A,k) result(res)
42+
${t1}$, intent(in) :: A(:,:)
43+
integer, intent(in) :: k
44+
${t1}$ :: res(minval(shape(A))-abs(k))
45+
end function diag_${t1[0]}$${k1}$_mat_k
46+
#:endfor
47+
end interface
48+
49+
! Matrix trace
50+
interface trace
51+
#:for k1, t1 in RCI_KINDS_TYPES
52+
module procedure trace_${t1[0]}$${k1}$
53+
#:endfor
54+
end interface
55+
56+
contains
57+
58+
function eye(n) result(res)
59+
integer, intent(in) :: n
60+
integer(int8) :: res(n, n)
61+
integer :: i
62+
res = 0
63+
do i = 1, n
64+
res(i, i) = 1
65+
end do
66+
end function eye
67+
68+
69+
#:for k1, t1 in RCI_KINDS_TYPES
70+
function trace_${t1[0]}$${k1}$(A) result(res)
71+
${t1}$, intent(in) :: A(:,:)
72+
${t1}$ :: res
73+
integer :: i
74+
res = 0
75+
do i = 1, minval(shape(A))
76+
res = res + A(i,i)
77+
end do
78+
end function trace_${t1[0]}$${k1}$
79+
#:endfor
80+
end module

‎src/stdlib_experimental_linalg.md

+156
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,156 @@
1+
# Linear Algebra
2+
3+
* [`diag` - Create a diagonal array or extract the diagonal elements of an array](#diag---create-a-diagonal-array-or-extract-the-diagonal-elements-of-an-array)
4+
* [`eye` - Construct the identity matrix](#eye---construct-the-identity-matrix)
5+
* [`trace` - Trace of a matrix](#trace---trace-of-a-matrix)
6+
7+
## `diag` - Create a diagonal array or extract the diagonal elements of an array
8+
9+
### Description
10+
11+
Create a diagonal array or extract the diagonal elements of an array
12+
13+
### Syntax
14+
15+
`d = diag(a [, k])`
16+
17+
### Arguments
18+
19+
`a`: Shall be a rank-1 or or rank-2 array. If `a` is a rank-1 array (i.e. a vector) then `diag` returns a rank-2 array with the elements of `a` on the diagonal. If `a` is a rank-2 array (i.e. a matrix) then `diag` returns a rank-1 array of the diagonal elements.
20+
21+
`k` (optional): Shall be a scalar of type `integer` and specifies the diagonal. The default `k = 0` represents the main diagonal, `k > 0` are diagonals above the main diagonal, `k < 0` are diagonals below the main diagonal.
22+
23+
### Return value
24+
25+
Returns a diagonal array or a vector with the extracted diagonal elements.
26+
27+
### Example
28+
29+
```fortran
30+
program demo_diag1
31+
use stdlib_experimental_linalg, only: diag
32+
implicit none
33+
real, allocatable :: A(:,:)
34+
integer :: i
35+
A = diag([(1,i=1,10)]) ! creates a 10 by 10 identity matrix
36+
end program demo_diag1
37+
```
38+
39+
```fortran
40+
program demo_diag2
41+
use stdlib_experimental_linalg, only: diag
42+
implicit none
43+
real :: v(:)
44+
real, allocatable :: A(:,:)
45+
integer :: i
46+
v = [1,2,3,4,5]
47+
A = diag(v) ! creates a 5 by 5 matrix with elements of v on the diagonal
48+
end program demo_diag2
49+
```
50+
51+
```fortran
52+
program demo_diag3
53+
use stdlib_experimental_linalg, only: diag
54+
implicit none
55+
integer, parameter :: n = 10
56+
real :: c(n), ul(n-1)
57+
real :: A(n,n)
58+
integer :: i
59+
c = 2
60+
ul = -1
61+
A = diag(ul,-1) + diag(c) + diag(ul,1) ! Gil Strang's favorite matrix
62+
end program demo_diag3
63+
```
64+
65+
```fortran
66+
program demo_diag4
67+
use stdlib_experimental_linalg, only: diag
68+
implicit none
69+
integer, parameter :: n = 12
70+
real :: A(n,n)
71+
real :: v(n)
72+
integer :: i
73+
call random_number(A)
74+
v = diag(A) ! v contains diagonal elements of A
75+
end program demo_diag4
76+
```
77+
78+
```fortran
79+
program demo_diag5
80+
use stdlib_experimental_linalg, only: diag
81+
implicit none
82+
integer, parameter :: n = 3
83+
real :: A(n,n)
84+
real, allocatable :: v(:)
85+
integer :: i
86+
A = reshape([1,2,3,4,5,6,7,8,9],[n,n])
87+
v = diag(A,-1) ! v is [2,6]
88+
v = diag(A,1) ! v is [4,8]
89+
end program demo_diag5
90+
```
91+
92+
## `eye` - Construct the identity matrix
93+
94+
### Description
95+
96+
Construct the identity matrix
97+
98+
## Syntax
99+
100+
`I = eye(n)`
101+
102+
### Arguments
103+
104+
`n`: Shall be a scalar of default type `integer`.
105+
106+
### Return value
107+
108+
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)`.
109+
110+
### Example
111+
112+
```fortran
113+
program demo_eye1
114+
use stdlib_experimental_linalg, only: eye
115+
implicit none
116+
real :: a(3,3)
117+
A = eye(3)
118+
end program demo_eye1
119+
```
120+
121+
```fortran
122+
program demo_eye2
123+
use stdlib_experimental_linalg, only: eye, diag
124+
implicit none
125+
print *, all(eye(4) == diag([1,1,1,1])) ! prints .true.
126+
end program demo_eye2
127+
```
128+
129+
## `trace` - Trace of a matrix
130+
131+
### Description
132+
133+
Trace of a matrix (rank-2 array)
134+
135+
### Syntax
136+
137+
`result = trace(A)`
138+
139+
### Arguments
140+
141+
`A`: Shall be a rank-2 array. If `A` is not square, then `trace(A)` will return the sum of diagonal values from the square sub-section of `A`.
142+
143+
### Return value
144+
145+
Returns the trace of the matrix, i.e. the sum of diagonal elements.
146+
147+
### Example
148+
```fortran
149+
program demo_trace
150+
use stdlib_experimental_linalg, only: trace
151+
implicit none
152+
real :: A(3,3)
153+
A = reshape([1,2,3,4,5,6,7,8,9],[3,3])
154+
print *, trace(A) ! 1 + 5 + 9
155+
end program demo_trace
156+
```

‎src/stdlib_experimental_linalg_diag.fypp

+80
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
#:include "common.fypp"
2+
#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
3+
submodule (stdlib_experimental_linalg) stdlib_experimental_linalg_diag
4+
5+
implicit none
6+
7+
contains
8+
9+
#:for k1, t1 in RCI_KINDS_TYPES
10+
function diag_${t1[0]}$${k1}$(v) result(res)
11+
${t1}$, intent(in) :: v(:)
12+
${t1}$ :: res(size(v),size(v))
13+
integer :: i
14+
res = 0
15+
do i = 1, size(v)
16+
res(i,i) = v(i)
17+
end do
18+
end function diag_${t1[0]}$${k1}$
19+
#:endfor
20+
21+
22+
#:for k1, t1 in RCI_KINDS_TYPES
23+
function diag_${t1[0]}$${k1}$_k(v,k) result(res)
24+
${t1}$, intent(in) :: v(:)
25+
integer, intent(in) :: k
26+
${t1}$ :: res(size(v)+abs(k),size(v)+abs(k))
27+
integer :: i, sz
28+
sz = size(v)
29+
res = 0
30+
if (k > 0) then
31+
do i = 1, sz
32+
res(i,k+i) = v(i)
33+
end do
34+
else if (k < 0) then
35+
do i = 1, sz
36+
res(i+abs(k),i) = v(i)
37+
end do
38+
else
39+
do i = 1, sz
40+
res(i,i) = v(i)
41+
end do
42+
end if
43+
end function diag_${t1[0]}$${k1}$_k
44+
#:endfor
45+
46+
#:for k1, t1 in RCI_KINDS_TYPES
47+
function diag_${t1[0]}$${k1}$_mat(A) result(res)
48+
${t1}$, intent(in) :: A(:,:)
49+
${t1}$ :: res(minval(shape(A)))
50+
integer :: i
51+
do i = 1, minval(shape(A))
52+
res(i) = A(i,i)
53+
end do
54+
end function diag_${t1[0]}$${k1}$_mat
55+
#:endfor
56+
57+
#:for k1, t1 in RCI_KINDS_TYPES
58+
function diag_${t1[0]}$${k1}$_mat_k(A,k) result(res)
59+
${t1}$, intent(in) :: A(:,:)
60+
integer, intent(in) :: k
61+
${t1}$ :: res(minval(shape(A))-abs(k))
62+
integer :: i, sz
63+
sz = minval(shape(A))-abs(k)
64+
if (k > 0) then
65+
do i = 1, sz
66+
res(i) = A(i,k+i)
67+
end do
68+
else if (k < 0) then
69+
do i = 1, sz
70+
res(i) = A(i+abs(k),i)
71+
end do
72+
else
73+
do i = 1, sz
74+
res(i) = A(i,i)
75+
end do
76+
end if
77+
end function diag_${t1[0]}$${k1}$_mat_k
78+
#:endfor
79+
80+
end submodule

‎src/tests/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ endmacro(ADDTEST)
88

99
add_subdirectory(ascii)
1010
add_subdirectory(io)
11+
add_subdirectory(linalg)
1112
add_subdirectory(optval)
1213
add_subdirectory(stats)
1314
add_subdirectory(system)

‎src/tests/linalg/CMakeLists.txt

+2
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
ADDTEST(linalg)
2+

0 commit comments

Comments
 (0)