1+ ! >
2+ ! ! @file m_lapack_example.f90
3+ ! ! @brief Contains module m_lapack_example
4+
5+ ! > @brief This module demonstrates the use of LAPACK in MFC post_process.
6+ ! ! It provides example routines that show how to use LAPACK for
7+ ! ! common linear algebra operations like solving linear systems.
8+ module m_lapack_example
9+
10+ use m_global_parameters ! < Global parameters for the code
11+ use m_mpi_proxy ! < Message passing interface (MPI) module proxy
12+
13+ implicit none
14+
15+ private ; public :: s_lapack_example_solve_linear_system, &
16+ s_lapack_example_eigenvalues
17+
18+ contains
19+
20+ ! > @brief Example subroutine demonstrating LAPACK usage for solving
21+ ! ! a linear system Ax = b using DGESV/SGESV
22+ ! ! This routine shows how to use LAPACK with MFC's precision system
23+ impure subroutine s_lapack_example_solve_linear_system ()
24+
25+ ! Local variables for the linear system
26+ integer , parameter :: n = 3 ! Size of the system
27+ real (wp), dimension (n, n) :: A ! Coefficient matrix
28+ real (wp), dimension (n) :: b ! Right-hand side vector
29+ real (wp), dimension (n) :: x ! Solution vector
30+
31+ ! LAPACK variables
32+ integer , dimension (n) :: ipiv ! Pivot indices
33+ integer :: info ! Return status
34+ integer :: nrhs = 1 ! Number of right-hand sides
35+
36+ ! Only run on the root process to avoid duplicate output
37+ if (proc_rank /= 0 ) return
38+
39+ ! Set up a simple 3x3 linear system: Ax = b
40+ ! Example:
41+ ! 2x + y + z = 8
42+ ! x + 3y + z = 10
43+ ! x + y + 4z = 16
44+ A(1 , :) = [2.0_wp , 1.0_wp , 1.0_wp ]
45+ A(2 , :) = [1.0_wp , 3.0_wp , 1.0_wp ]
46+ A(3 , :) = [1.0_wp , 1.0_wp , 4.0_wp ]
47+
48+ b = [8.0_wp , 10.0_wp , 16.0_wp ]
49+
50+ print * , " === LAPACK Linear System Solver Example ==="
51+ print * , " Solving the system Ax = b where:"
52+ print * , " A = [2 1 1; 1 3 1; 1 1 4]"
53+ print * , " b = [8; 10; 16]"
54+ print * , " "
55+
56+ ! Copy b to x (LAPACK will overwrite the right-hand side with solution)
57+ x = b
58+
59+ ! Call appropriate LAPACK routine based on precision
60+ #ifdef MFC_SINGLE_PRECISION
61+ call sgesv(n, nrhs, A, n, ipiv, x, n, info)
62+ print * , " Using single precision LAPACK (SGESV)"
63+ #else
64+ call dgesv(n, nrhs, A, n, ipiv, x, n, info)
65+ print * , " Using double precision LAPACK (DGESV)"
66+ #endif
67+
68+ ! Check for success
69+ if (info == 0 ) then
70+ print * , " Linear system solved successfully!"
71+ print * , " Solution: x = [" , x(1 ), " , " , x(2 ), " , " , x(3 ), " ]"
72+ print * , " Expected: x = [1, 2, 3]"
73+ else if (info < 0 ) then
74+ print * , " LAPACK error: argument " , - info, " had an illegal value"
75+ else
76+ print * , " LAPACK error: matrix is singular, solution could not be computed"
77+ end if
78+
79+ print * , " === End LAPACK Example ==="
80+ print * , " "
81+
82+ end subroutine s_lapack_example_solve_linear_system
83+
84+ ! > @brief Example subroutine demonstrating LAPACK usage for computing
85+ ! ! eigenvalues of a symmetric matrix using DSYEV/SSYEV
86+ impure subroutine s_lapack_example_eigenvalues ()
87+
88+ ! Local variables for eigenvalue computation
89+ integer , parameter :: n = 3 ! Size of the matrix
90+ real (wp), dimension (n, n) :: A ! Symmetric matrix
91+ real (wp), dimension (n) :: w ! Eigenvalues
92+ real (wp), dimension (3 * n) :: work ! Work array
93+ integer :: lwork = 3 * n ! Size of work array
94+ integer :: info ! Return status
95+ character :: jobz = ' N' ! Compute eigenvalues only
96+ character :: uplo = ' U' ! Upper triangular part of A
97+
98+ ! Only run on the root process to avoid duplicate output
99+ if (proc_rank /= 0 ) return
100+
101+ ! Set up a simple symmetric 3x3 matrix
102+ A(1 , :) = [4.0_wp , 1.0_wp , 1.0_wp ]
103+ A(2 , :) = [1.0_wp , 4.0_wp , 1.0_wp ]
104+ A(3 , :) = [1.0_wp , 1.0_wp , 4.0_wp ]
105+
106+ print * , " === LAPACK Eigenvalue Example ==="
107+ print * , " Computing eigenvalues of symmetric matrix:"
108+ print * , " A = [4 1 1; 1 4 1; 1 1 4]"
109+ print * , " "
110+
111+ ! Call appropriate LAPACK routine based on precision
112+ #ifdef MFC_SINGLE_PRECISION
113+ call ssyev(jobz, uplo, n, A, n, w, work, lwork, info)
114+ print * , " Using single precision LAPACK (SSYEV)"
115+ #else
116+ call dsyev(jobz, uplo, n, A, n, w, work, lwork, info)
117+ print * , " Using double precision LAPACK (DSYEV)"
118+ #endif
119+
120+ ! Check for success
121+ if (info == 0 ) then
122+ print * , " Eigenvalues computed successfully!"
123+ print * , " Eigenvalues: [" , w(1 ), " , " , w(2 ), " , " , w(3 ), " ]"
124+ print * , " Expected: [2, 5, 5] (approximately)"
125+ else if (info < 0 ) then
126+ print * , " LAPACK error: argument " , - info, " had an illegal value"
127+ else
128+ print * , " LAPACK error: algorithm failed to converge"
129+ end if
130+
131+ print * , " === End LAPACK Eigenvalue Example ==="
132+ print * , " "
133+
134+ end subroutine s_lapack_example_eigenvalues
135+
136+ end module m_lapack_example
0 commit comments