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