GCC Code Coverage Report

Directory:	src/fortran/lib/
File:	mod_misc_linalg.f90
Date:	2025-04-05 12:17:58

	Exec	Total	Coverage
Lines:	51	51	100.0%
Functions:	0	0	-%
Branches:	38	42	90.5%

  
      Line
      Branch
      Exec
      Source
    
      module raffle__misc_linalg
    
        !! Module contains various linear algebra functions and subroutines.
    
        use raffle__constants, only: real32, pi
    
        implicit none
    
        private
    
        public :: modu, cross
    
        public :: get_distance, get_angle, get_dihedral_angle
    
        public :: get_improper_dihedral_angle
    
        public :: inverse_3x3
    
        interface get_angle
    
           procedure get_angle_from_points, get_angle_from_vectors
    
        end interface get_angle
    
        interface get_dihedral_angle
    
           procedure get_dihedral_angle_from_points, get_dihedral_angle_from_vectors
    
        end interface get_dihedral_angle
    
        interface get_improper_dihedral_angle
    
           procedure get_improper_dihedral_angle_from_points, &
    
                get_improper_dihedral_angle_from_vectors
    
        end interface get_improper_dihedral_angle
    
      contains
    
      !###############################################################################
    
        1/2✓ Branch 0 taken 422656602 times.
✗ Branch 1 not taken.

      422656602
        pure function modu(vector)
    
          !! Return the magnitude of a vector of any size.
    
          implicit none
    
          ! Arguments
    
          real(real32),dimension(:), intent(in) :: vector
    
          !! Input vector.
    
          real(real32) :: modu
    
          !! Output magnitude.
    
        2/2✓ Branch 0 taken 1267969806 times.
✓ Branch 1 taken 422656602 times.

      1690626408
          modu = abs( sqrt( sum(vector(:)**2) ) )
    
      422656602
        end function modu
    
      !###############################################################################
    
      !###############################################################################
    
        1/2✓ Branch 0 taken 239223063 times.
✗ Branch 1 not taken.

      239223063
        pure function cross(a,b)
    
          !! Return the cross product of two vectors.
    
          implicit none
    
          ! Arguments
    
          real(real32), dimension(3), intent(in) :: a, b
    
          !! Input vectors.
    
          real(real32), dimension(3) :: cross
    
          !! Output cross product.
    
      239223063
          cross(1) = a(2) * b(3) - a(3) * b(2)
    
      239223063
          cross(2) = a(3) * b(1) - a(1) * b(3)
    
      239223063
          cross(3) = a(1) * b(2) - a(2) * b(1)
    
        end function cross
    
      !###############################################################################
    
      !###############################################################################
    
      1
        pure function get_distance(point1, point2) result(distance)
    
          !! Return the distance between two points.
    
          implicit none
    
          ! Arguments
    
          real(real32), dimension(3), intent(in) :: point1, point2
    
          !! Input points.
    
          real(real32) :: distance
    
          !! Output distance.
    
        2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.

      4
          distance = modu( point1 - point2 )
    
      1
          return
    
        end function get_distance
    
      !###############################################################################
    
      !###############################################################################
    
      119453303
        pure function get_angle_from_vectors(vector1, vector2) result(angle)
    
          !! Return the angle between two vectors.
    
          implicit none
    
          ! Arguments
    
          real(real32), dimension(3), intent(in) :: vector1, vector2
    
          !! Input vectors.
    
          real(real32) :: angle
    
          !! Output angle.
    
          angle =  dot_product(vector1,vector2) / &
    
        2/2✓ Branch 0 taken 358359909 times.
✓ Branch 1 taken 119453303 times.

      477813212
               ( modu(vector1) * modu(vector2) )
    
        2/2✓ Branch 0 taken 1785068 times.
✓ Branch 1 taken 117668235 times.

      119453303
          if(angle .ge. 1._real32)then
    
      1785068
             angle = 0._real32
    
        2/2✓ Branch 0 taken 1750023 times.
✓ Branch 1 taken 115918212 times.

      117668235
          elseif(angle .le. -1._real32)then
    
      1750023
             angle = pi
    
          else
    
      115918212
             angle = acos(angle)
    
          end if
    
      119453303
        end function get_angle_from_vectors
    
      !###############################################################################
    
      !###############################################################################
    
      3958532
        pure function get_angle_from_points(point1, point2, point3) result(angle)
    
          !! Return the angle formed by three points.
    
          !!
    
          !! The angle is formed by the path point1 -> point2 -> point3.
    
          implicit none
    
          ! Arguments
    
          real(real32), dimension(3), intent(in) :: point1, point2, point3
    
          !! Input points.
    
          real(real32) :: angle
    
          !! Output angle.
    
          angle = dot_product( point1 - point2, point3 - point2 ) / &
    
        6/6✓ Branch 0 taken 11875596 times.
✓ Branch 1 taken 3958532 times.
✓ Branch 2 taken 11875596 times.
✓ Branch 3 taken 3958532 times.
✓ Branch 4 taken 11875596 times.
✓ Branch 5 taken 3958532 times.

      39585320
               ( modu( point1 - point2 ) * modu( point3 - point2 ) )
    
        2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 3958520 times.

      3958532
          if(angle .ge. 1._real32)then
    
      12
             angle = 0._real32
    
        2/2✓ Branch 0 taken 2689 times.
✓ Branch 1 taken 3955831 times.

      3958520
          elseif(angle .le. -1._real32)then
    
      2689
             angle = pi
    
          else
    
      3955831
             angle = acos(angle)
    
          end if
    
      3958532
        end function get_angle_from_points
    
      !###############################################################################
    
      !###############################################################################
    
      1
        pure function get_dihedral_angle_from_vectors( &
    
             vector1, vector2, vector3) result(angle)
    
          !! Return the dihedral angle between two planes.
    
          !!
    
          !! The dihedral angle is the angle between the plane defined by the cross
    
          !! product of two vectors and a third vector.
    
          !! i.e. ( vector1 x vector2 ) . vector3
    
          implicit none
    
          ! Arguments
    
          real(real32), dimension(3), intent(in) :: vector1, vector2, vector3
    
          !! Input vectors.
    
          real(real32) :: angle
    
          !! Output angle.
    
      1
          angle = get_angle(cross(vector1, vector2), vector3)
    
      1
        end function get_dihedral_angle_from_vectors
    
      !###############################################################################
    
      !###############################################################################
    
      1
        pure function get_dihedral_angle_from_points( &
    
             point1, point2, point3, point4 &
    
        ) result(angle)
    
          !! Return the dihedral angle between two planes.
    
          !!
    
          !! The dihedral angle is the angle between the plane defined by four points.
    
          !! i.e. ( point2 - point1 ) x ( point3 - point2 ) . ( point4 - point2 )
    
          !! alt. angle between plane point1point2point3 and vector point2point4
    
          implicit none
    
          real(real32), dimension(3), intent(in) :: point1, point2, point3, point4
    
          real(real32) :: angle
    
        6/6✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 1 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 1 times.

      10
          angle = get_angle(cross(point2 - point1, point3 - point2), point4 - point2)
    
      1
        end function get_dihedral_angle_from_points
    
      !###############################################################################
    
      !###############################################################################
    
      278094
        pure function get_improper_dihedral_angle_from_vectors( &
    
             vector1, vector2, vector3 &
    
        ) result(angle)
    
          !! Return the improper dihedral angle between two planes.
    
          !!
    
          !! The improper dihedral angle is the angle between two planes made by
    
          !! three vectors.
    
          !! i.e. ( vector1 x vector2 ) . ( vector2 x vector3 )
    
          !! alt. angle between plane vector1vector2 and vector2vector3
    
          implicit none
    
          real(real32), dimension(3), intent(in) :: vector1, vector2, vector3
    
          real(real32) :: angle
    
          angle = get_angle( &
    
               cross(vector1, vector2), &
    
               cross(vector2, vector3) &
    
      278094
          )
    
          !! map angle back into the range [0, pi]
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 278094 times.

      278094
          if(angle .gt. pi) angle = 2._real32 * pi - angle
    
      278094
        end function get_improper_dihedral_angle_from_vectors
    
      !###############################################################################
    
      !###############################################################################
    
      119170665
        pure function get_improper_dihedral_angle_from_points( &
    
             point1, point2, point3, point4 &
    
        ) result(angle)
    
          !! Return the improper dihedral angle between two planes.
    
          !!
    
          !! The dihedral angle is the angle between the plane defined by four points.
    
          !! i.e. ( point2 - point1 ) x ( point3 - point1 ) . 
    
          !! ( point4 - point2 ) x ( point3 - point1 )
    
          !! alt. angle between plane point1point2point3 and point1point3point4
    
          implicit none
    
          real(real32), dimension(3), intent(in) :: point1, point2, point3, point4
    
          real(real32) :: angle
    
          angle = get_angle( &
    
               cross(point2 - point1, point3 - point1), &
    
               cross(point3 - point1, point4 - point1) &
    
        8/8✓ Branch 0 taken 357511995 times.
✓ Branch 1 taken 119170665 times.
✓ Branch 2 taken 357511995 times.
✓ Branch 3 taken 119170665 times.
✓ Branch 5 taken 357511995 times.
✓ Branch 6 taken 119170665 times.
✓ Branch 7 taken 357511995 times.
✓ Branch 8 taken 119170665 times.

      1549218645
          )
    
      119170665
        end function get_improper_dihedral_angle_from_points
    
      !###############################################################################
    
      !###############################################################################
    
        1/2✓ Branch 0 taken 54267 times.
✗ Branch 1 not taken.

      54267
        pure function inverse_3x3(mat) result(output)
    
          implicit none
    
          real(real32) :: det
    
          real(real32), dimension(3,3) :: output
    
          real(real32), dimension(3,3), intent(in) :: mat
    
          det = &
    
               mat(1,1) * mat(2,2) * mat(3,3) - mat(1,1) * mat(2,3) * mat(3,2) - &
    
               mat(1,2) * mat(2,1) * mat(3,3) + mat(1,2) * mat(2,3) * mat(3,1) + &
    
      54267
               mat(1,3) * mat(2,1) * mat(3,2) - mat(1,3) * mat(2,2) * mat(3,1)
    
      54267
          output(1,1) = +1._real32 / det * (mat(2,2) * mat(3,3) - mat(2,3) * mat(3,2))
    
      54267
          output(2,1) = -1._real32 / det * (mat(2,1) * mat(3,3) - mat(2,3) * mat(3,1))
    
      54267
          output(3,1) = +1._real32 / det * (mat(2,1) * mat(3,2) - mat(2,2) * mat(3,1))
    
      54267
          output(1,2) = -1._real32 / det * (mat(1,2) * mat(3,3) - mat(1,3) * mat(3,2))
    
      54267
          output(2,2) = +1._real32 / det * (mat(1,1) * mat(3,3) - mat(1,3) * mat(3,1))
    
      54267
          output(3,2) = -1._real32 / det * (mat(1,1) * mat(3,2) - mat(1,2) * mat(3,1))
    
      54267
          output(1,3) = +1._real32 / det * (mat(1,2) * mat(2,3) - mat(1,3) * mat(2,2))
    
      54267
          output(2,3) = -1._real32 / det * (mat(1,1) * mat(2,3) - mat(1,3) * mat(2,1))
    
      54267
          output(3,3) = +1._real32 / det * (mat(1,1) * mat(2,2) - mat(1,2) * mat(2,1))
    
      54267
        end function inverse_3x3
    
      !###############################################################################
    
      end module raffle__misc_linalg

Line	Branch	Exec	Source
1			module raffle__misc_linalg
2			!! Module contains various linear algebra functions and subroutines.
3			use raffle__constants, only: real32, pi
4			implicit none
5
6
7			private
8
9			public :: modu, cross
10			public :: get_distance, get_angle, get_dihedral_angle
11			public :: get_improper_dihedral_angle
12			public :: inverse_3x3
13
14
15			interface get_angle
16			procedure get_angle_from_points, get_angle_from_vectors
17			end interface get_angle
18
19			interface get_dihedral_angle
20			procedure get_dihedral_angle_from_points, get_dihedral_angle_from_vectors
21			end interface get_dihedral_angle
22
23			interface get_improper_dihedral_angle
24			procedure get_improper_dihedral_angle_from_points, &
25			get_improper_dihedral_angle_from_vectors
26			end interface get_improper_dihedral_angle
27
28			contains
29
30			!###############################################################################
31	1/2 ✓ Branch 0 taken 422656602 times. ✗ Branch 1 not taken.	422656602	pure function modu(vector)
32			!! Return the magnitude of a vector of any size.
33			implicit none
34
35			! Arguments
36			real(real32),dimension(:), intent(in) :: vector
37			!! Input vector.
38			real(real32) :: modu
39			!! Output magnitude.
40
41	2/2 ✓ Branch 0 taken 1267969806 times. ✓ Branch 1 taken 422656602 times.	1690626408	modu = abs( sqrt( sum(vector(:)**2) ) )
42		422656602	end function modu
43			!###############################################################################
44
45
46			!###############################################################################
47	1/2 ✓ Branch 0 taken 239223063 times. ✗ Branch 1 not taken.	239223063	pure function cross(a,b)
48			!! Return the cross product of two vectors.
49			implicit none
50
51			! Arguments
52			real(real32), dimension(3), intent(in) :: a, b
53			!! Input vectors.
54			real(real32), dimension(3) :: cross
55			!! Output cross product.
56
57		239223063	cross(1) = a(2) * b(3) - a(3) * b(2)
58		239223063	cross(2) = a(3) * b(1) - a(1) * b(3)
59		239223063	cross(3) = a(1) * b(2) - a(2) * b(1)
60
61			end function cross
62			!###############################################################################
63
64
65			!###############################################################################
66		1	pure function get_distance(point1, point2) result(distance)
67			!! Return the distance between two points.
68			implicit none
69
70			! Arguments
71			real(real32), dimension(3), intent(in) :: point1, point2
72			!! Input points.
73			real(real32) :: distance
74			!! Output distance.
75
76	2/2 ✓ Branch 0 taken 3 times. ✓ Branch 1 taken 1 times.	4	distance = modu( point1 - point2 )
77
78		1	return
79			end function get_distance
80			!###############################################################################
81
82
83			!###############################################################################
84		119453303	pure function get_angle_from_vectors(vector1, vector2) result(angle)
85			!! Return the angle between two vectors.
86			implicit none
87
88			! Arguments
89			real(real32), dimension(3), intent(in) :: vector1, vector2
90			!! Input vectors.
91			real(real32) :: angle
92			!! Output angle.
93
94			angle = dot_product(vector1,vector2) / &
95	2/2 ✓ Branch 0 taken 358359909 times. ✓ Branch 1 taken 119453303 times.	477813212	( modu(vector1) * modu(vector2) )
96	2/2 ✓ Branch 0 taken 1785068 times. ✓ Branch 1 taken 117668235 times.	119453303	if(angle .ge. 1._real32)then
97		1785068	angle = 0._real32
98	2/2 ✓ Branch 0 taken 1750023 times. ✓ Branch 1 taken 115918212 times.	117668235	elseif(angle .le. -1._real32)then
99		1750023	angle = pi
100			else
101		115918212	angle = acos(angle)
102			end if
103		119453303	end function get_angle_from_vectors
104			!###############################################################################
105
106
107			!###############################################################################
108		3958532	pure function get_angle_from_points(point1, point2, point3) result(angle)
109			!! Return the angle formed by three points.
110			!!
111			!! The angle is formed by the path point1 -> point2 -> point3.
112			implicit none
113
114			! Arguments
115			real(real32), dimension(3), intent(in) :: point1, point2, point3
116			!! Input points.
117			real(real32) :: angle
118			!! Output angle.
119
120			angle = dot_product( point1 - point2, point3 - point2 ) / &
121	6/6 ✓ Branch 0 taken 11875596 times. ✓ Branch 1 taken 3958532 times. ✓ Branch 2 taken 11875596 times. ✓ Branch 3 taken 3958532 times. ✓ Branch 4 taken 11875596 times. ✓ Branch 5 taken 3958532 times.	39585320	( modu( point1 - point2 ) * modu( point3 - point2 ) )
122	2/2 ✓ Branch 0 taken 12 times. ✓ Branch 1 taken 3958520 times.	3958532	if(angle .ge. 1._real32)then
123		12	angle = 0._real32
124	2/2 ✓ Branch 0 taken 2689 times. ✓ Branch 1 taken 3955831 times.	3958520	elseif(angle .le. -1._real32)then
125		2689	angle = pi
126			else
127		3955831	angle = acos(angle)
128			end if
129		3958532	end function get_angle_from_points
130			!###############################################################################
131
132
133			!###############################################################################
134		1	pure function get_dihedral_angle_from_vectors( &
135			vector1, vector2, vector3) result(angle)
136			!! Return the dihedral angle between two planes.
137			!!
138			!! The dihedral angle is the angle between the plane defined by the cross
139			!! product of two vectors and a third vector.
140			!! i.e. ( vector1 x vector2 ) . vector3
141			implicit none
142
143			! Arguments
144			real(real32), dimension(3), intent(in) :: vector1, vector2, vector3
145			!! Input vectors.
146			real(real32) :: angle
147			!! Output angle.
148
149		1	angle = get_angle(cross(vector1, vector2), vector3)
150
151		1	end function get_dihedral_angle_from_vectors
152			!###############################################################################
153
154
155			!###############################################################################
156		1	pure function get_dihedral_angle_from_points( &
157			point1, point2, point3, point4 &
158			) result(angle)
159			!! Return the dihedral angle between two planes.
160			!!
161			!! The dihedral angle is the angle between the plane defined by four points.
162			!! i.e. ( point2 - point1 ) x ( point3 - point2 ) . ( point4 - point2 )
163			!! alt. angle between plane point1point2point3 and vector point2point4
164			implicit none
165			real(real32), dimension(3), intent(in) :: point1, point2, point3, point4
166			real(real32) :: angle
167
168	6/6 ✓ Branch 0 taken 3 times. ✓ Branch 1 taken 1 times. ✓ Branch 2 taken 3 times. ✓ Branch 3 taken 1 times. ✓ Branch 5 taken 3 times. ✓ Branch 6 taken 1 times.	10	angle = get_angle(cross(point2 - point1, point3 - point2), point4 - point2)
169
170		1	end function get_dihedral_angle_from_points
171			!###############################################################################
172
173
174			!###############################################################################
175		278094	pure function get_improper_dihedral_angle_from_vectors( &
176			vector1, vector2, vector3 &
177			) result(angle)
178			!! Return the improper dihedral angle between two planes.
179			!!
180			!! The improper dihedral angle is the angle between two planes made by
181			!! three vectors.
182			!! i.e. ( vector1 x vector2 ) . ( vector2 x vector3 )
183			!! alt. angle between plane vector1vector2 and vector2vector3
184			implicit none
185			real(real32), dimension(3), intent(in) :: vector1, vector2, vector3
186			real(real32) :: angle
187
188			angle = get_angle( &
189			cross(vector1, vector2), &
190			cross(vector2, vector3) &
191		278094	)
192			!! map angle back into the range [0, pi]
193	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 278094 times.	278094	if(angle .gt. pi) angle = 2._real32 * pi - angle
194
195
196		278094	end function get_improper_dihedral_angle_from_vectors
197			!###############################################################################
198
199
200			!###############################################################################
201		119170665	pure function get_improper_dihedral_angle_from_points( &
202			point1, point2, point3, point4 &
203			) result(angle)
204			!! Return the improper dihedral angle between two planes.
205			!!
206			!! The dihedral angle is the angle between the plane defined by four points.
207			!! i.e. ( point2 - point1 ) x ( point3 - point1 ) .
208			!! ( point4 - point2 ) x ( point3 - point1 )
209			!! alt. angle between plane point1point2point3 and point1point3point4
210			implicit none
211			real(real32), dimension(3), intent(in) :: point1, point2, point3, point4
212			real(real32) :: angle
213
214			angle = get_angle( &
215			cross(point2 - point1, point3 - point1), &
216			cross(point3 - point1, point4 - point1) &
217	8/8 ✓ Branch 0 taken 357511995 times. ✓ Branch 1 taken 119170665 times. ✓ Branch 2 taken 357511995 times. ✓ Branch 3 taken 119170665 times. ✓ Branch 5 taken 357511995 times. ✓ Branch 6 taken 119170665 times. ✓ Branch 7 taken 357511995 times. ✓ Branch 8 taken 119170665 times.	1549218645	)
218
219		119170665	end function get_improper_dihedral_angle_from_points
220			!###############################################################################
221
222
223			!###############################################################################
224	1/2 ✓ Branch 0 taken 54267 times. ✗ Branch 1 not taken.	54267	pure function inverse_3x3(mat) result(output)
225			implicit none
226			real(real32) :: det
227			real(real32), dimension(3,3) :: output
228			real(real32), dimension(3,3), intent(in) :: mat
229
230			det = &
231			mat(1,1) * mat(2,2) * mat(3,3) - mat(1,1) * mat(2,3) * mat(3,2) - &
232			mat(1,2) * mat(2,1) * mat(3,3) + mat(1,2) * mat(2,3) * mat(3,1) + &
233		54267	mat(1,3) * mat(2,1) * mat(3,2) - mat(1,3) * mat(2,2) * mat(3,1)
234
235		54267	output(1,1) = +1._real32 / det * (mat(2,2) * mat(3,3) - mat(2,3) * mat(3,2))
236		54267	output(2,1) = -1._real32 / det * (mat(2,1) * mat(3,3) - mat(2,3) * mat(3,1))
237		54267	output(3,1) = +1._real32 / det * (mat(2,1) * mat(3,2) - mat(2,2) * mat(3,1))
238		54267	output(1,2) = -1._real32 / det * (mat(1,2) * mat(3,3) - mat(1,3) * mat(3,2))
239		54267	output(2,2) = +1._real32 / det * (mat(1,1) * mat(3,3) - mat(1,3) * mat(3,1))
240		54267	output(3,2) = -1._real32 / det * (mat(1,1) * mat(3,2) - mat(1,2) * mat(3,1))
241		54267	output(1,3) = +1._real32 / det * (mat(1,2) * mat(2,3) - mat(1,3) * mat(2,2))
242		54267	output(2,3) = -1._real32 / det * (mat(1,1) * mat(2,3) - mat(1,3) * mat(2,1))
243		54267	output(3,3) = +1._real32 / det * (mat(1,1) * mat(2,2) - mat(1,2) * mat(2,1))
244
245		54267	end function inverse_3x3
246			!###############################################################################
247
248			end module raffle__misc_linalg
249