Choral
algebra_lin.F90
Go to the documentation of this file.
1 
8 
9 module algebra_lin
10 
11  use real_type, only: rp
12  use basic_tools
13 
14  implicit none
15  private
16 
17  ! %----------------------------------------%
18  ! | |
19  ! | PUBLIC DATA |
20  ! | |
21  ! %----------------------------------------%
22 
23  !!
24  !! LINEAR ALGEBRA
25  !!
26  !! GENERIC
27  public :: transpose ! tested
28  public :: matvecprod ! tested
29  !!
30  !! NON GENERIC
31  public :: matmatprod ! tested
32  public :: invtrisup ! tested
33 
34  !!
35  !! LINEAR ALGEBRA in R**3
36  !!
37  !!
38  public :: nrm2_3d ! tested
39  public :: matvecprod_3x3 ! tested
40  public :: matvecprod_3x2 ! tested
41  public :: matvecprod_3x1 ! tested
42  public :: solve_3x3 ! tested
43  public :: det_3x3 ! tested
44  public :: crossprod_3d ! tested
45  public :: cplmtmat_3x2 ! tested
46  public :: cplmtmat_3x3 ! tested
47  public :: det_3x2 ! tested
48 
49 
50  ! %----------------------------------------%
51  ! | |
52  ! | GENERIc SUBROUTINES |
53  ! | |
54  ! %----------------------------------------%
55 
56  interface matvecprod
57  module procedure algebra_lin_matvecprod !! TESTED
58  end interface matvecprod
59 
60  interface transpose
61  module procedure array_transpose !! TESTED
62  end interface transpose
63 
64 contains
65 
66 
70  subroutine algebra_lin_matvecprod(y,A,X)
71  real(RP), dimension(:) , intent(in) :: X
72  real(RP), dimension(:) , intent(out):: y
73  real(RP), dimension(size(y,1), size(x,1)), intent(in) :: A
74 
75  integer :: ii
76 
77 #if DBG>0
78  if ( size(a,2)/= size(x,1) ) &
79  & call quit( &
80  & 'algebra_lin: algebra_lin_matVecProd: 1' )
81  if ( size(a,1)/= size(y,1) ) &
82  & call quit( &
83  & 'algebra_lin: algebra_lin_matVecProd: 2' )
84 #endif
85 
86  y = x(1)*a(:,1)
87  do ii=2, size(x,1)
88  y = y + x(ii)*a(:,ii)
89  end do
90 
91  end subroutine algebra_lin_matvecprod
92 
94  subroutine matmatprod(M, A, B)
95  real(RP), dimension(:,:), intent(out) :: M
96  real(RP), dimension(:,:), intent(in) :: A, B
97 
98  integer :: ii, jj
99 
100 #if DBG>0
101  if ( size(a,2)/= size(b,1) ) &
102  & call quit( &
103  & 'algebra_lin: matMatProd: 1' )
104  if ( size(a,1)/= size(m,1) ) &
105  & call quit( &
106  & 'algebra_lin: matMatProd: 2' )
107  if ( size(b,2)/= size(m,2) ) &
108  & call quit( &
109  & 'algebra_lin: matMatProd: 3' )
110 #endif
111 
112  do ii=1, size(a,1)
113  do jj=1, size(b,2)
114  m(ii,jj) = sum(a(ii,:)*b(:,jj))
115  end do
116  end do
117 
118  end subroutine matmatprod
119 
121  subroutine array_transpose(M, A)
122  real(RP), dimension(:,:), intent(out) :: M
123  real(RP), dimension(:,:), intent(in) :: A
124 
125  integer :: ii, jj
126 
127 #if DBG>0
128  if ( size(a,1)/= size(m,2) ) call quit( &
129  & 'algebra_lin: array_transpose: 1' )
130  if ( size(a,2)/= size(m,1) ) call quit( &
131  & 'algebra_lin: array_transpose: 2' )
132 #endif
133 
134  do ii=1, size(a,1)
135  do jj=1, size(a,2)
136  m(jj,ii) = a(ii,jj)
137  end do
138  end do
139 
140  end subroutine array_transpose
141 
142 
145  subroutine invtrisup(res,M,vec)
147  real(RP), dimension(:,:), intent(in) :: M
148  real(RP), dimension(:) , intent(in) :: vec
149  real(RP), dimension(:) , intent(out) :: res
150 
151  integer :: ii,ji,ni
152 
153 #if DBG>0
154  if ( size(m,1)/= size(m,2) ) call quit( &
155  & 'algebra_lin: invtrisup: 1' )
156  if ( size(m,1)/= size(vec,1) ) call quit( &
157  & 'algebra_lin: invtrisup: 2' )
158  if ( size(m,1)/= size(res,1) ) call quit( &
159  & 'algebra_lin: invtrisup: 3' )
160 #endif
161 
162  ni=size(m,1)
163  res(ni)=vec(ni)/m(ni,ni)
164 
165  do ii=ni-1,1,-1
166  res(ii)=vec(ii)
167  do ji=ii+1,ni
168  res(ii)=res(ii)-m(ii,ji)*res(ji)
169  end do
170  res(ii)=res(ii)/m(ii,ii)
171  end do
172 
173  end subroutine invtrisup
174 
177  function crossprod_3d(u,v) result(res)
178  real(RP), dimension(3), intent(in) :: u
179  real(RP), dimension(3), intent(in) :: v
180  real(RP), dimension(3) :: res
181 
182  res(1) = u(2)*v(3)-u(3)*v(2)
183  res(2) = -u(1)*v(3)+u(3)*v(1)
184  res(3) = u(1)*v(2)-u(2)*v(1)
185 
186  end function crossprod_3d
187 
190  function det_3x3(M) result(d)
191  real(RP), dimension(3,3), intent(in) :: M
192  real(RP) :: d
193 
194  real(RP) :: u, v, w
195 
196  u = m(2,2)*m(3,3)-m(3,2)*m(2,3)
197  v = m(1,2)*m(3,3)-m(1,3)*m(3,2)
198  w = m(1,2)*m(2,3)-m(1,3)*m(2,2)
199 
200  d = m(1,1)*u - m(2,1)*v + m(3,1)*w
201 
202  end function det_3x3
203 
204 
209  function det_3x2(M) result(d)
210  real(RP), dimension(3,2), intent(in) :: M
211  real(RP) :: d
212 
213  real(RP) :: u
214 
215  u = m(2,1)*m(3,2)-m(3,1)*m(2,2)
216  d = u**2
217  u = -m(1,1)*m(3,2)+m(3,1)*m(1,2)
218  d = d + u**2
219  u = m(1,1)*m(2,2)-m(2,1)*m(1,2)
220  d = d + u**2
221  d = sqrt(d)
222 
223  end function det_3x2
224 
225 
226 
232  function cplmtmat_3x3(M) result(N)
233  real(RP), dimension(3,3) :: N
234  real(RP), dimension(3,3), intent(in) :: M
235 
236  n(:,1) = crossprod_3d(m(:,2), m(:,3))
237  n(:,2) = crossprod_3d(m(:,3), m(:,1))
238  n(:,3) = crossprod_3d(m(:,1), m(:,2))
239 
240  end function cplmtmat_3x3
241 
246  function cplmtmat_3x2(M) result(N)
247  real(RP), dimension(3,2) :: N
248  real(RP), dimension(3,2), intent(in) :: M
249 
250  real(RP), dimension(3) :: e1, e2, e3
251  real(RP):: xB, xC, yC
252 
253  e1 = m(:,1)
254  e2 = m(:,2)
255  e3 = crossprod_3d( e1, e2 )
256  e2 = crossprod_3d(e3, e1)
257  e1 = e1 / nrm2_3d(e1)
258  e2 = e2 / nrm2_3d(e2)
259  xb = dot_product(m(:,1), e1)
260  xc = dot_product(m(:,2), e1)
261  yc = dot_product(m(:,2), e2)
262 
263  n(:,1) = yc*e1 - xc*e2
264  n(:,2) = xb*e2
265 
266  end function cplmtmat_3x2
267 
270  function nrm2_3d(v) result(x)
271  real(RP), dimension(3) :: v
272  real(RP) :: x
273 
274  x = sum(v*v)
275  x = sqrt(x)
276 
277  end function nrm2_3d
278 
281  subroutine matvecprod_3x3(y,A,X)
282  real(RP), dimension(3) , intent(out):: y
283  real(RP), dimension(3) , intent(in) :: X
284  real(RP), dimension(3,3), intent(in) :: A
285 
286  y(1:3) = a(1:3, 1) * x(1)
287  y(1:3) = y(1:3) + a(1:3, 2) * x(2)
288  y(1:3) = y(1:3) + a(1:3, 3) * x(3)
289 
290  end subroutine matvecprod_3x3
291 
292 
293 
296  subroutine matvecprod_3x2(y,A,X)
297  real(RP), dimension(3) , intent(out):: y
298  real(RP), dimension(2) , intent(in) :: X
299  real(RP), dimension(3,2), intent(in) :: A
300 
301  y(1:3) = a(1:3, 1) * x(1)
302  y(1:3) = y(1:3) + a(1:3, 2) * x(2)
303 
304  end subroutine matvecprod_3x2
305 
306 
309  subroutine matvecprod_3x1(y,A,X)
310  real(RP), dimension(3) , intent(out):: y
311  real(RP), dimension(1) , intent(in) :: X
312  real(RP), dimension(3,1), intent(in) :: A
313 
314  y(1:3) = a(1:3, 1) * x(1)
315 
316  end subroutine matvecprod_3x1
317 
318 
321  subroutine solve_3x3(x, A, b)
322  real(RP), dimension(3) , intent(out) :: x
323  real(RP), dimension(3,3), intent(in) :: A
324  real(RP), dimension(3) , intent(in) :: b
325 
326  real(RP), dimension(3) :: v
327 
328  v = crossprod_3d(a(:,2),a(:,3))
329  x(1) = dot_product(b,v) / dot_product(a(:,1),v)
330 
331  v = crossprod_3d(a(:,1),a(:,3))
332  x(2) = dot_product(b,v) / dot_product(a(:,2),v)
333 
334  v = crossprod_3d(a(:,1),a(:,2))
335  x(3) = dot_product(b,v) / dot_product(a(:,3),v)
336 
337  end subroutine solve_3x3
338 
339 end module algebra_lin
340 
341 
subroutine algebra_lin_matvecprod(y, A, X)
Matrix Vector Product v <= m*u.
Definition: algebra_lin.F90:71
LINEAR ALGEBRA TOOLS
Definition: algebra_lin.F90:9
real(rp) function, dimension(3, 2), public cplmtmat_3x2(M)
N = &#39;complementary&#39; matrix of M.
BASIC TOOLS
Definition: basic_tools.F90:8
subroutine, public matvecprod_3x2(y, A, X)
matrix vector product R2 –> R3
subroutine, public matvecprod_3x3(y, A, X)
matrix vector product R3 –> R3
integer, parameter, public rp
real(kind=RP) = real precision in the code REAL_TOL = epsilon to test real equality ...
Definition: real_type.F90:90
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine, public solve_3x3(x, A, b)
Solve A x = b for 3x3 matrix.
subroutine array_transpose(M, A)
M <– transpose(A)
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
subroutine, public invtrisup(res, M, vec)
Calcul de x=M^{-1}vec avec M triangulaire superieure.
subroutine, public matvecprod_3x1(y, A, X)
matrix vector product R1 –> R3
real(rp) function, dimension(3, 3), public cplmtmat_3x3(M)
N = complementary matrix of M = transpose of the cofactor matrix.
real(rp) function, public det_3x3(M)
determinant of 3x3 matrix
real(rp) function, dimension(3), public crossprod_3d(u, v)
Cross product in R**3.
real(rp) function, public nrm2_3d(v)
Euclidian norm.
real(rp) function, public det_3x2(M)
&#39;determinant&#39; of 3x2 matrix
subroutine, public matmatprod(M, A, B)
M <– A * B.
Definition: algebra_lin.F90:95