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
78 if (
size(a,2)/=
size(x,1) ) &
80 &
'algebra_lin: algebra_lin_matVecProd: 1' )
81 if (
size(a,1)/=
size(y,1) ) &
83 &
'algebra_lin: algebra_lin_matVecProd: 2' )
95 real(RP),
dimension(:,:),
intent(out) :: M
96 real(RP),
dimension(:,:),
intent(in) :: A, B
101 if (
size(a,2)/=
size(b,1) ) &
103 &
'algebra_lin: matMatProd: 1' )
104 if (
size(a,1)/=
size(m,1) ) &
106 &
'algebra_lin: matMatProd: 2' )
107 if (
size(b,2)/=
size(m,2) ) &
109 &
'algebra_lin: matMatProd: 3' )
114 m(ii,jj) = sum(a(ii,:)*b(:,jj))
122 real(RP),
dimension(:,:),
intent(out) :: M
123 real(RP),
dimension(:,:),
intent(in) :: A
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' )
147 real(RP),
dimension(:,:),
intent(in) :: M
148 real(RP),
dimension(:) ,
intent(in) :: vec
149 real(RP),
dimension(:) ,
intent(out) :: res
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' )
163 res(ni)=vec(ni)/m(ni,ni)
168 res(ii)=res(ii)-m(ii,ji)*res(ji)
170 res(ii)=res(ii)/m(ii,ii)
178 real(RP),
dimension(3),
intent(in) :: u
179 real(RP),
dimension(3),
intent(in) :: v
180 real(RP),
dimension(3) :: res
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)
191 real(RP),
dimension(3,3),
intent(in) :: M
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)
200 d = m(1,1)*u - m(2,1)*v + m(3,1)*w
210 real(RP),
dimension(3,2),
intent(in) :: M
215 u = m(2,1)*m(3,2)-m(3,1)*m(2,2)
217 u = -m(1,1)*m(3,2)+m(3,1)*m(1,2)
219 u = m(1,1)*m(2,2)-m(2,1)*m(1,2)
233 real(RP),
dimension(3,3) :: N
234 real(RP),
dimension(3,3),
intent(in) :: M
247 real(RP),
dimension(3,2) :: N
248 real(RP),
dimension(3,2),
intent(in) :: M
250 real(RP),
dimension(3) :: e1, e2, e3
251 real(RP):: xB, xC, yC
259 xb = dot_product(m(:,1), e1)
260 xc = dot_product(m(:,2), e1)
261 yc = dot_product(m(:,2), e2)
263 n(:,1) = yc*e1 - xc*e2
271 real(RP),
dimension(3) :: v
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
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)
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
301 y(1:3) = a(1:3, 1) * x(1)
302 y(1:3) = y(1:3) + a(1:3, 2) * x(2)
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
314 y(1:3) = a(1:3, 1) * x(1)
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
326 real(RP),
dimension(3) :: v
329 x(1) = dot_product(b,v) / dot_product(a(:,1),v)
332 x(2) = dot_product(b,v) / dot_product(a(:,2),v)
335 x(3) = dot_product(b,v) / dot_product(a(:,3),v)
subroutine algebra_lin_matvecprod(y, A, X)
Matrix Vector Product v <= m*u.
real(rp) function, dimension(3, 2), public cplmtmat_3x2(M)
N = 'complementary' matrix of M.
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 ...
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
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)
'determinant' of 3x2 matrix
subroutine, public matmatprod(M, A, B)
M <– A * B.