Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module mtx_lib
21 :
22 : use const_def, only: dp
23 :
24 : implicit none
25 :
26 : contains
27 :
28 : ! mesa includes sources for a subset of BLAS and dble.
29 : ! you can use those, or, better yet, you can use a package optimized
30 : ! for your machine such as GotoBLAS or Intel's MKL.
31 : ! see utils/makefile_header for details.
32 :
33 : ! see mtx/blas_src for the subset of BLAS routines included in mtx_lib
34 : ! see mtx/dble_src for the subset of dble routines included in mtx_lib
35 :
36 : ! subroutines for dense and banded matrix decompositions and solves
37 :
38 : include "mtx_dble_decsol.dek" ! dble versions
39 :
40 : !> Wraps the lapack DGBSVX routine for banded matrices to be used with block tridiagonal matrices.
41 : !!
42 : !! @params ublk The upper blocks. Shape is (nvars, nvars, nblocks-1).
43 : !! @params lblk The lower blocks. Shape is (nvars, nvars, nblocks-1).
44 : !! @params dblk The diagonal blocks. Shape is (nvars, nvars, nblocks).
45 : !! @params nblocks The number of blocks on the diagonal.
46 : !! @params nvars The width and height of each block (e.g. each block is nvars by nvars).
47 : !! @params x Array of shape (nblocks,nvar). Stores the result of solving Matrix.x = b.
48 : !! @params b Array of shape (nblocks,nvar). The right-hand side from the definition of x.
49 : !! @params ierr Integer-valued error code.
50 0 : subroutine DGBSVX_block_tridiagonal_padded(ublk, lblk, dblk, x, b, nblocks, nvar, rcond, ierr)
51 : use DGBSVX_wrapper, only: DGBSVX_tridiagonal_wrapper
52 :
53 : ! Inputs
54 : real(dp), dimension(:, :, :), intent(in) :: ublk, lblk, dblk
55 : real(dp), dimension(:, :), intent(in) :: b
56 : integer, intent(in) :: nblocks, nvar
57 :
58 : ! Intermediates
59 0 : real(dp) :: pre_conditioner(nvar, nblocks)
60 :
61 : ! Outputs
62 : real(dp), dimension(:, :), intent(out) :: x
63 : real(dp), intent(out) :: rcond
64 : integer, intent(out) :: ierr
65 :
66 0 : call DGBSVX_tridiagonal_wrapper(ublk(:,:,1:nblocks-1), lblk(:,:,2:nblocks), dblk, pre_conditioner, x, b, nblocks, nvar, rcond, ierr)
67 :
68 0 : end subroutine DGBSVX_block_tridiagonal_padded
69 :
70 : !> Wraps the lapack DGBSVX routine for banded matrices to be used with block tridiagonal matrices.
71 : !! This version returns the pre-conditioning array which was used to rescale each row.
72 : !!
73 : !! @params ublk The upper blocks. Shape is (nvars, nvars, nblocks-1).
74 : !! @params lblk The lower blocks. Shape is (nvars, nvars, nblocks-1).
75 : !! @params dblk The diagonal blocks. Shape is (nvars, nvars, nblocks).
76 : !! @params nblocks The number of blocks on the diagonal.
77 : !! @params nvars The width and height of each block (e.g. each block is nvars by nvars).
78 : !! @params x Array of shape (nblocks,nvar). Stores the result of solving Matrix.x = b.
79 : !! @params b Array of shape (nblocks,nvar). The right-hand side from the definition of x.
80 : !! @params pre_conditioner Pre-conditioning vector. Computed in this routine.
81 : !! @params ierr Integer-valued error code.
82 0 : subroutine DGBSVX_block_tridiagonal_pc(ublk, lblk, dblk, pre_conditioner, x, b, nblocks, nvar, ierr)
83 : use DGBSVX_wrapper, only: DGBSVX_tridiagonal_wrapper
84 :
85 : ! Inputs
86 : real(dp), dimension(:, :, :), intent(in) :: ublk, lblk, dblk
87 : real(dp), dimension(:, :), intent(in) :: b
88 : integer, intent(in) :: nblocks, nvar
89 :
90 : ! Intermediates
91 : real(dp) :: rcond
92 :
93 : ! Outputs
94 : real(dp), dimension(:, :), intent(out) :: pre_conditioner
95 : real(dp), dimension(:, :), intent(out) :: x
96 : integer, intent(out) :: ierr
97 :
98 0 : call DGBSVX_tridiagonal_wrapper(ublk, lblk, dblk, pre_conditioner, x, b, nblocks, nvar, rcond, ierr)
99 :
100 0 : end subroutine DGBSVX_block_tridiagonal_pc
101 :
102 : !> Wraps the lapack DGBSVX routine for banded matrices to be used with block tridiagonal matrices.
103 : !!
104 : !! @params ublk The upper blocks. Shape is (nvars, nvars, nblocks-1).
105 : !! @params lblk The lower blocks. Shape is (nvars, nvars, nblocks-1).
106 : !! @params dblk The diagonal blocks. Shape is (nvars, nvars, nblocks).
107 : !! @params nblocks The number of blocks on the diagonal.
108 : !! @params nvars The width and height of each block (e.g. each block is nvars by nvars).
109 : !! @params x Array of shape (nblocks,nvar). Stores the result of solving Matrix.x = b.
110 : !! @params b Array of shape (nblocks,nvar). The right-hand side from the definition of x.
111 : !! @params ierr Integer-valued error code.
112 0 : subroutine DGBSVX_block_tridiagonal(ublk, lblk, dblk, x, b, nblocks, nvar, ierr)
113 : use DGBSVX_wrapper, only: DGBSVX_tridiagonal_wrapper
114 :
115 : ! Inputs
116 : real(dp), dimension(:, :, :), intent(in) :: ublk, lblk, dblk
117 : real(dp), dimension(:, :), intent(in) :: b
118 : integer, intent(in) :: nblocks, nvar
119 :
120 : ! Intermediates
121 0 : real(dp) :: pre_conditioner(nvar, nblocks)
122 0 : real(dp) :: rcond
123 :
124 : ! Outputs
125 : real(dp), dimension(:, :), intent(out) :: x
126 : integer, intent(out) :: ierr
127 :
128 0 : call DGBSVX_tridiagonal_wrapper(ublk, lblk, dblk, pre_conditioner, x, b, nblocks, nvar, rcond, ierr)
129 :
130 0 : end subroutine DGBSVX_block_tridiagonal
131 :
132 : !> Wraps the lapack DGBSVX routine for banded matrices with a preconditioner
133 : !! This version returns the pre-conditioning array which was used to rescale each row.
134 : !!
135 : !! @param matrix_size The number of rows or columns in the square matrix.
136 : !! @params n_upper_bands The number of bands above the diagonal.
137 : !! @params n_lower_bands The number of bands below the diagonal.
138 : !! @params bands Array of shape (n_upper_bands + n_lower_bands + 1, matrix_size). Stores the bands going from upper-most to lower-most, aligned so the first entry in each band has index 1.
139 : !! @params x Array of shape (matrix_size). Stores the result of solving Matrix.x = b.
140 : !! @params b Array of shape (matrix_size). The right-hand side from the definition of x.
141 : !! @params pre_conditioner Pre-conditioning vector. Computed in this routine.
142 : !! @params ierr Integer-valued error code.
143 0 : subroutine DGBSVX_banded_pc(matrix_size, n_upper_bands, n_lower_bands, bands, x, b, pre_conditioner, ierr)
144 : use DGBSVX_wrapper, only: DGBSVX_banded_wrapper
145 :
146 : ! Inputs
147 : real(dp), dimension(:, :), intent(in) :: bands
148 : real(dp), dimension(:), intent(in) :: b
149 : integer, intent(in) :: matrix_size, n_upper_bands, n_lower_bands
150 :
151 : ! Outputs
152 : real(dp), dimension(:), intent(out) :: x
153 : real(dp), dimension(:), intent(out) :: pre_conditioner
154 : integer, intent(out) :: ierr
155 :
156 0 : call DGBSVX_banded_wrapper(matrix_size, n_upper_bands, n_lower_bands, bands, x, b, pre_conditioner, ierr)
157 :
158 0 : end subroutine DGBSVX_banded_pc
159 :
160 : !> Wraps the lapack DGBSVX routine for banded matrices with a preconditioner
161 : !!
162 : !! @param matrix_size The number of rows or columns in the square matrix.
163 : !! @params n_upper_bands The number of bands above the diagonal.
164 : !! @params n_lower_bands The number of bands below the diagonal.
165 : !! @params bands Array of shape (n_upper_bands + n_lower_bands + 1, matrix_size). Stores the bands going from upper-most to lower-most, aligned so the first entry in each band has index 1.
166 : !! @params x Array of shape (matrix_size). Stores the result of solving Matrix.x = b.
167 : !! @params b Array of shape (matrix_size). The right-hand side from the definition of x.
168 : !! @params ierr Integer-valued error code.
169 0 : subroutine DGBSVX_banded(matrix_size, n_upper_bands, n_lower_bands, bands, x, b, ierr)
170 : use DGBSVX_wrapper, only: DGBSVX_banded_wrapper
171 :
172 : ! Inputs
173 : real(dp), dimension(:, :), intent(in) :: bands
174 : real(dp), dimension(:), intent(in) :: b
175 : integer, intent(in) :: matrix_size, n_upper_bands, n_lower_bands
176 :
177 : ! Intermediates
178 0 : real(dp) :: pre_conditioner(matrix_size)
179 :
180 : ! Outputs
181 : real(dp), dimension(:), intent(out) :: x
182 : integer, intent(out) :: ierr
183 :
184 0 : call DGBSVX_banded_wrapper(matrix_size, n_upper_bands, n_lower_bands, bands, x, b, pre_conditioner, ierr)
185 :
186 0 : end subroutine DGBSVX_banded
187 :
188 : ! sometimes you just need a null version of a routine
189 : include "mtx_null_decsol.dek"
190 :
191 : ! sometimes you need to debug a jacobian by saving it to plotting data files
192 : include "mtx_debug_decsol.dek"
193 :
194 : ! sparse matrices come in many formats.
195 : ! for example, compressed row sparse format is used by SPARSKIT,
196 : ! while compressed column sparse format is used by Super_LU.
197 : ! here are conversion routines for these two options.
198 : include "mtx_formats.dek"
199 :
200 0 : subroutine mtx_write_hbcode1(iounit, n, nnzero, values, rowind, colptr, ierr)
201 : use mtx_support, only: write_hbcode1
202 : integer, intent(in) :: iounit, n, nnzero
203 : integer :: rowind(:) ! (nnzero)
204 : integer :: colptr(:) ! (n+1)
205 : real(dp) :: values(:) ! (nnzero)
206 : integer, intent(out) :: ierr
207 0 : call write_hbcode1(iounit, n, n, nnzero, values, rowind, colptr, ierr)
208 0 : end subroutine mtx_write_hbcode1
209 :
210 0 : subroutine mtx_write_block_tridiagonal(iounit, nvar, nblk, lblk, dblk, ublk, ierr)
211 : use mtx_support, only: write_block_tridiagonal
212 : integer, intent(in) :: iounit, nvar, nblk
213 : real(dp), intent(in), dimension(:, :, :) :: lblk, dblk, ublk
214 : integer, intent(out) :: ierr
215 0 : call write_block_tridiagonal(iounit, nvar, nblk, lblk, dblk, ublk, ierr)
216 0 : end subroutine mtx_write_block_tridiagonal
217 :
218 1 : subroutine mtx_read_block_tridiagonal(iounit, nvar, nblk, lblk1, dblk1, ublk1, ierr)
219 : use mtx_support, only: read_block_tridiagonal
220 : integer, intent(in) :: iounit
221 : integer, intent(out) :: nvar, nblk
222 : real(dp), pointer, dimension(:) :: lblk1, dblk1, ublk1 ! =(nvar,nvar,nblk) will be allocated
223 : integer, intent(out) :: ierr
224 1 : call read_block_tridiagonal(iounit, nvar, nblk, lblk1, dblk1, ublk1, ierr)
225 1 : end subroutine mtx_read_block_tridiagonal
226 :
227 : ! BCYCLIC multi-thread block tridiagonal
228 : include "mtx_bcyclic_dble_decsol.dek"
229 : ! S.P.Hirshman, K.S.Perumalla, V.E.Lynch, & R.Sanchez,
230 : ! BCYCLIC: A parallel block tridiagonal matrix cyclic solver,
231 : ! J. Computational Physics, 229 (2010) 6392-6404.
232 :
233 1 : subroutine block_dble_mv(nvar, nz, lblk, dblk, ublk, b, prod)
234 : ! set prod = A*b with A = block tridiagonal given by lblk, dblk, ublk
235 : use mtx_support, only: do_block_dble_mv
236 : integer, intent(in) :: nvar, nz
237 : real(dp), pointer, dimension(:, :, :), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz)
238 : real(dp), pointer, dimension(:, :), intent(in) :: b ! (nvar,nz)
239 : real(dp), pointer, dimension(:, :), intent(inout) :: prod ! (nvar,nz)
240 1 : call do_block_dble_mv(nvar, nz, lblk, dblk, ublk, b, prod)
241 1 : end subroutine block_dble_mv
242 :
243 0 : subroutine multiply_xa(n, A1, x, b)
244 : ! calculates b = x*A
245 : use mtx_support, only: do_multiply_xa
246 : integer, intent(in) :: n
247 : real(dp), pointer, intent(in) :: A1(:) ! =(n, n)
248 : real(dp), pointer, intent(in) :: x(:) ! (n)
249 : real(dp), pointer, intent(inout) :: b(:) ! (n)
250 0 : call do_multiply_xa(n, A1, x, b)
251 0 : end subroutine multiply_xa
252 :
253 0 : subroutine block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1)
254 : ! calculates b = x*A
255 : use mtx_support, only: do_block_multiply_xa
256 : integer, intent(in) :: nvar, nz
257 : real(dp), dimension(:), intent(in), pointer :: lblk1, dblk1, ublk1 ! =(nvar,nvar,nz)
258 : real(dp), intent(in), pointer :: x1(:) ! =(nvar,nz)
259 : real(dp), intent(inout), pointer :: b1(:) ! =(nvar,nz)
260 0 : call do_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1)
261 0 : end subroutine block_multiply_xa
262 :
263 20 : subroutine band_multiply_xa(n, kl, ku, ab1, ldab, x, b)
264 : ! calculates b = x*a = transpose(a)*x
265 : use mtx_support, only: do_band_multiply_xa
266 : integer, intent(in) :: n
267 : ! the number of linear equations, i.e., the order of the
268 : ! matrix a. n >= 0.
269 : integer, intent(in) :: kl
270 : ! the number of subdiagonals within the band of a. kl >= 0.
271 : integer, intent(in) :: ku
272 : ! the number of superdiagonals within the band of a. ku >= 0.
273 : integer, intent(in) :: ldab
274 : ! the leading dimension of the array ab. ldab >= kl+ku+1.
275 : real(dp), intent(in), pointer :: ab1(:) ! =(ldab, n)
276 : ! the matrix a in band storage, in rows 1 to kl+ku+1;
277 : ! the j-th column of a is stored in the j-th column of the
278 : ! array ab as follows:
279 : ! ab(ku+1+i-j, j) = a(i, j) for max(1, j-ku)<=i<=min(n, j+kl)
280 : real(dp), intent(in), pointer :: x(:) ! (n)
281 : ! the input vector to be multiplied by the matrix.
282 : real(dp), intent(inout), pointer :: b(:) ! (n)
283 : ! on exit, set to matrix product of x*a = b
284 20 : call do_band_multiply_xa(n, kl, ku, ab1, ldab, x, b)
285 20 : end subroutine band_multiply_xa
286 :
287 : include "mtx_lapack95.dek"
288 :
289 : ! utilities for working with jacobians
290 : include "mtx_jac.dek"
291 :
292 : ! the following call dble routines to estimate matrix condition numbers.
293 : include "mtx_rcond.dek"
294 :
295 4 : integer function decsol_option(which_decsol_option, ierr)
296 : use mtx_def
297 : character(len=*), intent(in) :: which_decsol_option
298 : integer, intent(out) :: ierr
299 : character(len=64) :: option
300 4 : ierr = 0
301 4 : option = which_decsol_option
302 :
303 4 : if (option == 'lapack') then
304 : decsol_option = lapack
305 :
306 0 : else if (option == 'bcyclic_dble') then
307 : decsol_option = bcyclic_dble
308 :
309 : else
310 0 : ierr = -1
311 0 : decsol_option = -1
312 : end if
313 4 : end function decsol_option
314 :
315 4 : subroutine decsol_option_str(which_decsol_option, decsol_option, ierr)
316 : use mtx_def
317 : integer, intent(in) :: which_decsol_option
318 : character(len=*), intent(out) :: decsol_option
319 : integer, intent(out) :: ierr
320 4 : ierr = 0
321 :
322 4 : if (which_decsol_option == lapack) then
323 1 : decsol_option = 'lapack'
324 3 : else if (which_decsol_option == bcyclic_dble) then
325 1 : decsol_option = 'bcyclic_dble'
326 :
327 : else
328 2 : ierr = -1
329 2 : decsol_option = ''
330 : end if
331 :
332 4 : end subroutine decsol_option_str
333 :
334 0 : logical function is_block_tridiagonal_decsol(which_decsol_option)
335 : use mtx_def
336 : integer, intent(in) :: which_decsol_option
337 0 : is_block_tridiagonal_decsol = (which_decsol_option == bcyclic_dble)
338 0 : end function is_block_tridiagonal_decsol
339 :
340 : end module mtx_lib
|
