Cookbook

Short answers to "how do I ...?". Each recipe is a whole program with its real output; the reference has the details.

• Building vectors
  • A vector from its components
  • Every component set to one number
• Arithmetic
  • Sum, difference, negation
  • A vector times a number
  • A number added to every component
  • The product and the quotient of two vectors, component by component
• Products, lengths, angles
  • Dot, cross and triple products
  • The length of a vector
  • A unit vector
  • The angle between two vectors
  • The components parallel and orthogonal to a vector
  • A velocity reflected by a wall
• Faces
  • The normal and the area of a triangle
  • The normal and the area of a quadrilateral
• Points, lines, planes
  • The distance of a point from a line
  • The distance of a point from a plane
  • Are three points on one line?
  • Are four points on one circle?
• Rotations and mirrors
  • A rotation about an axis
  • A rotation matrix applied many times
  • The mirror image across a plane
  • A mirror matrix applied many times
  • Any matrix applied to a vector
• Comparisons
  • Comparing two vectors
  • Comparing the length of a vector with a number
• Arrays and kinds
  • Whole arrays of vectors at once
  • Vectors of each precision
  • True quadruple precision
  • A method or a free function
• Input and output
  • Printing a vector
  • A binary file of vectors
  • A text file of vectors
  • The n-th vector of a stream file
  • A direct access file of vectors
  • An I/O error caught
• On the GPU
  • The device-callable routines in an OpenACC loop

Building vectors

A vector from its components

program from_components
!< Cookbook: a vector from its components.
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: a, b, c, d

a%x = 1._R8P ; a%y = 2._R8P ; a%z = 3._R8P   ! component by component
b = 1 * ex + 2 * ey + 3 * ez                  ! from the versors
d = vector(1._R8P, 2._R8P, 3._R8P)            ! the structure constructor: x, y, z
call a%printf(prefix='a = ')
call b%printf(prefix='b = ')
call d%printf(prefix='d = ')
call c%printf(prefix='c = ')                  ! never set: (0, 0, 0)
endprogram from_components

$ from_components
a = +1.0, +2.0, +3.0
b = +1.0, +2.0, +3.0
d = +1.0, +2.0, +3.0
c = +0.0, +0.0, +0.0

Three spellings of the same vector: component by component, from the versors, with the structure constructor vector(x, y, z). A vector never set is the origin.

Every component set to one number

program from_number
!< Cookbook: every component set to one number, of any integer or real kind.
use penf, only : I1P, I2P, I4P, I8P, R4P, R8P
use vecfor, only : vector
implicit none
type(vector) :: v

v = 2.5_R8P  ; call v%printf(prefix='real(R8P)    ')
v = 0.1_R4P  ; call v%printf(prefix='real(R4P)    ')
v = 7_I8P    ; call v%printf(prefix='integer(I8P) ')
v = -3_I4P   ; call v%printf(prefix='integer(I4P) ')
v = 4_I2P    ; call v%printf(prefix='integer(I2P) ')
v = 1_I1P    ; call v%printf(prefix='integer(I1P) ')
endprogram from_number

$ from_number
real(R8P)    +2.5, +2.5, +2.5
real(R4P)    +0.10000000149011612, +0.10000000149011612, +0.10000000149011612
integer(I8P) +7.0, +7.0, +7.0
integer(I4P) -3.0, -3.0, -3.0
integer(I2P) +4.0, +4.0, +4.0
integer(I1P) +1.0, +1.0, +1.0

Any integer or real kind. The number is converted to the kind of the vector, with its own precision: 0.1_R4P is the single precision number nearest to 0.1, not the double precision one.

Arithmetic

Sum, difference, negation

program sum_difference
!< Cookbook: sum, difference, negation.
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: a, b, v

a = ex + 2 * ey
b = 3 * ex - ez
v = a + b  ; call v%printf(prefix='a + b = ')
v = a - b  ; call v%printf(prefix='a - b = ')
v = -a     ; call v%printf(prefix='-a    = ')
v = +a     ; call v%printf(prefix='+a    = ')
endprogram sum_difference

$ sum_difference
a + b = +4.0, +2.0, -1.0
a - b = -2.0, +2.0, +1.0
-a    = -1.0, -2.0, -0.0
+a    = +1.0, +2.0, +0.0

-a of a zero component is -0.0, the IEEE negative zero, equal to +0.0.

A vector times a number

program scaling
!< Cookbook: a vector times or divided by a number, of any integer or real kind, on either side of *.
use penf, only : I8P, R4P, R8P
use vecfor, only : ex, ey, vector
implicit none
type(vector) :: a, v

a = 2 * ex + 4 * ey
v = 3 * a        ; call v%printf(prefix='3 * a      = ')
v = a * 0.5_R8P  ; call v%printf(prefix='a * 0.5    = ')
v = a / 4        ; call v%printf(prefix='a / 4      = ')
v = 2_I8P * a    ; call v%printf(prefix='2_I8P * a  = ')
v = a * 0.1_R4P  ; call v%printf(prefix='a * 0.1_R4P = ')
v = a * 0.1_R8P  ; call v%printf(prefix='a * 0.1_R8P = ')
endprogram scaling

$ scaling
3 * a      = +6.0, +12.0, +0.0
a * 0.5    = +1.0, +2.0, +0.0
a / 4      = +0.5, +1.0, +0.0
2_I8P * a  = +4.0, +8.0, +0.0
a * 0.1_R4P = +0.20000000298023224, +0.4000000059604645, +0.0
a * 0.1_R8P = +0.2, +0.4, +0.0

* takes the number on either side, / on the right; any integer or real kind, converted to the kind of the vector before the product. 0.1_R4P carries its single precision error into a double precision vector: write the literals with the kind of the vector.

A number added to every component

program plus_number
!< Cookbook: a number added to or subtracted from every component.
use penf, only : R8P
use vecfor, only : ex, ey, vector
implicit none
type(vector) :: a, v

a = ex + 2 * ey
v = a + 1         ; call v%printf(prefix='a + 1   = ')
v = 1 + a         ; call v%printf(prefix='1 + a   = ')
v = a - 0.5_R8P   ; call v%printf(prefix='a - 0.5 = ')
v = 0.5_R8P - a   ; call v%printf(prefix='0.5 - a = ')
endprogram plus_number

$ plus_number
a + 1   = +2.0, +3.0, +1.0
1 + a   = +2.0, +3.0, +1.0
a - 0.5 = +0.5, +1.5, -0.5
0.5 - a = -0.5, -1.5, +0.5

a + 1 adds 1 to x, y and z alike. To move a point along an axis, add a vector: a + ex.

The product and the quotient of two vectors, component by component

program componentwise
!< Cookbook: the product and the quotient of two vectors, component by component.
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: a, s, v

a = 2 * ex + 3 * ey + 4 * ez
s = 2 * ex + ey + ez                ! a scale factor for each axis
v = a * s ; call v%printf(prefix='a * s = ')
v = a / s ; call v%printf(prefix='a / s = ')
endprogram componentwise

$ componentwise
a * s = +4.0, +3.0, +4.0
a / s = +1.0, +3.0, +4.0

a * s is $(a_x{s}_x,a_y{s}_y,a_z{s}_z)$: a different scale on each axis. It is not the dot product, .dot., nor the cross product, .cross..

Products, lengths, angles

Dot, cross and triple products

program products
!< Cookbook: dot product, cross product, triple product.
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: a, b, c, v

a = ex + 2 * ey
b = 3 * ex + ey
c = ez
print '(A,F6.2)', 'a . b         = ', a .dot. b
v = a .cross. b ; call v%printf(prefix='a x b         = ')
v = b .cross. a ; call v%printf(prefix='b x a         = ')
v = a .cross. a ; call v%printf(prefix='a x a         = ')
print '(A,F6.2)', 'a . (b x c)   = ', a .dot. (b .cross. c)
endprogram products

$ products
a . b         =   5.00
a x b         = +0.0, +0.0, -5.0
b x a         = +0.0, +0.0, +5.0
a x a         = +0.0, +0.0, +0.0
a . (b x c)   =  -5.00

The cross product is anti-commutative, zero for parallel vectors; a .dot. (b .cross. c) is the signed volume of the parallelepiped of the three vectors. Parenthesise: user-defined binary operators such as .dot. bind more loosely than any intrinsic one, so a .dot. b + c is a .dot. (b + c).

The length of a vector

program lengths
!< Cookbook: the length and the squared length of a vector.
use vecfor, only : ex, ey, ez, normL2, sq_norm, vector
implicit none
type(vector) :: a

a = 2 * ex + 3 * ey + 6 * ez
print '(A,F6.2)', 'a%normL2()  = ', a%normL2()
print '(A,F6.2)', 'normL2(a)   = ', normL2(a)
print '(A,F6.2)', 'a%sq_norm() = ', a%sq_norm()
print '(A,F6.2)', 'sq_norm(a)  = ', sq_norm(a)
endprogram lengths

$ lengths
a%normL2()  =   7.00
normL2(a)   =   7.00
a%sq_norm() =  49.00
sq_norm(a)  =  49.00

sq_norm is the squared length, without a square root: enough to compare lengths.

A unit vector

program unit_vector
!< Cookbook: a unit vector, as a new vector or in place; the zero vector.
use vecfor, only : ex, ey, ez, normalized, vector
implicit none
type(vector) :: a, u, z

a = 2 * ex + 3 * ey + 6 * ez
u = a%normalized()   ; call u%printf(prefix='a%normalized()  = ')
u = normalized(a)    ; call u%printf(prefix='normalized(a)   = ')
call a%printf(prefix='a, unchanged    = ')
call a%normalize     ; call a%printf(prefix='after normalize = ')
call z%normalize     ; call z%printf(prefix='zero, normalized = ')
endprogram unit_vector

$ unit_vector
a%normalized()  = +0.2857142857142857, +0.42857142857142855, +0.8571428571428571
normalized(a)   = +0.2857142857142857, +0.42857142857142855, +0.8571428571428571
a, unchanged    = +2.0, +3.0, +6.0
after normalize = +0.2857142857142857, +0.42857142857142855, +0.8571428571428571
zero, normalized = +0.0, +0.0, +0.0

normalized returns a copy, normalize changes the vector. The zero vector stays zero: a length below tiny(1._R8P) is replaced by that length plus tiny, so the division is never by zero.

The angle between two vectors

program angles
!< Cookbook: the angle between two vectors, in radians and degrees.
use penf, only : R8P
use vecfor, only : angle, ex, ey, vector
implicit none
type(vector) :: a, b
real(R8P)    :: deg

deg = 180._R8P / acos(-1._R8P)
a = ex
b = ex + ey
print '(A,F8.4)', 'a%angle(b)         [rad] ', a%angle(b)
print '(A,F8.4)', 'angle(a, b)        [deg] ', angle(a, b) * deg
print '(A,F8.4)', 'a and 3a           [deg] ', angle(a, 3 * a) * deg
print '(A,F8.4)', 'a and -a           [deg] ', angle(a, -a) * deg
print '(A,F8.4)', 'a and ey           [deg] ', angle(a, ey) * deg
endprogram angles

$ angles
a%angle(b)         [rad]   0.7854
angle(a, b)        [deg]  45.0000
a and 3a           [deg]   0.0000
a and -a           [deg] 180.0000
a and ey           [deg]  90.0000

In radians, between 0 and π, from atan2: accurate near 0 and π too. With the zero vector the angle is 0.

The components parallel and orthogonal to a vector

program paral_ortho
!< Cookbook: the components of a vector parallel and orthogonal to another one.
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: a, d, z, v

a = 2 * ex + 3 * ey + ez
d = 5 * ex                         ! only the direction of d matters
v = a .paral. d ; call v%printf(prefix='a .paral. d = ')
v = a .ortho. d ; call v%printf(prefix='a .ortho. d = ')
v = (a .paral. d) + (a .ortho. d) ; call v%printf(prefix='their sum   = ')
v = a .paral. z ; call v%printf(prefix='a .paral. 0 = ')
endprogram paral_ortho

$ paral_ortho
a .paral. d = +2.0, +0.0, +0.0
a .ortho. d = +0.0, +3.0, +1.0
their sum   = +2.0, +3.0, +1.0
a .paral. 0 = NaN, NaN, NaN

Only the direction of the right operand matters. The zero vector has none: .paral. and .ortho. give NaN.

A velocity reflected by a wall

program reflect
!< Cookbook: a velocity reflected by a wall, by components or by a mirror.
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: v, n, r1, r2

v = 2 * ex - 3 * ey + ez
n = ey                                   ! the normal of the wall, any length
r1 = v - 2 * (v .paral. n)
r2 = v
call r2%mirror(normal=n)
call r1%printf(prefix='v - 2 (v .paral. n) = ')
call r2%printf(prefix='v mirrored          = ')
endprogram reflect

$ reflect
v - 2 (v .paral. n) = +2.0, +3.0, +1.0
v mirrored          = +2.0, +3.0, +1.0

The same reflection, by components or by a mirror across the plane of the wall (see tutorial, chapter 6).

Faces

The normal and the area of a triangle

program triangle
!< Cookbook: the normal and the area of a triangle.
use vecfor, only : ex, ey, ez, face_normal3, vector
implicit none
type(vector) :: a, b, c, n

a = 0 * ex ; b = 4 * ex ; c = 3 * ey                         ! counter-clockwise seen from +z
n = face_normal3(pt1=a, pt2=b, pt3=c)          ; call n%printf(prefix='normal        = ')
print '(A,F6.2)', 'area          = ', n%normL2()
n = face_normal3(pt1=a, pt2=b, pt3=c, norm='y') ; call n%printf(prefix='unit normal   = ')
n = face_normal3(pt1=a, pt2=c, pt3=b, norm='y') ; call n%printf(prefix='b and c swapped = ')
n = n%face_normal3(pt1=a, pt2=b, pt3=c)          ; call n%printf(prefix='type-bound    = ')
endprogram triangle

$ triangle
normal        = +0.0, +0.0, +6.0
area          =   6.00
unit normal   = +0.0, +0.0, +1.0
b and c swapped = +0.0, +0.0, -1.0
type-bound    = +0.0, +0.0, +6.0

The length of the normal is the area; norm='y' makes it a unit vector; the order of the points sets its side. n%face_normal3(...) does not use n: the method is a free function bound to the type.

The normal and the area of a quadrilateral

program quadrilateral
!< Cookbook: the normal and the area of a quadrilateral.
use vecfor, only : ex, ey, ez, face_normal4, vector
implicit none
type(vector) :: q(4), n

q = [0 * ex, 4 * ex, 4 * ex + 3 * ey, 3 * ey]          ! a 4 x 3 rectangle, counter-clockwise seen from +z
n = face_normal4(pt1=q(1), pt2=q(2), pt3=q(3), pt4=q(4))           ; call n%printf(prefix='normal      = ')
print '(A,F6.2)', 'area        = ', n%normL2()
n = face_normal4(pt1=q(1), pt2=q(2), pt3=q(3), pt4=q(4), norm='y') ; call n%printf(prefix='unit normal = ')
q(3) = q(3) + ez                                       ! one corner lifted: no longer flat
n = face_normal4(pt1=q(1), pt2=q(2), pt3=q(3), pt4=q(4), norm='y') ; call n%printf(prefix='warped      = ')
endprogram quadrilateral

$ quadrilateral
normal      = +0.0, -0.0, +12.0
area        =  12.00
unit normal = +0.0, -0.0, +1.0
warped      = -0.12237255246720063, -0.16316340328960083, +0.9789804197376051

From the cross product of the diagonals, $\frac{1}{2}(p_3-p_1)\times (p_4-p_2)$: the exact area of a flat quadrilateral. For a warped one it is the vector area of its four edges, the same for every surface they bound: the normal to use for a flux through the face.

Points, lines, planes

The distance of a point from a line

program to_line
!< Cookbook: the distance of a point from a line.
use penf, only : R8P
use vecfor, only : distance_to_line, ex, ey, ez, vector
implicit none
type(vector) :: a, b, p

a = 0 * ex
b = 2 * ex                                     ! the x axis
p = 5 * ex + 3 * ey + 4 * ez
print '(A,F6.2)', 'p%distance_to_line(...) = ', p%distance_to_line(pt1=a, pt2=b)
print '(A,F6.2)', 'distance_to_line(p,...) = ', distance_to_line(p, pt1=a, pt2=b)
endprogram to_line

$ to_line
p%distance_to_line(...) =   5.00
distance_to_line(p,...) =   5.00

The line through the two points, infinite: the point (5, 3, 4) is 5 from the x axis, beyond the segment of the two points.

The distance of a point from a plane

program to_plane
!< Cookbook: the signed distance of a point from a plane, and the vector from the plane to the point.
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: a, b, c, p, q, d

a = ez ; b = ex + ez ; c = ey + ez                       ! the plane z = 1, normal +z
p = 2 * ex + 3 * ey + 4 * ez
print '(A,F6.2)', 'above:            ', p%distance_to_plane(pt1=a, pt2=b, pt3=c)
q = 2 * ex - 3 * ez
print '(A,F6.2)', 'below:            ', q%distance_to_plane(pt1=a, pt2=b, pt3=c)
print '(A,F6.2)', 'points swapped:   ', p%distance_to_plane(pt1=a, pt2=c, pt3=b)
d = p%distance_vectorial_to_plane(pt1=a, pt2=b, pt3=c) ; call d%printf(prefix='vectorial:        ')
d = p%projection_onto_plane(pt1=a, pt2=b, pt3=c)       ; call d%printf(prefix='projection:       ')
endprogram to_plane

$ to_plane
above:              3.00
below:             -4.00
points swapped:    -3.00
vectorial:        +0.0, +0.0, +3.0
projection:       +2.0, +3.0, +1.0

Signed: positive on the side of face_normal3(pt1, pt2, pt3), which swapping two points reverses. The vectorial distance goes from the plane to the point; the projection is the foot of the point on the plane.

Are three points on one line?

program collinear
!< Cookbook: are three points on one line? Exact by default, within a tolerance if asked.
use penf, only : R8P
use vecfor, only : ex, ey, is_collinear, vector
implicit none
type(vector) :: a, b, p

a = ex ; b = ey
p = (a + b) / 2
print '(A,L1)', 'midpoint:                   ', p%is_collinear(pt1=a, pt2=b)
p = 0.1_R8P * a + 0.9_R8P * b
print '(A,L1)', '0.1 a + 0.9 b:              ', p%is_collinear(pt1=a, pt2=b)
print '(A,L1)', '0.1 a + 0.9 b, within 1e-12: ', is_collinear(p, pt1=a, pt2=b, tolerance=1.e-12_R8P)
p = 2 * a                                                  ! off the line
print '(A,L1)', '2 a, within 1e-12:          ', p%is_collinear(pt1=a, pt2=b, tolerance=1.e-12_R8P)
endprogram collinear

$ collinear
midpoint:                   T
0.1 a + 0.9 b:              F
0.1 a + 0.9 b, within 1e-12: T
2 a, within 1e-12:          F

Exact by default, so computed points need a tolerance, a length in the units of the geometry.

Are four points on one circle?

program concyclic_points
!< Cookbook: are four points on one circle? The order matters.
use penf, only : R8P
use vecfor, only : ex, ey, ez, is_concyclic, vector
implicit none
type(vector) :: a, b, c, d

a = ex ; b = ey ; c = -1 * ex ; d = -1 * ey               ! the unit circle, counter-clockwise
print '(A,L1)', 'a, b, c, d in order:         ', a%is_concyclic(pt1=b, pt2=c, pt3=d)
print '(A,L1)', 'a, c, b, d out of order:     ', a%is_concyclic(pt1=c, pt2=b, pt3=d)
d = -1 * ey + 0.1_R8P * ez                                 ! lifted off the plane
print '(A,L1)', 'd lifted by 0.1:             ', is_concyclic(a, pt1=b, pt2=c, pt3=d)
d = cos(4._R8P) * ex + sin(4._R8P) * ey                    ! on the circle, up to the rounding of cos and sin
print '(A,L1)', 'a, b, c, (cos 4, sin 4):     ', a%is_concyclic(pt1=b, pt2=c, pt3=d)
print '(A,L1)', 'the same, within 1e-12:      ', a%is_concyclic(pt1=b, pt2=c, pt3=d, tolerance=1.e-12_R8P)
endprogram concyclic_points

$ concyclic_points
a, b, c, d in order:         T
a, c, b, d out of order:     F
d lifted by 0.1:             F
a, b, c, (cos 4, sin 4):     F
the same, within 1e-12:      T

The four points must be given in their order around the circle. Without tolerance the test succeeds or fails on the rounding of the last digit: for computed points pass a tolerance.

Rotations and mirrors

A rotation about an axis

program rotate_axis
!< Cookbook: a rotation about an axis through the origin; the angle in radians, counter-clockwise about the axis.
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: v
real(R8P)    :: pi

pi = acos(-1._R8P)
v = ex
call v%rotate(axis=ez, angle=pi / 2)
print '(A,3F8.4)', 'ex, 90 deg about z      ', v%x, v%y, v%z
v = ex
call v%rotate(axis=-3 * ez, angle=pi / 2)                  ! the axis need not be a unit vector
print '(A,3F8.4)', 'ex, 90 deg about -z     ', v%x, v%y, v%z
v = ex + ez
call v%rotate(axis=ez, angle=pi)
print '(A,3F8.4)', 'ex + ez, 180 deg about z', v%x, v%y, v%z
endprogram rotate_axis

$ rotate_axis
ex, 90 deg about z        0.0000  1.0000  0.0000
ex, 90 deg about -z       0.0000 -1.0000  0.0000
ex + ez, 180 deg about z -1.0000  0.0000  1.0000

Counter-clockwise seen from the tip of the axis, the angle in radians; the axis passes through the origin and need not be a unit vector.

A rotation matrix applied many times

program rotate_matrix
!< Cookbook: a rotation matrix, computed once and applied many times; two rotations composed.
use penf, only : R8P
use vecfor, only : ex, ey, ez, rotation_matrix, vector
implicit none
type(vector) :: v(3), w
real(R8P)    :: r(3,3), rz(3,3), rx(3,3), pi
integer      :: i

pi = acos(-1._R8P)
r = rotation_matrix(axis=ex + ey + ez, angle=2 * pi / 3)
v = [ex, ey, ez]
do i = 1, 3
   w = r .matrix. v(i)                                    ! a new vector
   print '(A,I0,A,3F8.4)', 'v', i, ' rotated  ', w%x, w%y, w%z
   call v(i)%rotate(matrix=r)                             ! in place: the same
enddo
print '(A,3F8.4)', 'v1 in place ', v(1)%x, v(1)%y, v(1)%z
rz = rotation_matrix(axis=ez, angle=pi / 2)
rx = rotation_matrix(axis=ex, angle=pi / 2)
w = matmul(rx, rz) .matrix. ex                            ! first about z, then about x
print '(A,3F8.4)', 'ex, z then x', w%x, w%y, w%z
endprogram rotate_matrix

$ rotate_matrix
v1 rotated    0.0000  1.0000  0.0000
v2 rotated    0.0000  0.0000  1.0000
v3 rotated    1.0000  0.0000  0.0000
v1 in place   0.0000  1.0000  0.0000
ex, z then x  0.0000  0.0000  1.0000

r .matrix. v returns a new vector, call v%rotate(matrix=r) changes v. Matrices compose with matmul, the rightmost first.

The mirror image across a plane

program mirror_plane
!< Cookbook: the mirror image of a point across a plane, through the origin or not.
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: p, n, o, m

p = 3 * ex + ey + 2 * ez
n = ex + ey                                    ! the plane x + y = 0
m = p
call m%mirror(normal=n)
print '(A,3F8.4)', 'across x + y = 0  ', m%x, m%y, m%z
o = 2 * ez                                     ! a point of the plane z = 2
m = p - o
call m%mirror(normal=ez)
m = m + o
print '(A,3F8.4)', 'across z = 2      ', m%x, m%y, m%z
endprogram mirror_plane

$ mirror_plane
across x + y = 0   -1.0000 -3.0000  2.0000
across z = 2        3.0000  1.0000  2.0000

mirror reflects across the plane through the origin orthogonal to normal. For a plane through another point, move the point there first, and back after.

A mirror matrix applied many times

program mirror_matrix_reuse
!< Cookbook: a mirror matrix, computed once and applied many times.
use penf, only : R8P
use vecfor, only : ex, ey, ez, mirror_matrix, vector
implicit none
type(vector) :: v(2)
real(R8P)    :: m(3,3)
integer      :: i

m = mirror_matrix(normal=ex + ey + ez)
print '(A,3F8.4)', 'M  ', m(1,:)
print '(A,3F8.4)', '   ', m(2,:)
print '(A,3F8.4)', '   ', m(3,:)
v = [ex + ey + ez, ex - ey]
do i = 1, 2
   call v(i)%mirror(matrix=m)
   print '(A,I0,A,3F8.4)', 'v', i, ' mirrored ', v(i)%x, v(i)%y, v(i)%z
enddo
endprogram mirror_matrix_reuse

$ mirror_matrix_reuse
M    0.3333 -0.6667 -0.6667
    -0.6667  0.3333 -0.6667
    -0.6667 -0.6667  0.3333
v1 mirrored  -1.0000 -1.0000 -1.0000
v2 mirrored   1.0000 -1.0000  0.0000

The Householder matrix $I-2\hat{n}{\hat{n}}^T$: symmetric, its own inverse. A vector along the normal is reversed, a vector in the plane stays.

Any matrix applied to a vector

program matrix_product
!< Cookbook: any 3 x 3 matrix applied to a vector.
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: v, w
real(R8P)    :: shear(3,3)

shear = reshape([1._R8P, 0._R8P, 0._R8P, &      ! by columns: x gains half of y
                 0.5_R8P, 1._R8P, 0._R8P, &
                 0._R8P, 0._R8P, 1._R8P], [3, 3])
v = ex + 2 * ey + 3 * ez
w = shear .matrix. v
call w%printf(prefix='shear .matrix. v = ')
endprogram matrix_product

$ matrix_product
shear .matrix. v = +2.0, +2.0, +3.0

.matrix. takes any 3 x 3 real matrix of the kind of the vector, on the left.

Comparisons

Comparing two vectors

program compare
!< Cookbook: comparing two vectors; by length, and == also by direction.
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: a, b, c

a = 3 * ex
b = 3 * ey                               ! as long as a, another direction
c = 2 * ex + 2 * ey + 2 * ez             ! longer
print '(A,L1)', 'a == b  ', a == b
print '(A,L1)', 'a /= b  ', a /= b
print '(A,L1)', 'a <= b  ', a <= b
print '(A,L1)', 'a >= b  ', a >= b
print '(A,L1)', 'a <  c  ', a < c
print '(A,L1)', 'c >  b  ', c > b
print '(A,L1)', 'a == 3 ex ', a == 3 * ex
endprogram compare

$ compare
a == b  F
a /= b  T
a <= b  T
a >= b  T
a <  c  T
c >  b  T
a == 3 ex T

<, <=, >, >= compare lengths: a and b, as long as each other, are both <= and >=. == and /= compare lengths and directions: a == b is false. Both are exact, without a tolerance.

Comparing the length of a vector with a number

program compare_number
!< Cookbook: comparing the length of a vector with a number.
use penf, only : I4P, R4P, R8P
use vecfor, only : ex, ey, vector
implicit none
type(vector) :: v

v = 3 * ex + 4 * ey                      ! its length is 5
print '(A,L1)', 'v == 5        ', v == 5
print '(A,L1)', 'v <  5.5_R8P  ', v < 5.5_R8P
print '(A,L1)', '4_I4P < v     ', 4_I4P < v
print '(A,L1)', 'v >= 6._R4P   ', v >= 6._R4P
endprogram compare_number

$ compare_number
v == 5        T
v <  5.5_R8P  T
4_I4P < v     T
v >= 6._R4P   F

The number may be of any integer or real kind, on either side.

Arrays and kinds

Whole arrays of vectors at once

program arrays
!< Cookbook: whole arrays of vectors at once.
use penf, only : R8P
use vecfor, only : ex, ey, ez, normL2, vector
implicit none
type(vector) :: p(5), c
integer      :: i

do i = 1, 5
   p(i) = i * ex + (6 - i) * ey               ! five points on a line
enddo
print '(A,5F6.2)', 'x of the points   ', p%x                       ! a component of every vector
print '(A,5F6.2)', 'their lengths     ', normL2(p)                 ! elemental functions: one result each
print '(A,5F6.2)', 'distances from ez ', p%distance_to_line(pt1=0 * ex, pt2=ez)
p = p + ez                                                          ! elemental operators
print '(A,5F6.2)', 'z after + ez      ', p%z
print '(A,I0)',    'longer than 5     ', count(p > 5)
c = sum(p%x) * ex + sum(p%y) * ey + sum(p%z) * ez                   ! no sum(p): sum the components
call c%printf(prefix='their sum         ')
endprogram arrays

$ arrays
x of the points     1.00  2.00  3.00  4.00  5.00
their lengths       5.10  4.47  4.24  4.47  5.10
distances from ez   5.10  4.47  4.24  4.47  5.10
z after + ez        1.00  1.00  1.00  1.00  1.00
longer than 5     2
their sum         +15.0, +15.0, +5.0

Operators and functions are elemental; p%x is the array of the x components. Fortran has no sum of an array of a derived type: sum the components.

Vectors of each precision

program kinds
!< Cookbook: vectors of each precision, their versors and their free functions.
use penf, only : R4P, R8P
use vecfor, only : angle_R4P, ex_R4P, ey_R4P, ex_R8P, normL2_R8P, vector, vector_R4P, vector_R8P, vector_R16P
implicit none
type(vector_R4P)  :: a      ! real(R4P) components
type(vector_R8P)  :: b      ! real(R8P) components
type(vector_R16P) :: c      ! real(R16P) components
type(vector)      :: d      ! the default: real(R8P) components, but not the type vector_R8P

a = ex_R4P + ey_R4P
b = 3 * ex_R8P
print '(A,ES14.7)', 'angle_R4P(a, ex_R4P) = ', angle_R4P(a, ex_R4P)
print '(A,ES14.7)', 'normL2_R8P(b)        = ', normL2_R8P(b)
print '(A,3I4)', 'kinds of a, b, c:      ', kind(a%x), kind(b%x), kind(c%x)
d%x = b%x ; d%y = b%y ; d%z = b%z             ! d = b does not compile: two different types
call d%printf(prefix='d = ')
endprogram kinds

$ kinds
angle_R4P(a, ex_R4P) =  7.8539819E-01
normL2_R8P(b)        =  3.0000000E+00
kinds of a, b, c:         4   8   8
d = +3.0, +0.0, +0.0

Each kind has its type, its versors and its free functions, all exported by vecfor. vector and vector_R8P have the same kind but are two types: copy the components. vector_R16P has the kind 8 here, double precision: see the next recipe.

True quadruple precision

program quad_precision
!< Cookbook: true quadruple precision, with VecFor built with -DPENF_R16P.
use penf, only : R16P
use vecfor, only : ex_R16P, ey_R16P, vector_R16P
implicit none
type(vector_R16P) :: v

v = ex_R16P + ey_R16P
print '(A,I0)', 'digits of vector_R16P: ', precision(v%x)
print '(A,F38.34)', 'its length:  ', v%normL2()
print '(A,F38.34)', 'sqrt(2):     ', sqrt(2._R16P)
v = v / 3
call v%printf(prefix='(ex + ey) / 3 = ')
endprogram quad_precision

$ quad_precision
digits of vector_R16P: 33
its length:    1.4142135623730950488016887242096980
sqrt(2):       1.4142135623730950488016887242096980
(ex + ey) / 3 = +0.3333333333333333333333333333333333, +0.3333333333333333333333333333333333, +0.0

Built against VecFor and PENF compiled with -DPENF_R16P (see Precision and kinds): 33 significant digits.

A method or a free function

program two_spellings
!< Cookbook: a method of the vector, or a free function: the same result.
use penf, only : R8P
use vecfor, only : angle, distance_to_line, ex, ey, ez, face_normal3, is_collinear, normL2, normalized, sq_norm, vector
implicit none
type(vector) :: v, u, w

v = 3 * ex + 4 * ey
print '(A,2F8.4)', 'normL2:           ', v%normL2(), normL2(v)
print '(A,2F8.4)', 'sq_norm:          ', v%sq_norm(), sq_norm(v)
print '(A,2F8.4)', 'angle:            ', v%angle(ex), angle(v, ex)
print '(A,2F8.4)', 'distance_to_line: ', ez%distance_to_line(pt1=ex, pt2=ey), distance_to_line(ez, pt1=ex, pt2=ey)
print '(A,2L2)',   'is_collinear:     ', v%is_collinear(pt1=ex, pt2=ey), is_collinear(v, pt1=ex, pt2=ey)
u = v%normalized() ; w = normalized(v)
print '(A,L2)',    'normalized:       ', u == w
u = v%face_normal3(pt1=ex, pt2=ey, pt3=ez) ; w = face_normal3(pt1=ex, pt2=ey, pt3=ez)
print '(A,L2)',    'face_normal3:     ', u == w
endprogram two_spellings

$ two_spellings
normL2:             5.0000  5.0000
sq_norm:           25.0000 25.0000
angle:              0.9273  0.9273
distance_to_line:   1.2247  1.2247
is_collinear:      F F
normalized:        T
face_normal3:      T

The methods that return a value are free functions too, with the vector as first argument: use them on expressions (normL2(a - b)), where (a - b)%normL2() is not Fortran.

Input and output

Printing a vector

program printing
!< Cookbook: printing a vector, with a prefix, a separator and a suffix, on any unit.
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: v
integer      :: u, ios

v = 0.5_R8P * ex - 2 * ey + ez / 3
call v%printf
call v%printf(prefix='v = (', sep=', ', suffix=')')
call v%printf(prefix='', sep=' ')
print '(A,3F8.3)', 'a format of your own:', v%x, v%y, v%z
open(newunit=u, file='vector.txt', status='replace')
call v%printf(unit=u, prefix='in a file: ', iostat=ios)
close(u)
endprogram printing

$ printing && cat vector.txt
+0.5, -2.0, +0.3333333333333333
v = (+0.5, -2.0, +0.3333333333333333)
+0.5 -2.0 +0.3333333333333333
a format of your own:   0.500  -2.000   0.333
in a file: +0.5, -2.0, +0.3333333333333333

prefix, sep and suffix frame the components; unit writes to any open unit. For a fixed format, print the components.

A binary file of vectors

program binary_file
!< Cookbook: vectors saved to and loaded from an unformatted file: the exact bits.
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: v, w
integer      :: u

v = ex / 3 + ey / 7 - ez / 11
open(newunit=u, file='v.bin', form='unformatted', status='replace')
call v%save_into_file(unit=u)
rewind(u)
call w%load_from_file(unit=u)
close(u)
print '(A,L1)', 'loaded == saved: ', w == v
call w%printf(prefix='w = ')
endprogram binary_file

$ binary_file
loaded == saved: T
w = +0.3333333333333333, +0.14285714285714285, -0.09090909090909091

Unformatted: the exact bits.

A text file of vectors

program text_file
!< Cookbook: vectors saved to and loaded from a text file, with a format.
use penf, only : R8P
use vecfor, only : ex, ey, ez, normL2, vector
implicit none
type(vector) :: p(3), q
integer      :: u, i

p = [ex, 2 * ey, ex + ey + ez / 3]
open(newunit=u, file='points.txt', status='replace')
do i = 1, 3
   call p(i)%save_into_file(unit=u, fmt='(3F10.5)')
enddo
rewind(u)
do i = 1, 3
   call q%load_from_file(unit=u, fmt='(3F10.5)')
   print '(A,I0,A,L1,A,ES8.1)', 'point ', i, ' read back equal: ', q == p(i), ', difference ', normL2(q - p(i))
enddo
close(u)
endprogram text_file

$ text_file && cat points.txt
point 1 read back equal: T, difference  0.0E+00
point 2 read back equal: T, difference  0.0E+00
point 3 read back equal: F, difference  3.3E-06
   1.00000   0.00000   0.00000
   0.00000   2.00000   0.00000
   1.00000   1.00000   0.33333

fmt is a format for three reals; the file keeps the digits of the format, no more: 1/3 comes back $3.3\times {10}^{-6}$ off.

The n-th vector of a stream file

program stream_file
!< Cookbook: the n-th vector of a stream file, read straight away.
use, intrinsic :: iso_fortran_env, only : int64
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: v
integer      :: u, i
integer      :: bytes

open(newunit=u, file='points.raw', access='stream', form='unformatted', status='replace')
do i = 1, 10
   v = i * ex + i**2 * ey
   call v%save_into_file(unit=u)
enddo
bytes = v%iolen()                                   ! the file storage units of one vector
print '(A,I0)', 'bytes per vector: ', bytes
call v%load_from_file(unit=u, pos=6 * int(bytes, int64) + 1)
call v%printf(prefix='vector 7: ')
close(u)
endprogram stream_file

$ stream_file
bytes per vector: 24
vector 7: +7.0, +49.0, +0.0

pos is in file storage units, iolen the storage units of one vector; vector n starts at (n - 1) * iolen + 1.

A direct access file of vectors

program direct_access
!< Cookbook: a direct access file of vectors, its record length from iolen.
use penf, only : R8P
use vecfor, only : ex, ey, iolen, vector
implicit none
type(vector) :: v
integer      :: u, i

open(newunit=u, file='points.dat', access='direct', form='unformatted', recl=iolen(v), status='replace')
do i = 1, 10
   v = i * ex - i * ey
   write(u, rec=i) v%x, v%y, v%z                    ! save_into_file has no rec=: write the components
enddo
read(u, rec=4) v%x, v%y, v%z
call v%printf(prefix='record 4: ')
close(u)
endprogram direct_access

$ direct_access
record 4: +4.0, -4.0, +0.0

iolen is the record length. save_into_file and load_from_file have no rec=: read and write the components.

An I/O error caught

program io_error
!< Cookbook: an I/O error caught, instead of stopping the program.
use vecfor, only : vector
implicit none
type(vector)       :: v
integer            :: u, ios
character(len=200) :: msg

open(newunit=u, file='bad.txt', status='replace')
write(u, '(A)') '    1.0000    2.0000      oops'
rewind(u)
call v%load_from_file(unit=u, fmt='(3F10.4)', iostat=ios, iomsg=msg)
print '(A,I0)', 'iostat: ', ios
print '(2A)',   'iomsg:  ', trim(msg)
call v%load_from_file(unit=u, fmt='(3F10.4)', iostat=ios, iomsg=msg)
print '(A,I0,2A)', 'at the end: iostat ', ios, ', ', trim(msg)
close(u)
endprogram io_error

$ io_error
iostat: 5010
iomsg:  Bad value during floating point read
at the end: iostat -1, End of file

Pass iostat (and iomsg) to every save_into_file and load_from_file: without them an error is ignored, silently.

On the GPU

The device-callable routines in an OpenACC loop

program device_routines
!< Cookbook: the seven device-callable routines, inside an OpenACC loop; the same results as the operators.
use penf, only : R8P
use vecfor, only : assign_vector_oac, crossproduct_oac, dotproduct_oac, ex, ey, ez, R8P_mul_vector_oac, vector, &
                   vector_mul_R8P_oac, vector_sub_vector_oac, vector_sum_vector_oac
implicit none
integer, parameter :: n = 1000
type(vector)       :: a(n), b(n), c(n), d(n), e(n), f(n), g(n), h(n)
real(R8P)          :: dots(n)
integer            :: i

do i = 1, n
   a(i) = i * ex + ey / i
   b(i) = ez - i * ey
enddo
!$acc parallel loop copyin(a, b) copyout(c, d, e, f, g, h, dots)
do i = 1, n
   call vector_sum_vector_oac(a(i), b(i), c(i))       ! c = a + b
   call vector_sub_vector_oac(a(i), b(i), d(i))       ! d = a - b
   call crossproduct_oac(a(i), b(i), e(i))            ! e = a .cross. b
   call R8P_mul_vector_oac(2._R8P, a(i), f(i))        ! f = 2 * a
   call vector_mul_R8P_oac(a(i), 0.5_R8P, g(i))       ! g = a * 0.5
   call assign_vector_oac(h(i), b(i))                 ! h = b
   dots(i) = dotproduct_oac(a(i), b(i))               ! a .dot. b
enddo
print '(A,L1)', 'a + b       ', all(c == a + b)
print '(A,L1)', 'a - b       ', all(d == a - b)
print '(A,L1)', 'a .cross. b ', all(e == (a .cross. b))
print '(A,L1)', '2 * a       ', all(f == 2._R8P * a)
print '(A,L1)', 'a * 0.5     ', all(g == a * 0.5_R8P)
print '(A,L1)', 'h = b       ', all(h == b)
print '(A,L1)', 'a .dot. b   ', all(dots == (a .dot. b))
endprogram device_routines

$ device_routines
a + b       T
a - b       T
a .cross. b T
2 * a       T
a * 0.5     T
h = b       T
a .dot. b   T

The seven _oac routines and the operators they stand for; this output is the gfortran build, without -fopenacc. See tutorial, chapter 11 and Device-callable API.