forked from ESCOMP/CLUBB_CESM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
gmres_cache.F90
179 lines (137 loc) · 6.67 KB
/
gmres_cache.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
!----------------------------------------------------------------------------
! $Id$
!==============================================================================
module gmres_cache
#ifdef MKL
use clubb_precision, only: &
dp ! double precision
! Description:
! This module contains cache data structures for the GMRES wrapper class.
!
! This is mostly to allow us to get around some...odd errors when it was
! integrated into the gmres_wrap module. The cache variables are public, as
! they will need to be passed in whenever gmres_solve is called.
implicit none
public :: gmres_cache_matrix_init, gmres_cache_soln, &
gmres_cache_temp_init
private ! Default scope
real( kind = dp ), public, allocatable, dimension(:,:) :: &
gmres_prev_soln, & ! Stores the previous solution vectors from earlier
! GMRES solve runs. The first dimension is for the
! actual vector; the second dimension is to determine
! which cache to access--this is done via the GMRES
! indices for each of the different matrices.
gmres_prev_precond_a ! Stores the previous preconditioner matrix from
! earlier GMRES solve runs. The first dimension is
! for the a-array itself; the second dimension is to
! determine which cached array to access--this is
! done via the GMRES indices for each of the
! different matrices.
!$omp threadprivate( gmres_prev_soln, gmres_prev_precond_a )
real( kind = dp ), public, allocatable, dimension(:) :: &
gmres_temp_intlc, & ! Temporary array that stores GMRES internal values
! for the interlaced matrices (2 x gr%nz grid
! levels)
gmres_temp_norm ! Temporary array that stores GMRES internal values
! for the non-interlaced matrices (gr%nz grid
! levels)
!$omp threadprivate( gmres_tmp_intlc, gmres_temp_norm )
integer, public :: &
gmres_tempsize_norm, & ! Size of the temporary array for
! non-interlaced matrices
gmres_tempsize_intlc ! Size of the temporary array for
! interlaced matrices
!$omp threadprivate( gmres_tempsize_norm, gmres_tempsize_intlc )
integer, public, parameter :: &
maximum_gmres_idx = 1 ! Maximum number of different types of solves the
! wrapper can keep memory for. If new matrices are
! added that GMRES is to be used for, increase this
! number and add a public parameter corresponding to
! the matrix below:
integer, public, parameter :: &
gmres_idx_wp2wp3 = 1 ! GMRES wrapper index for the wp2_wp3 matrices
logical, public, dimension(maximum_gmres_idx) :: &
l_gmres_soln_ok ! Stores if the current solution is "okay"--that is, if an
! initial solution has been passed in for that particular
! cache index. This defaults to false and is set to true
! when a solution is updated.
!$omp threadprivate(l_gmres_soln_ok)
contains
subroutine gmres_cache_temp_init(numeqns) ! Intent(in)
! Description:
! Initialization subroutine for the temporary arrays for GMRES
!
! This subroutine initializes the temporary arrays that are used to work
! the GMRES solver.
!
! These temporary arrays are used for all GMRES solves.
!
! References:
! None
implicit none
! Input Variables
integer, intent(in) :: &
numeqns ! Number of equations for non-interlaced matrices (gr%nz)
integer :: &
numeqns_intlc ! Number of equations for interlaced matrices
numeqns_intlc = numeqns * 2
! Figure out the sizes of the temporary arrays
! The equations were lifted from the Intel documentation of dfgmres:
! http://www.intel.com/software/products/mkl/docs/webhelp/ssr/functn_rci_dfgmres.html
! All of the ipar(15)s have been replaced with "numeqns", as the code
! examples seemed to use N (numeqns) in place of ipar(15).
gmres_tempsize_norm = ((((2*numeqns + 1)*numeqns) &
+ (numeqns*(numeqns+9))/2) + 1) ! Known magic number
gmres_tempsize_intlc = ((((2*numeqns_intlc + 1)*numeqns_intlc) &
+ (numeqns_intlc*(numeqns_intlc+9))/2) + 1) ! Known magic number
! Allocate the temporary arrays
allocate( gmres_temp_intlc(1:gmres_tempsize_intlc), &
gmres_temp_norm(1:gmres_tempsize_norm) )
end subroutine gmres_cache_temp_init
subroutine gmres_cache_matrix_init(max_numeqns, max_elements, & ! Intent(in)
max_gmres_idx) ! Intent(in)
! Description:
! Initialization subroutine for the caches for GMRES.
!
! This initializes the cache that stores the previous solution and
! previous preconditioner values for all GMRES solves.
!
! References:
! None
implicit none
! Input Variables
integer, intent(in) :: &
max_numeqns, & ! Maximum number of equations for a matrix that will be
! solved with GMRES
max_elements, & ! Maximum number of non-zero elements for a matrix that
! will be solved with GMRES
max_gmres_idx ! Maximum number of distinct matrices that will be solved
! with GMRES
allocate( gmres_prev_soln(1:max_numeqns,1:max_gmres_idx), &
gmres_prev_precond_a(1:max_elements,1:max_gmres_idx) )
l_gmres_soln_ok = .false.
end subroutine gmres_cache_matrix_init
subroutine gmres_cache_soln(numeqns, gmres_idx, solution) ! Intent(in)
! Description:
! Subroutine that caches a previous solution for a particular GMRES-solved
! matrix.
!
! Stores the current solution in the cache so it can be referenced for
! the next GMRES solve. This subroutine will also set the solution_ok
! flag for that particular GMRES index.
!
! References:
! None
use clubb_precision, only: &
core_rknd ! Variable(s)
implicit none
integer, intent(in) :: &
numeqns, & ! The number of equations in the solution vector
gmres_idx ! The index for the particular matrix solved by GMRES
real( kind = core_rknd ), dimension(numeqns), intent(in) :: &
solution ! The solution vector to be cached
gmres_prev_soln(1:numeqns,gmres_idx) = solution
l_gmres_soln_ok(gmres_idx) = .true.
end subroutine gmres_cache_soln
#endif /* MKL */
end module gmres_cache