Line data Source code
1 : #if 0
2 : ! This file is part of ELPA.
3 : !
4 : ! The ELPA library was originally created by the ELPA consortium,
5 : ! consisting of the following organizations:
6 : !
7 : ! - Max Planck Computing and Data Facility (MPCDF), formerly known as
8 : ! Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
9 : ! - Bergische Universität Wuppertal, Lehrstuhl für angewandte
10 : ! Informatik,
11 : ! - Technische Universität München, Lehrstuhl für Informatik mit
12 : ! Schwerpunkt Wissenschaftliches Rechnen ,
13 : ! - Fritz-Haber-Institut, Berlin, Abt. Theorie,
14 : ! - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
15 : ! Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
16 : ! and
17 : ! - IBM Deutschland GmbH
18 : !
19 : !
20 : ! More information can be found here:
21 : ! http://elpa.mpcdf.mpg.de/
22 : !
23 : ! ELPA is free software: you can redistribute it and/or modify
24 : ! it under the terms of the version 3 of the license of the
25 : ! GNU Lesser General Public License as published by the Free
26 : ! Software Foundation.
27 : !
28 : ! ELPA is distributed in the hope that it will be useful,
29 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
30 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 : ! GNU Lesser General Public License for more details.
32 : !
33 : ! You should have received a copy of the GNU Lesser General Public License
34 : ! along with ELPA. If not, see <http://www.gnu.org/licenses/>
35 : !
36 : ! ELPA reflects a substantial effort on the part of the original
37 : ! ELPA consortium, and we ask you to respect the spirit of the
38 : ! license that we chose: i.e., please contribute any changes you
39 : ! may have back to the original ELPA library distribution, and keep
40 : ! any derivatives of ELPA under the same license that we chose for
41 : ! the original distribution, the GNU Lesser General Public License.
42 : !
43 : !
44 : ! --------------------------------------------------------------------------------------------------
45 : !
46 : ! This file contains the compute intensive kernels for the Householder transformations.
47 : ! It should be compiled with the highest possible optimization level.
48 : !
49 : ! On Intel use -O3 -xSSE4.2 (or the SSE level fitting to your CPU)
50 : !
51 : ! Copyright of the original code rests with the authors inside the ELPA
52 : ! consortium. The copyright of any additional modifications shall rest
53 : ! with their original authors, but shall adhere to the licensing terms
54 : ! distributed along with the original code in the file "COPYING".
55 : !
56 : ! --------------------------------------------------------------------------------------------------
57 : #endif
58 :
59 : ! the Intel compiler creates a temp array copy of array q!
60 : ! this should be prevented, if possible without using assumed size arrays
61 : #ifdef DOUBLE_PRECISION_COMPLEX
62 1320960 : subroutine single_hh_trafo_complex_generic_double(q, hh, nb, nq, ldq)
63 : #else
64 660480 : subroutine single_hh_trafo_complex_generic_single(q, hh, nb, nq, ldq)
65 : #endif
66 : use precision
67 : use elpa_abstract_impl
68 : implicit none
69 :
70 : ! class(elpa_abstract_impl_t), intent(inout) :: obj
71 : integer(kind=ik), intent(in) :: nb, nq, ldq
72 : #ifdef USE_ASSUMED_SIZE
73 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(ldq,*)
74 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(*)
75 : #else
76 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(1:ldq,1:nb)
77 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(1:nb)
78 : #endif
79 :
80 : integer(kind=ik) :: i
81 :
82 : !#ifdef DOUBLE_PRECISION_COMPLEX
83 : !
84 : ! call obj%timer%start("kernel generic: single_hh_trafo_complex_generic_double")
85 : !
86 : !#else
87 : !
88 : ! call obj%timer%start("kernel generic: single_hh_trafo_complex_generic_single")
89 : !
90 : !#endif
91 : ! Safety only:
92 :
93 1981440 : if(mod(ldq,4) /= 0) STOP 'double_hh_trafo: ldq not divisible by 4!'
94 :
95 : ! Do the Householder transformations
96 :
97 : ! Always a multiple of 4 Q-rows is transformed, even if nq is smaller
98 :
99 7430400 : do i=1,nq-8,12
100 : #ifdef DOUBLE_PRECISION_COMPLEX
101 :
102 : #ifdef USE_ASSUMED_SIZE
103 1816320 : call hh_trafo_complex_kernel_12_double(q(i,1),hh, nb, ldq)
104 : #else
105 1816320 : call hh_trafo_complex_kernel_12_double(q(i:ldq,1:nb),hh(1:nb), nb, ldq)
106 : #endif
107 :
108 : #else
109 :
110 : #ifdef USE_ASSUMED_SIZE
111 908160 : call hh_trafo_complex_kernel_12_single(q(i,1),hh, nb, ldq)
112 : #else
113 908160 : call hh_trafo_complex_kernel_12_single(q(i:ldq,1:nb),hh(1:nb), nb, ldq)
114 : #endif
115 :
116 : #endif
117 : enddo
118 :
119 : ! i > nq-8 now, i.e. at most 8 rows remain
120 :
121 1981440 : if(nq-i+1 > 4) then
122 : #ifdef DOUBLE_PRECISION_COMPLEX
123 :
124 : #ifdef USE_ASSUMED_SIZE
125 165120 : call hh_trafo_complex_kernel_8_double(q(i,1),hh, nb, ldq)
126 : #else
127 165120 : call hh_trafo_complex_kernel_8_double(q(i:ldq,1:nb),hh(1:nb), nb, ldq)
128 : #endif
129 :
130 : #else
131 :
132 : #ifdef USE_ASSUMED_SIZE
133 82560 : call hh_trafo_complex_kernel_8_single(q(i,1),hh, nb, ldq)
134 : #else
135 82560 : call hh_trafo_complex_kernel_8_single(q(i:ldq,1:nb),hh(1:nb), nb, ldq)
136 : #endif
137 :
138 : #endif
139 1486080 : else if(nq-i+1 > 0) then
140 : #ifdef DOUBLE_PRECISION_COMPLEX
141 :
142 : #ifdef USE_ASSUMED_SIZE
143 495360 : call hh_trafo_complex_kernel_4_double(q(i,1),hh, nb, ldq)
144 : #else
145 495360 : call hh_trafo_complex_kernel_4_double(q(i:ldq,1:nb),hh(1:nb), nb, ldq)
146 : #endif
147 :
148 : #else
149 :
150 : #ifdef USE_ASSUMED_SIZE
151 247680 : call hh_trafo_complex_kernel_4_single(q(i,1),hh, nb, ldq)
152 : #else
153 247680 : call hh_trafo_complex_kernel_4_single(q(i:ldq,1:nb),hh(1:nb), nb, ldq)
154 : #endif
155 :
156 : #endif
157 : endif
158 :
159 : !#ifdef DOUBLE_PRECISION_COMPLEX
160 : !
161 : ! call obj%timer%stop("kernel generic: single_hh_trafo_complex_generic_double")
162 : !
163 : !#else
164 : !
165 : ! call obj%timer%stop("kernel generic: single_hh_trafo_complex_generic_single")
166 : !
167 : !#endif
168 :
169 : #ifdef DOUBLE_PRECISION_COMPLEX
170 1320960 : end subroutine single_hh_trafo_complex_generic_double
171 : #else
172 660480 : end subroutine single_hh_trafo_complex_generic_single
173 : #endif
174 : ! --------------------------------------------------------------------------------------------------
175 :
176 : #ifdef DOUBLE_PRECISION_COMPLEX
177 0 : subroutine double_hh_trafo_complex_generic_double(q, hh, nb, nq, ldq, ldh)
178 : #else
179 0 : subroutine double_hh_trafo_complex_generic_single(q, hh, nb, nq, ldq, ldh)
180 : #endif
181 : use precision
182 : implicit none
183 :
184 : integer(kind=ik), intent(in) :: nb, nq, ldq, ldh
185 : #ifdef USE_ASSUMED_SIZE
186 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(ldq,*)
187 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(ldh,*)
188 : #else
189 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(1:ldq,1:nb+1)
190 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(1:ldh,1:2)
191 : #endif
192 : complex(kind=COMPLEX_DATATYPE) :: s
193 :
194 : integer(kind=ik) :: i
195 :
196 : ! Safety only:
197 : #ifdef DOUBLE_PRECISION_COMPLEX
198 :
199 : ! call obj%timer%start("kernel generic: double_hh_trafo_complex_generic_double")
200 :
201 : #else
202 :
203 : ! call obj%timer%start("kernel generic: double_hh_trafo_complex_generic_single")
204 :
205 : #endif
206 0 : if(mod(ldq,4) /= 0) STOP 'double_hh_trafo: ldq not divisible by 4!'
207 :
208 : ! Calculate dot product of the two Householder vectors
209 :
210 0 : s = conjg(hh(2,2)*1)
211 0 : do i=3,nb
212 0 : s = s+(conjg(hh(i,2))*hh(i-1,1))
213 : enddo
214 :
215 : ! Do the Householder transformations
216 :
217 : ! Always a multiple of 4 Q-rows is transformed, even if nq is smaller
218 :
219 0 : do i=1,nq,4
220 : #ifdef DOUBLE_PRECISION_COMPLEX
221 :
222 : #ifdef USE_ASSUMED_SIZE
223 0 : call hh_trafo_complex_kernel_4_2hv_double(q(i,1),hh, nb, ldq, ldh, s)
224 : #else
225 0 : call hh_trafo_complex_kernel_4_2hv_double(q(i:ldq,1:nb+1),hh(1:ldh,1:2), nb, ldq, ldh, s)
226 : #endif
227 :
228 : #else
229 :
230 : #ifdef USE_ASSUMED_SIZE
231 0 : call hh_trafo_complex_kernel_4_2hv_single(q(i,1),hh, nb, ldq, ldh, s)
232 : #else
233 0 : call hh_trafo_complex_kernel_4_2hv_single(q(i:ldq,1:nb+1),hh(1:ldh,1:2), nb, ldq, ldh, s)
234 : #endif
235 :
236 : #endif
237 : enddo
238 :
239 : !do i=1,nq-8,12
240 : #ifdef DOUBLE_PRECISION_COMPLEX
241 :
242 : #ifdef USE_ASSUMED_SIZE
243 : ! call hh_trafo_complex_kernel_12_2hv(q(i,1),hh, nb, ldq, ldh, s)
244 : #else
245 : ! call hh_trafo_complex_kernel_12_2hv(q(i:ldq,1:nb+1),hh(1:ldh,1:2), nb, ldq, ldh, s)
246 : #endif
247 :
248 : #else
249 :
250 : #ifdef USE_ASSUMED_SIZE
251 : ! call hh_trafo_complex_kernel_12_2hv(q(i,1),hh, nb, ldq, ldh, s)
252 : #else
253 : ! call hh_trafo_complex_kernel_12_2hv(q(i:ldq,1:nb+1),hh(1:ldh,1:2), nb, ldq, ldh, s)
254 : #endif
255 :
256 : #endif
257 : !enddo
258 :
259 : ! i > nq-8 now, i.e. at most 8 rows remain
260 :
261 : !if(nq-i+1 > 4) then
262 : #ifdef DOUBLE_PRECISION_COMPLEX
263 :
264 : #ifdef USE_ASSUMED_SIZE
265 : ! call hh_trafo_complex_kernel_8_2hv(q(i,1),hh, nb, ldq, ldh, s)
266 : #else
267 : ! call hh_trafo_complex_kernel_8_2hv(q(i:ldq,1:nb+1),hh(1:ldh,1:2), nb, ldq, ldh, s)
268 : #endif
269 :
270 : #else
271 :
272 : #ifdef USE_ASSUMED_SIZE
273 : ! call hh_trafo_complex_kernel_8_2hv(q(i,1),hh, nb, ldq, ldh, s)
274 : #else
275 : ! call hh_trafo_complex_kernel_8_2hv(q(i:ldq,1:nb+1),hh(1:ldh,1:2), nb, ldq, ldh, s)
276 : #endif
277 :
278 : #endif
279 : !else if(nq-i+1 > 0) then
280 : #ifdef DOUBLE_PRECISION_COMPLEX
281 :
282 : #ifdef USE_ASSUMED_SIZE
283 : ! call hh_trafo_complex_kernel_4_2hv(q(i:ldq,1:nb+1),hh(1:ldh,1:2), nb, ldq, ldh, s)
284 : #else
285 :
286 : #endif
287 :
288 : #else
289 :
290 : #ifdef USE_ASSUMED_SIZE
291 : ! call hh_trafo_complex_kernel_4_2hv(q(i:ldq,1:nb+1),hh(1:ldh,1:2), nb, ldq, ldh, s)
292 : #else
293 :
294 : #endif
295 :
296 : #endif
297 : !endif
298 :
299 : #ifdef DOUBLE_PRECISION_COMPLEX
300 :
301 : ! call obj%timer%stop("kernel generic: double_hh_trafo_complex_generic_double")
302 :
303 : #else
304 :
305 : ! call obj%timer%stop("kernel generic: double_hh_trafo_complex_generic_single")
306 :
307 : #endif
308 :
309 : #ifdef DOUBLE_PRECISION_COMPLEX
310 0 : end subroutine double_hh_trafo_complex_generic_double
311 : #else
312 0 : end subroutine double_hh_trafo_complex_generic_single
313 : #endif
314 : ! --------------------------------------------------------------------------------------------------
315 :
316 : #ifdef DOUBLE_PRECISION_COMPLEX
317 3632640 : subroutine hh_trafo_complex_kernel_12_double(q, hh, nb, ldq)
318 : #else
319 1816320 : subroutine hh_trafo_complex_kernel_12_single(q, hh, nb, ldq)
320 : #endif
321 : use precision
322 : implicit none
323 :
324 : integer(kind=ik), intent(in) :: nb, ldq
325 : #ifdef USE_ASSUMED_SIZE
326 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(ldq,*)
327 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(*)
328 : #else
329 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(:,:)
330 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(1:nb)
331 : #endif
332 : complex(kind=COMPLEX_DATATYPE) :: x1, x2, x3, x4, x5, x6, x7, x8, x9, xa, xb, xc
333 : complex(kind=COMPLEX_DATATYPE) :: h1, tau1
334 : integer(kind=ik) :: i
335 :
336 : #ifdef DOUBLE_PRECISION_COMPLEX
337 :
338 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_12_double")
339 :
340 : #else
341 :
342 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_12_single")
343 :
344 : #endif
345 5448960 : x1 = q(1,1)
346 5448960 : x2 = q(2,1)
347 5448960 : x3 = q(3,1)
348 5448960 : x4 = q(4,1)
349 5448960 : x5 = q(5,1)
350 5448960 : x6 = q(6,1)
351 5448960 : x7 = q(7,1)
352 5448960 : x8 = q(8,1)
353 5448960 : x9 = q(9,1)
354 5448960 : xa = q(10,1)
355 5448960 : xb = q(11,1)
356 5448960 : xc = q(12,1)
357 :
358 : !DEC$ VECTOR ALIGNED
359 174366720 : do i=2,nb
360 168917760 : h1 = conjg(hh(i))
361 168917760 : x1 = x1 + q(1,i)*h1
362 168917760 : x2 = x2 + q(2,i)*h1
363 168917760 : x3 = x3 + q(3,i)*h1
364 168917760 : x4 = x4 + q(4,i)*h1
365 168917760 : x5 = x5 + q(5,i)*h1
366 168917760 : x6 = x6 + q(6,i)*h1
367 168917760 : x7 = x7 + q(7,i)*h1
368 168917760 : x8 = x8 + q(8,i)*h1
369 168917760 : x9 = x9 + q(9,i)*h1
370 168917760 : xa = xa + q(10,i)*h1
371 168917760 : xb = xb + q(11,i)*h1
372 168917760 : xc = xc + q(12,i)*h1
373 : enddo
374 :
375 5448960 : tau1 = hh(1)
376 :
377 5448960 : h1 = -tau1
378 5448960 : x1 = x1*h1
379 5448960 : x2 = x2*h1
380 5448960 : x3 = x3*h1
381 5448960 : x4 = x4*h1
382 5448960 : x5 = x5*h1
383 5448960 : x6 = x6*h1
384 5448960 : x7 = x7*h1
385 5448960 : x8 = x8*h1
386 5448960 : x9 = x9*h1
387 5448960 : xa = xa*h1
388 5448960 : xb = xb*h1
389 5448960 : xc = xc*h1
390 :
391 5448960 : q(1,1) = q(1,1) + x1
392 5448960 : q(2,1) = q(2,1) + x2
393 5448960 : q(3,1) = q(3,1) + x3
394 5448960 : q(4,1) = q(4,1) + x4
395 5448960 : q(5,1) = q(5,1) + x5
396 5448960 : q(6,1) = q(6,1) + x6
397 5448960 : q(7,1) = q(7,1) + x7
398 5448960 : q(8,1) = q(8,1) + x8
399 5448960 : q(9,1) = q(9,1) + x9
400 5448960 : q(10,1) = q(10,1) + xa
401 5448960 : q(11,1) = q(11,1) + xb
402 5448960 : q(12,1) = q(12,1) + xc
403 :
404 : !DEC$ VECTOR ALIGNED
405 174366720 : do i=2,nb
406 168917760 : h1 = hh(i)
407 168917760 : q(1,i) = q(1,i) + x1*h1
408 168917760 : q(2,i) = q(2,i) + x2*h1
409 168917760 : q(3,i) = q(3,i) + x3*h1
410 168917760 : q(4,i) = q(4,i) + x4*h1
411 168917760 : q(5,i) = q(5,i) + x5*h1
412 168917760 : q(6,i) = q(6,i) + x6*h1
413 168917760 : q(7,i) = q(7,i) + x7*h1
414 168917760 : q(8,i) = q(8,i) + x8*h1
415 168917760 : q(9,i) = q(9,i) + x9*h1
416 168917760 : q(10,i) = q(10,i) + xa*h1
417 168917760 : q(11,i) = q(11,i) + xb*h1
418 168917760 : q(12,i) = q(12,i) + xc*h1
419 : enddo
420 :
421 : #ifdef DOUBLE_PRECISION_COMPLEX
422 :
423 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_12_double")
424 :
425 : #else
426 :
427 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_12_single")
428 :
429 : #endif
430 :
431 : #ifdef DOUBLE_PRECISION_COMPLEX
432 3632640 : end subroutine hh_trafo_complex_kernel_12_double
433 : #else
434 1816320 : end subroutine hh_trafo_complex_kernel_12_single
435 : #endif
436 : ! --------------------------------------------------------------------------------------------------
437 :
438 : #ifdef DOUBLE_PRECISION_COMPLEX
439 330240 : subroutine hh_trafo_complex_kernel_8_double(q, hh, nb, ldq)
440 : #else
441 165120 : subroutine hh_trafo_complex_kernel_8_single(q, hh, nb, ldq)
442 : #endif
443 : use precision
444 : implicit none
445 :
446 : integer(kind=ik), intent(in) :: nb, ldq
447 : #ifdef USE_ASSUMED_SIZE
448 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(ldq,*)
449 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(*)
450 : #else
451 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(:,:)
452 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(1:nb)
453 : #endif
454 : complex(kind=COMPLEX_DATATYPE) :: x1, x2, x3, x4, x5, x6, x7, x8
455 : complex(kind=COMPLEX_DATATYPE) :: h1, tau1
456 : integer(kind=ik) :: i
457 :
458 : #ifdef DOUBLE_PRECISION_COMPLEX
459 :
460 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_8_double")
461 :
462 : #else
463 :
464 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_8_single")
465 :
466 : #endif
467 495360 : x1 = q(1,1)
468 495360 : x2 = q(2,1)
469 495360 : x3 = q(3,1)
470 495360 : x4 = q(4,1)
471 495360 : x5 = q(5,1)
472 495360 : x6 = q(6,1)
473 495360 : x7 = q(7,1)
474 495360 : x8 = q(8,1)
475 :
476 : !DEC$ VECTOR ALIGNED
477 15851520 : do i=2,nb
478 15356160 : h1 = conjg(hh(i))
479 15356160 : x1 = x1 + q(1,i)*h1
480 15356160 : x2 = x2 + q(2,i)*h1
481 15356160 : x3 = x3 + q(3,i)*h1
482 15356160 : x4 = x4 + q(4,i)*h1
483 15356160 : x5 = x5 + q(5,i)*h1
484 15356160 : x6 = x6 + q(6,i)*h1
485 15356160 : x7 = x7 + q(7,i)*h1
486 15356160 : x8 = x8 + q(8,i)*h1
487 : enddo
488 :
489 495360 : tau1 = hh(1)
490 :
491 495360 : h1 = -tau1
492 495360 : x1 = x1*h1
493 495360 : x2 = x2*h1
494 495360 : x3 = x3*h1
495 495360 : x4 = x4*h1
496 495360 : x5 = x5*h1
497 495360 : x6 = x6*h1
498 495360 : x7 = x7*h1
499 495360 : x8 = x8*h1
500 :
501 495360 : q(1,1) = q(1,1) + x1
502 495360 : q(2,1) = q(2,1) + x2
503 495360 : q(3,1) = q(3,1) + x3
504 495360 : q(4,1) = q(4,1) + x4
505 495360 : q(5,1) = q(5,1) + x5
506 495360 : q(6,1) = q(6,1) + x6
507 495360 : q(7,1) = q(7,1) + x7
508 495360 : q(8,1) = q(8,1) + x8
509 :
510 : !DEC$ VECTOR ALIGNED
511 15851520 : do i=2,nb
512 15356160 : h1 = hh(i)
513 15356160 : q(1,i) = q(1,i) + x1*h1
514 15356160 : q(2,i) = q(2,i) + x2*h1
515 15356160 : q(3,i) = q(3,i) + x3*h1
516 15356160 : q(4,i) = q(4,i) + x4*h1
517 15356160 : q(5,i) = q(5,i) + x5*h1
518 15356160 : q(6,i) = q(6,i) + x6*h1
519 15356160 : q(7,i) = q(7,i) + x7*h1
520 15356160 : q(8,i) = q(8,i) + x8*h1
521 : enddo
522 :
523 : #ifdef DOUBLE_PRECISION_COMPLEX
524 :
525 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_8_double")
526 :
527 : #else
528 :
529 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_8_single")
530 :
531 : #endif
532 :
533 : #ifdef DOUBLE_PRECISION_COMPLEX
534 330240 : end subroutine hh_trafo_complex_kernel_8_double
535 : #else
536 165120 : end subroutine hh_trafo_complex_kernel_8_single
537 : #endif
538 : ! --------------------------------------------------------------------------------------------------
539 :
540 : #ifdef DOUBLE_PRECISION_COMPLEX
541 990720 : subroutine hh_trafo_complex_kernel_4_double(q, hh, nb, ldq)
542 : #else
543 495360 : subroutine hh_trafo_complex_kernel_4_single(q, hh, nb, ldq)
544 : #endif
545 : use precision
546 : implicit none
547 :
548 : integer(kind=ik), intent(in) :: nb, ldq
549 : #ifdef USE_ASSUMED_SIZE
550 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(ldq,*)
551 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(*)
552 : #else
553 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(:,:)
554 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(1:nb)
555 : #endif
556 : complex(kind=COMPLEX_DATATYPE) :: x1, x2, x3, x4
557 : complex(kind=COMPLEX_DATATYPE) :: h1, tau1
558 : integer(kind=ik) :: i
559 :
560 : #ifdef DOUBLE_PRECISION_COMPLEX
561 :
562 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_4_double")
563 :
564 : #else
565 :
566 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_4_single")
567 :
568 : #endif
569 :
570 1486080 : x1 = q(1,1)
571 1486080 : x2 = q(2,1)
572 1486080 : x3 = q(3,1)
573 1486080 : x4 = q(4,1)
574 :
575 : !DEC$ VECTOR ALIGNED
576 47554560 : do i=2,nb
577 46068480 : h1 = conjg(hh(i))
578 46068480 : x1 = x1 + q(1,i)*h1
579 46068480 : x2 = x2 + q(2,i)*h1
580 46068480 : x3 = x3 + q(3,i)*h1
581 46068480 : x4 = x4 + q(4,i)*h1
582 : enddo
583 :
584 1486080 : tau1 = hh(1)
585 :
586 1486080 : h1 = -tau1
587 1486080 : x1 = x1*h1
588 1486080 : x2 = x2*h1
589 1486080 : x3 = x3*h1
590 1486080 : x4 = x4*h1
591 :
592 1486080 : q(1,1) = q(1,1) + x1
593 1486080 : q(2,1) = q(2,1) + x2
594 1486080 : q(3,1) = q(3,1) + x3
595 1486080 : q(4,1) = q(4,1) + x4
596 :
597 : !DEC$ VECTOR ALIGNED
598 47554560 : do i=2,nb
599 46068480 : h1 = hh(i)
600 46068480 : q(1,i) = q(1,i) + x1*h1
601 46068480 : q(2,i) = q(2,i) + x2*h1
602 46068480 : q(3,i) = q(3,i) + x3*h1
603 46068480 : q(4,i) = q(4,i) + x4*h1
604 : enddo
605 :
606 : #ifdef DOUBLE_PRECISION_COMPLEX
607 :
608 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_4_double")
609 :
610 : #else
611 :
612 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_4_single")
613 :
614 : #endif
615 :
616 : #ifdef DOUBLE_PRECISION_COMPLEX
617 990720 : end subroutine hh_trafo_complex_kernel_4_double
618 : #else
619 495360 : end subroutine hh_trafo_complex_kernel_4_single
620 : #endif
621 : ! --------------------------------------------------------------------------------------------------
622 :
623 : #ifdef DOUBLE_PRECISION_COMPLEX
624 0 : subroutine hh_trafo_complex_kernel_4_2hv_double(q, hh, nb, ldq, ldh, s)
625 : #else
626 0 : subroutine hh_trafo_complex_kernel_4_2hv_single(q, hh, nb, ldq, ldh, s)
627 : #endif
628 : use precision
629 : implicit none
630 :
631 : integer(kind=ik), intent(in) :: nb, ldq, ldh
632 : #ifdef USE_ASSUMED_SIZE
633 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(ldq,*)
634 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(ldh,*)
635 : #else
636 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(:,:)
637 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(1:ldh,1:2)
638 : #endif
639 : complex(kind=COMPLEX_DATATYPE), intent(in) :: s
640 :
641 : complex(kind=COMPLEX_DATATYPE) :: x1, x2, x3, x4, y1, y2, y3, y4
642 : complex(kind=COMPLEX_DATATYPE) :: h1, h2, tau1, tau2
643 : integer(kind=ik) :: i
644 :
645 : #ifdef DOUBLE_PRECISION_COMPLEX
646 :
647 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_4_2hv_double")
648 :
649 : #else
650 :
651 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_4_2hv_single")
652 :
653 : #endif
654 0 : x1 = q(1,2)
655 0 : x2 = q(2,2)
656 0 : x3 = q(3,2)
657 0 : x4 = q(4,2)
658 :
659 0 : y1 = q(1,1) + q(1,2)*conjg(hh(2,2))
660 0 : y2 = q(2,1) + q(2,2)*conjg(hh(2,2))
661 0 : y3 = q(3,1) + q(3,2)*conjg(hh(2,2))
662 0 : y4 = q(4,1) + q(4,2)*conjg(hh(2,2))
663 :
664 : !DEC$ VECTOR ALIGNED
665 0 : do i=3,nb
666 0 : h1 = conjg(hh(i-1,1))
667 0 : h2 = conjg(hh(i,2))
668 0 : x1 = x1 + q(1,i)*h1
669 0 : y1 = y1 + q(1,i)*h2
670 0 : x2 = x2 + q(2,i)*h1
671 0 : y2 = y2 + q(2,i)*h2
672 0 : x3 = x3 + q(3,i)*h1
673 0 : y3 = y3 + q(3,i)*h2
674 0 : x4 = x4 + q(4,i)*h1
675 0 : y4 = y4 + q(4,i)*h2
676 : enddo
677 :
678 0 : x1 = x1 + q(1,nb+1)*conjg(hh(nb,1))
679 0 : x2 = x2 + q(2,nb+1)*conjg(hh(nb,1))
680 0 : x3 = x3 + q(3,nb+1)*conjg(hh(nb,1))
681 0 : x4 = x4 + q(4,nb+1)*conjg(hh(nb,1))
682 :
683 0 : tau1 = hh(1,1)
684 0 : tau2 = hh(1,2)
685 :
686 0 : h1 = -tau1
687 0 : x1 = x1*h1
688 0 : x2 = x2*h1
689 0 : x3 = x3*h1
690 0 : x4 = x4*h1
691 0 : h1 = -tau2
692 0 : h2 = -tau2*s
693 0 : y1 = y1*h1 + x1*h2
694 0 : y2 = y2*h1 + x2*h2
695 0 : y3 = y3*h1 + x3*h2
696 0 : y4 = y4*h1 + x4*h2
697 :
698 0 : q(1,1) = q(1,1) + y1
699 0 : q(2,1) = q(2,1) + y2
700 0 : q(3,1) = q(3,1) + y3
701 0 : q(4,1) = q(4,1) + y4
702 :
703 0 : q(1,2) = q(1,2) + x1 + y1*hh(2,2)
704 0 : q(2,2) = q(2,2) + x2 + y2*hh(2,2)
705 0 : q(3,2) = q(3,2) + x3 + y3*hh(2,2)
706 0 : q(4,2) = q(4,2) + x4 + y4*hh(2,2)
707 :
708 : !DEC$ VECTOR ALIGNED
709 0 : do i=3,nb
710 0 : h1 = hh(i-1,1)
711 0 : h2 = hh(i,2)
712 0 : q(1,i) = q(1,i) + x1*h1 + y1*h2
713 0 : q(2,i) = q(2,i) + x2*h1 + y2*h2
714 0 : q(3,i) = q(3,i) + x3*h1 + y3*h2
715 0 : q(4,i) = q(4,i) + x4*h1 + y4*h2
716 : enddo
717 :
718 0 : q(1,nb+1) = q(1,nb+1) + x1*hh(nb,1)
719 0 : q(2,nb+1) = q(2,nb+1) + x2*hh(nb,1)
720 0 : q(3,nb+1) = q(3,nb+1) + x3*hh(nb,1)
721 0 : q(4,nb+1) = q(4,nb+1) + x4*hh(nb,1)
722 :
723 : #ifdef DOUBLE_PRECISION_COMPLEX
724 :
725 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_4_2hv_double")
726 :
727 : #else
728 :
729 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_4_2hv_single")
730 :
731 : #endif
732 :
733 : #ifdef DOUBLE_PRECISION_COMPLEX
734 0 : end subroutine hh_trafo_complex_kernel_4_2hv_double
735 : #else
736 0 : end subroutine hh_trafo_complex_kernel_4_2hv_single
737 : #endif
738 :
739 : ! --------------------------------------------------------------------------------------------------
740 :
741 : #ifdef DOUBLE_PRECISION_COMPLEX
742 0 : subroutine hh_trafo_complex_kernel_8_2hv_double(q, hh, nb, ldq, ldh, s)
743 : #else
744 0 : subroutine hh_trafo_complex_kernel_8_2hv_single(q, hh, nb, ldq, ldh, s)
745 : #endif
746 : use precision
747 : implicit none
748 :
749 : integer(kind=ik), intent(in) :: nb, ldq, ldh
750 : #ifdef USE_ASSUMED_SIZE
751 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(ldq,*)
752 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(ldh,*)
753 : #else
754 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(:,:)
755 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(1:ldh,1:2)
756 : #endif
757 : complex(kind=COMPLEX_DATATYPE), intent(in) :: s
758 :
759 : complex(kind=COMPLEX_DATATYPE) :: x1, x2, x3, x4, x5, x6 ,x7, x8, y1, y2, y3, y4, y5, y6, y7, y8
760 : complex(kind=COMPLEX_DATATYPE) :: h1, h2, tau1, tau2
761 : integer(kind=ik) :: i
762 :
763 : #ifdef DOUBLE_PRECISION_COMPLEX
764 :
765 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_8_2hv_double")
766 :
767 : #else
768 :
769 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_8_2hv_single")
770 :
771 : #endif
772 :
773 0 : x1 = q(1,2)
774 0 : x2 = q(2,2)
775 0 : x3 = q(3,2)
776 0 : x4 = q(4,2)
777 0 : x5 = q(5,2)
778 0 : x6 = q(6,2)
779 0 : x7 = q(7,2)
780 0 : x8 = q(8,2)
781 :
782 0 : y1 = q(1,1) + q(1,2)*conjg(hh(2,2))
783 0 : y2 = q(2,1) + q(2,2)*conjg(hh(2,2))
784 0 : y3 = q(3,1) + q(3,2)*conjg(hh(2,2))
785 0 : y4 = q(4,1) + q(4,2)*conjg(hh(2,2))
786 0 : y5 = q(5,1) + q(5,2)*conjg(hh(2,2))
787 0 : y6 = q(6,1) + q(6,2)*conjg(hh(2,2))
788 0 : y7 = q(7,1) + q(7,2)*conjg(hh(2,2))
789 0 : y8 = q(8,1) + q(8,2)*conjg(hh(2,2))
790 :
791 : !DEC$ VECTOR ALIGNED
792 0 : do i=3,nb
793 0 : h1 = conjg(hh(i-1,1))
794 0 : h2 = conjg(hh(i,2))
795 0 : x1 = x1 + q(1,i)*h1
796 0 : y1 = y1 + q(1,i)*h2
797 0 : x2 = x2 + q(2,i)*h1
798 0 : y2 = y2 + q(2,i)*h2
799 0 : x3 = x3 + q(3,i)*h1
800 0 : y3 = y3 + q(3,i)*h2
801 0 : x4 = x4 + q(4,i)*h1
802 0 : y4 = y4 + q(4,i)*h2
803 0 : x5 = x5 + q(5,i)*h1
804 0 : y5 = y5 + q(5,i)*h2
805 0 : x6 = x6 + q(6,i)*h1
806 0 : y6 = y6 + q(6,i)*h2
807 0 : x7 = x7 + q(7,i)*h1
808 0 : y7 = y7 + q(7,i)*h2
809 0 : x8 = x8 + q(8,i)*h1
810 0 : y8 = y8 + q(8,i)*h2
811 : enddo
812 :
813 0 : x1 = x1 + q(1,nb+1)*conjg(hh(nb,1))
814 0 : x2 = x2 + q(2,nb+1)*conjg(hh(nb,1))
815 0 : x3 = x3 + q(3,nb+1)*conjg(hh(nb,1))
816 0 : x4 = x4 + q(4,nb+1)*conjg(hh(nb,1))
817 0 : x5 = x5 + q(5,nb+1)*conjg(hh(nb,1))
818 0 : x6 = x6 + q(6,nb+1)*conjg(hh(nb,1))
819 0 : x7 = x7 + q(7,nb+1)*conjg(hh(nb,1))
820 0 : x8 = x8 + q(8,nb+1)*conjg(hh(nb,1))
821 :
822 0 : tau1 = hh(1,1)
823 0 : tau2 = hh(1,2)
824 :
825 0 : h1 = -tau1
826 0 : x1 = x1*h1
827 0 : x2 = x2*h1
828 0 : x3 = x3*h1
829 0 : x4 = x4*h1
830 0 : x5 = x5*h1
831 0 : x6 = x6*h1
832 0 : x7 = x7*h1
833 0 : x8 = x8*h1
834 :
835 0 : h1 = -tau2
836 0 : h2 = -tau2*s
837 0 : y1 = y1*h1 + x1*h2
838 0 : y2 = y2*h1 + x2*h2
839 0 : y3 = y3*h1 + x3*h2
840 0 : y4 = y4*h1 + x4*h2
841 0 : y5 = y5*h1 + x5*h2
842 0 : y6 = y6*h1 + x6*h2
843 0 : y7 = y7*h1 + x7*h2
844 0 : y8 = y8*h1 + x8*h2
845 :
846 0 : q(1,1) = q(1,1) + y1
847 0 : q(2,1) = q(2,1) + y2
848 0 : q(3,1) = q(3,1) + y3
849 0 : q(4,1) = q(4,1) + y4
850 0 : q(5,1) = q(5,1) + y5
851 0 : q(6,1) = q(6,1) + y6
852 0 : q(7,1) = q(7,1) + y7
853 0 : q(8,1) = q(8,1) + y8
854 :
855 0 : q(1,2) = q(1,2) + x1 + y1*hh(2,2)
856 0 : q(2,2) = q(2,2) + x2 + y2*hh(2,2)
857 0 : q(3,2) = q(3,2) + x3 + y3*hh(2,2)
858 0 : q(4,2) = q(4,2) + x4 + y4*hh(2,2)
859 0 : q(5,2) = q(5,2) + x5 + y5*hh(2,2)
860 0 : q(6,2) = q(6,2) + x6 + y6*hh(2,2)
861 0 : q(7,2) = q(7,2) + x7 + y7*hh(2,2)
862 0 : q(8,2) = q(8,2) + x8 + y8*hh(2,2)
863 :
864 : !DEC$ VECTOR ALIGNED
865 0 : do i=3,nb
866 0 : h1 = hh(i-1,1)
867 0 : h2 = hh(i,2)
868 0 : q(1,i) = q(1,i) + x1*h1 + y1*h2
869 0 : q(2,i) = q(2,i) + x2*h1 + y2*h2
870 0 : q(3,i) = q(3,i) + x3*h1 + y3*h2
871 0 : q(4,i) = q(4,i) + x4*h1 + y4*h2
872 0 : q(5,i) = q(5,i) + x5*h1 + y5*h2
873 0 : q(6,i) = q(6,i) + x6*h1 + y6*h2
874 0 : q(7,i) = q(7,i) + x7*h1 + y7*h2
875 0 : q(8,i) = q(8,i) + x8*h1 + y8*h2
876 : enddo
877 :
878 0 : q(1,nb+1) = q(1,nb+1) + x1*hh(nb,1)
879 0 : q(2,nb+1) = q(2,nb+1) + x2*hh(nb,1)
880 0 : q(3,nb+1) = q(3,nb+1) + x3*hh(nb,1)
881 0 : q(4,nb+1) = q(4,nb+1) + x4*hh(nb,1)
882 0 : q(5,nb+1) = q(5,nb+1) + x5*hh(nb,1)
883 0 : q(6,nb+1) = q(6,nb+1) + x6*hh(nb,1)
884 0 : q(7,nb+1) = q(7,nb+1) + x7*hh(nb,1)
885 0 : q(8,nb+1) = q(8,nb+1) + x8*hh(nb,1)
886 :
887 : #ifdef DOUBLE_PRECISION_COMPLEX
888 :
889 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_8_2hv_double")
890 :
891 : #else
892 :
893 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_8_2hv_single")
894 :
895 : #endif
896 :
897 : #ifdef DOUBLE_PRECISION_COMPLEX
898 0 : end subroutine hh_trafo_complex_kernel_8_2hv_double
899 : #else
900 0 : end subroutine hh_trafo_complex_kernel_8_2hv_single
901 : #endif
902 : ! --------------------------------------------------------------------------------------------------
903 :
904 : #ifdef DOUBLE_PRECISION_COMPLEX
905 0 : subroutine hh_trafo_complex_kernel_12_2hv_double(q, hh, nb, ldq, ldh, s)
906 : #else
907 0 : subroutine hh_trafo_complex_kernel_12_2hv_single(q, hh, nb, ldq, ldh, s)
908 : #endif
909 : use precision
910 : implicit none
911 :
912 : integer(kind=ik), intent(in) :: nb, ldq, ldh
913 : #ifdef USE_ASSUMED_SIZE
914 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(ldq,*)
915 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(ldh,*)
916 : #else
917 : complex(kind=COMPLEX_DATATYPE), intent(inout) :: q(:,:)
918 : complex(kind=COMPLEX_DATATYPE), intent(in) :: hh(1:ldh,1:2)
919 : #endif
920 : complex(kind=COMPLEX_DATATYPE), intent(in) :: s
921 :
922 : complex(kind=COMPLEX_DATATYPE) :: x1, x2, x3, x4, x5, x6 ,x7, x8, x9, x10, x11, x12, y1, y2, y3, y4, y5, y6, &
923 : y7, y8, y9, y10, y11, y12
924 : complex(kind=COMPLEX_DATATYPE) :: h1, h2, tau1, tau2
925 : integer(kind=ik) :: i
926 :
927 : #ifdef DOUBLE_PRECISION_COMPLEX
928 :
929 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_12_2hv_double")
930 :
931 : #else
932 :
933 : ! call obj%timer%start("kernel generic: hh_trafo_complex_kernel_12_2hv_single")
934 :
935 : #endif
936 0 : x1 = q(1,2)
937 0 : x2 = q(2,2)
938 0 : x3 = q(3,2)
939 0 : x4 = q(4,2)
940 0 : x5 = q(5,2)
941 0 : x6 = q(6,2)
942 0 : x7 = q(7,2)
943 0 : x8 = q(8,2)
944 0 : x9 = q(9,2)
945 0 : x10 = q(10,2)
946 0 : x11 = q(11,2)
947 0 : x12 = q(12,2)
948 :
949 0 : y1 = q(1,1) + q(1,2)*conjg(hh(2,2))
950 0 : y2 = q(2,1) + q(2,2)*conjg(hh(2,2))
951 0 : y3 = q(3,1) + q(3,2)*conjg(hh(2,2))
952 0 : y4 = q(4,1) + q(4,2)*conjg(hh(2,2))
953 0 : y5 = q(5,1) + q(5,2)*conjg(hh(2,2))
954 0 : y6 = q(6,1) + q(6,2)*conjg(hh(2,2))
955 0 : y7 = q(7,1) + q(7,2)*conjg(hh(2,2))
956 0 : y8 = q(8,1) + q(8,2)*conjg(hh(2,2))
957 0 : y9 = q(9,1) + q(9,2)*conjg(hh(2,2))
958 0 : y10 = q(10,1) + q(10,2)*conjg(hh(2,2))
959 0 : y11 = q(11,1) + q(11,2)*conjg(hh(2,2))
960 0 : y12 = q(12,1) + q(12,2)*conjg(hh(2,2))
961 :
962 : !DEC$ VECTOR ALIGNED
963 0 : do i=3,nb
964 0 : h1 = conjg(hh(i-1,1))
965 0 : h2 = conjg(hh(i,2))
966 0 : x1 = x1 + q(1,i)*h1
967 0 : y1 = y1 + q(1,i)*h2
968 0 : x2 = x2 + q(2,i)*h1
969 0 : y2 = y2 + q(2,i)*h2
970 0 : x3 = x3 + q(3,i)*h1
971 0 : y3 = y3 + q(3,i)*h2
972 0 : x4 = x4 + q(4,i)*h1
973 0 : y4 = y4 + q(4,i)*h2
974 0 : x5 = x5 + q(5,i)*h1
975 0 : y5 = y5 + q(5,i)*h2
976 0 : x6 = x6 + q(6,i)*h1
977 0 : y6 = y6 + q(6,i)*h2
978 0 : x7 = x7 + q(7,i)*h1
979 0 : y7 = y7 + q(7,i)*h2
980 0 : x8 = x8 + q(8,i)*h1
981 0 : y8 = y8 + q(8,i)*h2
982 0 : x9 = x9 + q(9,i)*h1
983 0 : y9 = y9 + q(9,i)*h2
984 0 : x10 = x10 + q(10,i)*h1
985 0 : y10 = y10 + q(10,i)*h2
986 0 : x11 = x11 + q(11,i)*h1
987 0 : y11 = y11 + q(11,i)*h2
988 0 : x12 = x12 + q(12,i)*h1
989 0 : y12 = y12 + q(12,i)*h2
990 : enddo
991 :
992 0 : x1 = x1 + q(1,nb+1)*conjg(hh(nb,1))
993 0 : x2 = x2 + q(2,nb+1)*conjg(hh(nb,1))
994 0 : x3 = x3 + q(3,nb+1)*conjg(hh(nb,1))
995 0 : x4 = x4 + q(4,nb+1)*conjg(hh(nb,1))
996 0 : x5 = x5 + q(5,nb+1)*conjg(hh(nb,1))
997 0 : x6 = x6 + q(6,nb+1)*conjg(hh(nb,1))
998 0 : x7 = x7 + q(7,nb+1)*conjg(hh(nb,1))
999 0 : x8 = x8 + q(8,nb+1)*conjg(hh(nb,1))
1000 0 : x9 = x9 + q(9,nb+1)*conjg(hh(nb,1))
1001 0 : x10 = x10 + q(10,nb+1)*conjg(hh(nb,1))
1002 0 : x11 = x11 + q(11,nb+1)*conjg(hh(nb,1))
1003 0 : x12 = x12 + q(12,nb+1)*conjg(hh(nb,1))
1004 :
1005 0 : tau1 = hh(1,1)
1006 0 : tau2 = hh(1,2)
1007 :
1008 0 : h1 = -tau1
1009 0 : x1 = x1*h1
1010 0 : x2 = x2*h1
1011 0 : x3 = x3*h1
1012 0 : x4 = x4*h1
1013 0 : x5 = x5*h1
1014 0 : x6 = x6*h1
1015 0 : x7 = x7*h1
1016 0 : x8 = x8*h1
1017 0 : x9 = x9*h1
1018 0 : x10 = x10*h1
1019 0 : x11 = x11*h1
1020 0 : x12 = x12*h1
1021 0 : h1 = -tau2
1022 0 : h2 = -tau2*s
1023 0 : y1 = y1*h1 + x1*h2
1024 0 : y2 = y2*h1 + x2*h2
1025 0 : y3 = y3*h1 + x3*h2
1026 0 : y4 = y4*h1 + x4*h2
1027 0 : y5 = y5*h1 + x5*h2
1028 0 : y6 = y6*h1 + x6*h2
1029 0 : y7 = y7*h1 + x7*h2
1030 0 : y8 = y8*h1 + x8*h2
1031 0 : y9 = y9*h1 + x9*h2
1032 0 : y10 = y10*h1 + x10*h2
1033 0 : y11 = y11*h1 + x11*h2
1034 0 : y12 = y12*h1 + x12*h2
1035 :
1036 0 : q(1,1) = q(1,1) + y1
1037 0 : q(2,1) = q(2,1) + y2
1038 0 : q(3,1) = q(3,1) + y3
1039 0 : q(4,1) = q(4,1) + y4
1040 0 : q(5,1) = q(5,1) + y5
1041 0 : q(6,1) = q(6,1) + y6
1042 0 : q(7,1) = q(7,1) + y7
1043 0 : q(8,1) = q(8,1) + y8
1044 0 : q(9,1) = q(9,1) + y9
1045 0 : q(10,1) = q(10,1) + y10
1046 0 : q(11,1) = q(11,1) + y11
1047 0 : q(12,1) = q(12,1) + y12
1048 :
1049 0 : q(1,2) = q(1,2) + x1 + y1*hh(2,2)
1050 0 : q(2,2) = q(2,2) + x2 + y2*hh(2,2)
1051 0 : q(3,2) = q(3,2) + x3 + y3*hh(2,2)
1052 0 : q(4,2) = q(4,2) + x4 + y4*hh(2,2)
1053 0 : q(5,2) = q(5,2) + x5 + y5*hh(2,2)
1054 0 : q(6,2) = q(6,2) + x6 + y6*hh(2,2)
1055 0 : q(7,2) = q(7,2) + x7 + y7*hh(2,2)
1056 0 : q(8,2) = q(8,2) + x8 + y8*hh(2,2)
1057 0 : q(9,2) = q(9,2) + x9 + y9*hh(2,2)
1058 0 : q(10,2) = q(10,2) + x10 + y10*hh(2,2)
1059 0 : q(11,2) = q(11,2) + x11 + y11*hh(2,2)
1060 0 : q(12,2) = q(12,2) + x12 + y12*hh(2,2)
1061 :
1062 : !DEC$ VECTOR ALIGNED
1063 0 : do i=3,nb
1064 0 : h1 = hh(i-1,1)
1065 0 : h2 = hh(i,2)
1066 0 : q(1,i) = q(1,i) + x1*h1 + y1*h2
1067 0 : q(2,i) = q(2,i) + x2*h1 + y2*h2
1068 0 : q(3,i) = q(3,i) + x3*h1 + y3*h2
1069 0 : q(4,i) = q(4,i) + x4*h1 + y4*h2
1070 0 : q(5,i) = q(5,i) + x5*h1 + y5*h2
1071 0 : q(6,i) = q(6,i) + x6*h1 + y6*h2
1072 0 : q(7,i) = q(7,i) + x7*h1 + y7*h2
1073 0 : q(8,i) = q(8,i) + x8*h1 + y8*h2
1074 0 : q(9,i) = q(9,i) + x9*h1 + y9*h2
1075 0 : q(10,i) = q(10,i) + x10*h1 + y10*h2
1076 0 : q(11,i) = q(11,i) + x11*h1 + y11*h2
1077 0 : q(12,i) = q(12,i) + x12*h1 + y12*h2
1078 : enddo
1079 :
1080 0 : q(1,nb+1) = q(1,nb+1) + x1*hh(nb,1)
1081 0 : q(2,nb+1) = q(2,nb+1) + x2*hh(nb,1)
1082 0 : q(3,nb+1) = q(3,nb+1) + x3*hh(nb,1)
1083 0 : q(4,nb+1) = q(4,nb+1) + x4*hh(nb,1)
1084 0 : q(5,nb+1) = q(5,nb+1) + x5*hh(nb,1)
1085 0 : q(6,nb+1) = q(6,nb+1) + x6*hh(nb,1)
1086 0 : q(7,nb+1) = q(7,nb+1) + x7*hh(nb,1)
1087 0 : q(8,nb+1) = q(8,nb+1) + x8*hh(nb,1)
1088 0 : q(9,nb+1) = q(9,nb+1) + x9*hh(nb,1)
1089 0 : q(10,nb+1) = q(10,nb+1) + x10*hh(nb,1)
1090 0 : q(11,nb+1) = q(11,nb+1) + x11*hh(nb,1)
1091 0 : q(12,nb+1) = q(12,nb+1) + x12*hh(nb,1)
1092 :
1093 : #ifdef DOUBLE_PRECISION_COMPLEX
1094 :
1095 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_12_2hv_double")
1096 :
1097 : #else
1098 :
1099 : ! call obj%timer%stop("kernel generic: hh_trafo_complex_kernel_12_2hv_single")
1100 :
1101 : #endif
1102 :
1103 : #ifdef DOUBLE_PRECISION_COMPLEX
1104 0 : end subroutine hh_trafo_complex_kernel_12_2hv_double
1105 : #else
1106 0 : end subroutine hh_trafo_complex_kernel_12_2hv_single
1107 : #endif
|
