Line data Source code
1 : module pre_conditioners
2 :
3 : use const_def, only: dp
4 :
5 : implicit none
6 :
7 : private
8 : public :: compute_block_preconditioner, compute_band_preconditioner, compute_block_left_preconditioner
9 :
10 : contains
11 :
12 : !> Computes a pre-conditioning array for rows of a block triadiagonal matrix.
13 : !! Each row in the matrix can then be divided by the corresponding element in this array
14 : !! to form a usually better-conditioned matrix.
15 : !! This element is chosen to be the absolute value of the element in that row with the greatest magnitude.
16 : !!
17 : !! @params ublk The upper blocks. Shape is (nvars, nvars, nblocks-1).
18 : !! @params lblk The lower blocks. Shape is (nvars, nvars, nblocks-1).
19 : !! @params dblk The diagonal blocks. Shape is (nvars, nvars, nblocks).
20 : !! @params nblocks The number of blocks on the diagonal.
21 : !! @params nvars The width and height of each block (e.g. each block is nvars by nvars).
22 : !! @params pre_conditioner Pre-conditioning vector. Computed in this routine.
23 0 : subroutine compute_block_preconditioner(ublk, lblk, dblk, nblocks, nvar, pre_conditioner)
24 : real(dp), dimension(:,:,:), intent(in) :: ublk, lblk, dblk
25 : integer, intent(in) :: nblocks, nvar
26 :
27 : integer :: j, k
28 :
29 : real(dp), dimension(:,:), intent(out) :: pre_conditioner
30 :
31 0 : pre_conditioner = 0d0
32 0 : do j=1,nblocks
33 0 : do k=1,nvar
34 0 : pre_conditioner(k,j) = pre_conditioner(k,j) + sum(abs(dblk(k,:,j)))
35 :
36 0 : if (j < nblocks) then
37 0 : pre_conditioner(k,j) = pre_conditioner(k,j) + sum(abs(ublk(k,:,j)))
38 : end if
39 :
40 0 : if (j > 1) then
41 0 : pre_conditioner(k,j) = pre_conditioner(k,j) + sum(abs(lblk(k,:,j-1)))
42 : end if
43 : end do
44 : end do
45 :
46 0 : end subroutine compute_block_preconditioner
47 :
48 : !> Computes a pre-conditioning array for columns of a block triadiagonal matrix.
49 : !! Each column in the matrix can then be divided by the corresponding element in this array
50 : !! to form a usually better-conditioned matrix.
51 : !! This element is chosen to be the absolute value of the element in that column with the greatest magnitude.
52 : !! Note that the final answer to the linear solve must be divided element-wise by this pre-conditioner.
53 : !!
54 : !! @params ublk The upper blocks. Shape is (nvars, nvars, nblocks-1).
55 : !! @params lblk The lower blocks. Shape is (nvars, nvars, nblocks-1).
56 : !! @params dblk The diagonal blocks. Shape is (nvars, nvars, nblocks).
57 : !! @params nblocks The number of blocks on the diagonal.
58 : !! @params nvars The width and height of each block (e.g. each block is nvars by nvars).
59 : !! @params pre_conditioner Pre-conditioning vector. Computed in this routine.
60 0 : subroutine compute_block_left_preconditioner(ublk, lblk, dblk, nblocks, nvar, &
61 0 : right_pre_conditioner, pre_conditioner)
62 : real(dp), dimension(:,:,:), intent(in) :: ublk, lblk, dblk
63 : real(dp), dimension(:,:), intent(in) :: right_pre_conditioner
64 : integer, intent(in) :: nblocks, nvar
65 :
66 : integer :: j, k, k2
67 :
68 : real(dp), dimension(:,:), intent(out) :: pre_conditioner
69 :
70 0 : pre_conditioner = 0d0
71 0 : do j=1,nblocks
72 0 : do k=1,nvar
73 0 : do k2=1,nvar
74 0 : pre_conditioner(k,j) = pre_conditioner(k,j) + abs(dblk(k2,k,j) / right_pre_conditioner(k2,j))
75 :
76 0 : if (j < nblocks) then
77 0 : pre_conditioner(k,j) = pre_conditioner(k,j) + abs(lblk(k2,k,j) / right_pre_conditioner(k2,j+1))
78 : end if
79 :
80 0 : if (j > 1) then
81 0 : pre_conditioner(k,j) = pre_conditioner(k,j) + abs(ublk(k2,k,j-1) / right_pre_conditioner(k2,j-1))
82 : end if
83 : end do
84 : end do
85 : end do
86 :
87 0 : end subroutine compute_block_left_preconditioner
88 :
89 : !> Computes a pre-conditioning array for rows of a banded matrix.
90 : !! Each row in the matrix can then be divided by the corresponding element in this array
91 : !! to form a usually better-conditioned matrix.
92 : !! This element is chosen to be the absolute value of the element in that row with the greatest magnitude.
93 : !!
94 : !! @params n_upper_bands The number of bands above the diagonal.
95 : !! @params n_lower_bands The number of bands below the diagonal.
96 : !! @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.
97 : !! @params pre_conditioner Pre-conditioning vector. Computed in this routine.
98 0 : subroutine compute_band_preconditioner(matrix_size, n_upper_bands, n_lower_bands, bands, pre_conditioner)
99 : ! Inputs
100 : real(dp), dimension(:,:), intent(in) :: bands
101 : integer, intent(in) :: matrix_size, n_upper_bands, n_lower_bands
102 :
103 : ! Intermediates
104 : integer :: i, j, n_bands
105 :
106 : ! Outputs
107 : real(dp), dimension(:), intent(out) :: pre_conditioner
108 :
109 0 : n_bands = n_upper_bands + n_lower_bands + 1
110 :
111 0 : pre_conditioner = 0d0
112 :
113 : ! Upper bands
114 0 : do j=1,n_upper_bands
115 : ! The first upper band has index 1 in bands and runs 1...matrix_size-n_upper_bands.
116 0 : do i=1,matrix_size + j - n_upper_bands - 1
117 0 : pre_conditioner(i) = max(pre_conditioner(i), abs(bands(j,i)))
118 : end do
119 : end do
120 :
121 : ! Diagonal band
122 0 : do i=1,matrix_size
123 0 : pre_conditioner(i) = max(pre_conditioner(i), abs(bands(n_upper_bands + 1, i)))
124 : end do
125 :
126 : ! Lower bands
127 0 : do j=1,n_lower_bands
128 : ! The first lower band has index n_upper_bands+2 in bands and runs 1...matrix_size-1
129 0 : do i=j+1,matrix_size
130 0 : pre_conditioner(i) = max(pre_conditioner(i), abs(bands(j + n_upper_bands + 1, i-j)))
131 : end do
132 : end do
133 :
134 0 : end subroutine compute_band_preconditioner
135 :
136 : end module pre_conditioners
|
