Line data Source code
1 : ! This file is part of ELPA.
2 : !
3 : ! The ELPA library was originally created by the ELPA consortium,
4 : ! consisting of the following organizations:
5 : !
6 : ! - Max Planck Computing and Data Facility (MPCDF), formerly known as
7 : ! Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
8 : ! - Bergische Universität Wuppertal, Lehrstuhl für angewandte
9 : ! Informatik,
10 : ! - Technische Universität München, Lehrstuhl für Informatik mit
11 : ! Schwerpunkt Wissenschaftliches Rechnen ,
12 : ! - Fritz-Haber-Institut, Berlin, Abt. Theorie,
13 : ! - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
14 : ! Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
15 : ! and
16 : ! - IBM Deutschland GmbH
17 : !
18 : !
19 : ! More information can be found here:
20 : ! http://elpa.mpcdf.mpg.de/
21 : !
22 : ! ELPA is free software: you can redistribute it and/or modify
23 : ! it under the terms of the version 3 of the license of the
24 : ! GNU Lesser General Public License as published by the Free
25 : ! Software Foundation.
26 : !
27 : ! ELPA is distributed in the hope that it will be useful,
28 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : ! GNU Lesser General Public License for more details.
31 : !
32 : ! You should have received a copy of the GNU Lesser General Public License
33 : ! along with ELPA. If not, see <http://www.gnu.org/licenses/>
34 : !
35 : ! ELPA reflects a substantial effort on the part of the original
36 : ! ELPA consortium, and we ask you to respect the spirit of the
37 : ! license that we chose: i.e., please contribute any changes you
38 : ! may have back to the original ELPA library distribution, and keep
39 : ! any derivatives of ELPA under the same license that we chose for
40 : ! the original distribution, the GNU Lesser General Public License.
41 : !
42 : ! Author: A. Marek, MPCDF
43 :
44 : subroutine prepare_matrix_random_&
45 : &MATH_DATATYPE&
46 : &_&
47 10656 : &PRECISION&
48 10656 : & (na, myid, sc_desc, a, z, as)
49 :
50 :
51 : use test_util
52 : implicit none
53 : #include "../../src/general/precision_kinds.F90"
54 : integer(kind=ik), intent(in) :: myid, na, sc_desc(:)
55 : MATH_DATATYPE(kind=rck), intent(inout) :: z(:,:), a(:,:), as(:,:)
56 :
57 : #if COMPLEXCASE == 1
58 10368 : real(kind=rk) :: xr(size(a,dim=1), size(a,dim=2))
59 : #endif /* COMPLEXCASE */
60 :
61 :
62 10656 : integer, allocatable :: iseed(:)
63 : integer :: n
64 :
65 : ! for getting a hermitian test matrix A we get a random matrix Z
66 : ! and calculate A = Z + Z**H
67 :
68 : ! we want different random numbers on every process
69 : ! (otherwise A might get rank deficient):
70 :
71 10656 : call random_seed(size=n)
72 10656 : allocate(iseed(n))
73 10656 : iseed(:) = myid
74 10656 : call random_seed(put=iseed)
75 : #if REALCASE == 1
76 5472 : call random_number(z)
77 :
78 5472 : a(:,:) = z(:,:)
79 : #endif /* REALCASE */
80 :
81 : #if COMPLEXCASE == 1
82 5184 : call random_number(xr)
83 :
84 5184 : z(:,:) = xr(:,:)
85 5184 : call RANDOM_NUMBER(xr)
86 5184 : z(:,:) = z(:,:) + (0.0_rk,1.0_rk)*xr(:,:)
87 5184 : a(:,:) = z(:,:)
88 : #endif /* COMPLEXCASE */
89 :
90 10656 : if (myid == 0) then
91 7104 : print '(a)','| Random matrix block has been set up. (only processor 0 confirms this step)'
92 : endif
93 :
94 : #if REALCASE == 1
95 : #ifdef WITH_MPI
96 : call p&
97 : &BLAS_CHAR&
98 3648 : &tran(na, na, ONE, z, 1, 1, sc_desc, ONE, a, 1, 1, sc_desc) ! A = A + Z**T
99 : #else /* WITH_MPI */
100 1824 : a = a + transpose(z)
101 : #endif /* WITH_MPI */
102 : #endif /* REALCASE */
103 :
104 : #if COMPLEXCASE == 1
105 : #ifdef WITH_MPI
106 : call p&
107 : &BLAS_CHAR&
108 3456 : &tranc(na, na, ONE, z, 1, 1, sc_desc, ONE, a, 1, 1, sc_desc) ! A = A + Z**H
109 : #else /* WITH_MPI */
110 1728 : a = a + transpose(conjg(z))
111 : #endif /* WITH_MPI */
112 : #endif /* COMPLEXCASE */
113 :
114 :
115 10656 : if (myid == 0) then
116 7104 : print '(a)','| Random matrix block has been symmetrized'
117 : endif
118 :
119 : ! save original matrix A for later accuracy checks
120 :
121 10656 : as = a
122 :
123 10656 : deallocate(iseed)
124 :
125 10656 : end subroutine
126 :
127 : #if REALCASE == 1
128 : #ifdef DOUBLE_PRECISION_REAL
129 : !c> void prepare_matrix_random_real_double_f(int na, int myid, int na_rows, int na_cols,
130 : !c> int sc_desc[9],
131 : !c> double *a, double *z, double *as);
132 : #else
133 : !c> void prepare_matrix_random_real_single_f(int na, int myid, int na_rows, int na_cols,
134 : !c> int sc_desc[9],
135 : !c> float *a, float *z, float *as);
136 : #endif
137 : #endif /* REALCASE */
138 :
139 : #if COMPLEXCASE == 1
140 : #ifdef DOUBLE_PRECISION_COMPLEX
141 : !c> void prepare_matrix_random_complex_double_f(int na, int myid, int na_rows, int na_cols,
142 : !c> int sc_desc[9],
143 : !c> complex double *a, complex double *z, complex double *as);
144 : #else
145 : !c> void prepare_matrix_random_complex_single_f(int na, int myid, int na_rows, int na_cols,
146 : !c> int sc_desc[9],
147 : !c> complex float *a, complex float *z, complex float *as);
148 : #endif
149 : #endif /* COMPLEXCASE */
150 :
151 : subroutine prepare_matrix_random_&
152 : &MATH_DATATYPE&
153 : &_wrapper_&
154 3360 : &PRECISION&
155 3360 : & (na, myid, na_rows, na_cols, sc_desc, a, z, as) &
156 : bind(C, name="prepare_matrix_random_&
157 : &MATH_DATATYPE&
158 : &_&
159 : &PRECISION&
160 : &_f")
161 : use iso_c_binding
162 :
163 : implicit none
164 : #include "../../src/general/precision_kinds.F90"
165 :
166 : integer(kind=c_int) , value :: myid, na, na_rows, na_cols
167 : integer(kind=c_int) :: sc_desc(1:9)
168 : MATH_DATATYPE(kind=rck) :: z(1:na_rows,1:na_cols), a(1:na_rows,1:na_cols), &
169 : as(1:na_rows,1:na_cols)
170 : call prepare_matrix_random_&
171 : &MATH_DATATYPE&
172 : &_&
173 : &PRECISION&
174 3360 : & (na, myid, sc_desc, a, z, as)
175 3360 : end subroutine
176 :
177 :
178 : subroutine prepare_matrix_toeplitz_&
179 : &MATH_DATATYPE&
180 : &_&
181 2304 : &PRECISION&
182 2304 : & (na, diagonalElement, subdiagonalElement, d, sd, ds, sds, a, as, &
183 : nblk, np_rows, np_cols, my_prow, my_pcol)
184 : use test_util
185 : implicit none
186 : #include "../../src/general/precision_kinds.F90"
187 :
188 : integer, intent(in) :: na, nblk, np_rows, np_cols, my_prow, my_pcol
189 : MATH_DATATYPE(kind=rck) :: diagonalElement, subdiagonalElement
190 : MATH_DATATYPE(kind=rck) :: d(:), sd(:), ds(:), sds(:)
191 : MATH_DATATYPE(kind=rck) :: a(:,:), as(:,:)
192 :
193 : integer :: ii, rowLocal, colLocal
194 :
195 2304 : d(:) = diagonalElement
196 2304 : sd(:) = subdiagonalElement
197 2304 : a(:,:) = ZERO
198 :
199 : ! set up the diagonal and subdiagonals (for general solver test)
200 347904 : do ii=1, na ! for diagonal elements
201 345600 : if (map_global_array_index_to_local_index(ii, ii, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
202 230400 : a(rowLocal,colLocal) = diagonalElement
203 : endif
204 : enddo
205 345600 : do ii=1, na-1
206 343296 : if (map_global_array_index_to_local_index(ii, ii+1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
207 228864 : a(rowLocal,colLocal) = subdiagonalElement
208 : endif
209 : enddo
210 :
211 345600 : do ii=2, na
212 343296 : if (map_global_array_index_to_local_index(ii, ii-1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
213 228864 : a(rowLocal,colLocal) = subdiagonalElement
214 : endif
215 : enddo
216 :
217 2304 : ds = d
218 2304 : sds = sd
219 2304 : as = a
220 2304 : end subroutine
221 :
222 : subroutine prepare_matrix_toeplitz_mixed_complex&
223 : &_&
224 : &MATH_DATATYPE&
225 : &_&
226 1440 : &PRECISION&
227 : #if COMPLEXCASE == 1
228 1440 : & (na, diagonalElement, subdiagonalElement, d, sd, ds, sds, a, as, &
229 : nblk, np_rows, np_cols, my_prow, my_pcol)
230 : #endif
231 : #if REALCASE == 1
232 : & (na, diagonalElement, subdiagonalElement, d, sd, ds, sds, &
233 : nblk, np_rows, np_cols, my_prow, my_pcol)
234 : #endif
235 : use test_util
236 : implicit none
237 :
238 : integer, intent(in) :: na, nblk, np_rows, np_cols, my_prow, my_pcol
239 : real(kind=C_DATATYPE_KIND) :: diagonalElement, subdiagonalElement
240 :
241 : real(kind=C_DATATYPE_KIND) :: d(:), sd(:), ds(:), sds(:)
242 :
243 : #if COMPLEXCASE == 1
244 : complex(kind=C_DATATYPE_KIND) :: a(:,:), as(:,:)
245 : #endif
246 : #if REALCASE == 1
247 : #endif
248 :
249 : integer :: ii, rowLocal, colLocal
250 : #if COMPLEXCASE == 1
251 1440 : d(:) = diagonalElement
252 1440 : sd(:) = subdiagonalElement
253 :
254 : ! set up the diagonal and subdiagonals (for general solver test)
255 217440 : do ii=1, na ! for diagonal elements
256 216000 : if (map_global_array_index_to_local_index(ii, ii, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
257 144000 : a(rowLocal,colLocal) = diagonalElement
258 : endif
259 : enddo
260 216000 : do ii=1, na-1
261 214560 : if (map_global_array_index_to_local_index(ii, ii+1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
262 143040 : a(rowLocal,colLocal) = subdiagonalElement
263 : endif
264 : enddo
265 :
266 216000 : do ii=2, na
267 214560 : if (map_global_array_index_to_local_index(ii, ii-1, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
268 143040 : a(rowLocal,colLocal) = subdiagonalElement
269 : endif
270 : enddo
271 :
272 1440 : ds = d
273 1440 : sds = sd
274 1440 : as = a
275 : #endif
276 1440 : end subroutine
277 :
278 : subroutine prepare_matrix_frank_&
279 : &MATH_DATATYPE&
280 : &_&
281 1152 : &PRECISION&
282 1152 : & (na, a, z, as, nblk, np_rows, np_cols, my_prow, my_pcol)
283 : use test_util
284 : implicit none
285 :
286 : integer, intent(in) :: na, nblk, np_rows, np_cols, my_prow, my_pcol
287 :
288 : #if REALCASE == 1
289 : real(kind=C_DATATYPE_KIND) :: a(:,:), z(:,:), as(:,:)
290 : #endif
291 : #if COMPLEXCASE == 1
292 : complex(kind=C_DATATYPE_KIND) :: a(:,:), z(:,:), as(:,:)
293 : #endif
294 :
295 : integer :: i, j, rowLocal, colLocal
296 :
297 173952 : do i = 1, na
298 26092800 : do j = 1, na
299 25920000 : if (map_global_array_index_to_local_index(i, j, rowLocal, colLocal, nblk, np_rows, np_cols, my_prow, my_pcol)) then
300 17280000 : if (j .le. i) then
301 8697600 : a(rowLocal,colLocal) = real((na+1-i), kind=C_DATATYPE_KIND) / real(na, kind=C_DATATYPE_KIND)
302 : else
303 8582400 : a(rowLocal,colLocal) = real((na+1-j), kind=C_DATATYPE_KIND) / real(na, kind=C_DATATYPE_KIND)
304 : endif
305 : endif
306 : enddo
307 : enddo
308 :
309 1152 : z(:,:) = a(:,:)
310 1152 : as(:,:) = a(:,:)
311 :
312 1152 : end subroutine
313 :
314 :
315 : ! vim: syntax=fortran
|
