1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
module matrix_m
use monoid_m, only: derive_extended_monoid
use semiring_m, only: semiring
use unit_ring_m, only: unit_ring_only_minus, derive_unit_ring_from_minus
use field_m, only: field_only_division
implicit none
private
public :: matrix_tmpl
template matrix_tmpl(T, plus_t, zero_t, times_t, one_t, n)
require :: semiring(T, plus_t, zero_t, times_t, one_t)
instantiate derive_extended_monoid(T, plus_t, zero_t), only: sum => mconcat
integer :: n
private
public :: &
matrix, &
plus_matrix, &
times_matrix, &
zero, &
one, &
matrix_subtraction_tmpl
type :: matrix
type(T) :: elements(n, n)
end type
interface operator(+)
procedure :: plus_matrix
end interface
interface operator(*)
procedure times_matrix
end interface
template matrix_subtraction_tmpl(minus_t)
require :: unit_ring_only_minus(T, plus_t, zero_t, times_t, one_t, minus_t)
private
public :: operator(-), gaussian_solver_tmpl
interface operator(-)
procedure minus_matrix
end interface
template gaussian_solver_tmpl(div_t)
instantiate derive_unit_ring_from_minus(T, plus_t, zero_t, times_t, one_t, minus_t), only: negate
require :: field_only_division(T, plus_t, zero_t, times_t, one_t, minus_t, negate, div_t)
private
contains
elemental function div_matrix(x, y) result(quotient)
type(matrix), intent(in) :: x, y
type(matrix) :: quotient
quotient = back_substitute(row_eschelon(x), y)
end function
pure function row_eschelon(x) result(reduced)
type(matrix), intent(in) :: x
type(matrix) :: reduced
integer :: i, ii, j
type(T) :: r
reduced = x
do i = 1, n
! Assume pivot m(i,i) is not zero
do ii = i+1, n
r = div_t(reduced%elements(i,i), reduced%elements(ii,i))
reduced%elements(ii, i) = zero_t()
do j = i+1, n
reduced%elements(ii, j) = minus_t(reduced%elements(ii, j), times_t(reduced%elements(i, j), r))
end do
end do
end do
end function
pure function back_substitute(x, y) result(solved)
type(matrix), intent(in) :: x, y
type(matrix) :: solved
integer :: i, j
type(T) :: tmp(n)
solved = y
do i = n, 1, -1
tmp = zero_t()
do j = i+1, n
tmp = plus_t(tmp, times_t(x%elements(i,j), solved%elements(:,j)))
end do
solved%elements(:,i) = div_t(minus_t(solved%elements(:, i), tmp), x%elements(i,i))
end do
end function
end template
contains
elemental function minus_matrix(x, y) result(difference)
type(matrix), intent(in) :: x, y
type(matrix) :: difference
difference%elements = minus_t(x%elements, y%elements)
end function
end template
contains
elemental function plus_matrix(x, y) result(combined)
type(matrix), intent(in) :: x, y
type(matrix) :: combined
combined%elements = plus_t(x%elements, y%elements)
end function
pure function zero()
type(matrix) :: zero
zero%elements = zero_t()
end function
elemental function times_matrix(x, y) result(combined)
type(matrix), intent(in) :: x, y
type(matrix) :: combined
integer :: i, j, k
type(T) :: dot
do i = 1, n
do j = 1, n
combined%elements(i, j) = sum(times_t(x%elements(i,:), y%elements(:,j)))
end do
end do
end function
pure function one()
type(matrix) :: one
integer :: i
one%elements = zero_t()
do concurrent (i = 1:n)
one%elements(i, i) = one_t()
end do
end function
end template
end module
|
