120 SUBROUTINE simpsonrule_init(sr_env, xnodes, nnodes, a, b, shape_id, conv, weights, tnodes_restart)
122 INTEGER,
INTENT(inout) :: nnodes
123 COMPLEX(kind=dp),
DIMENSION(nnodes),
INTENT(out) :: xnodes
124 COMPLEX(kind=dp),
INTENT(in) :: a, b
125 INTEGER,
INTENT(in) :: shape_id
126 REAL(kind=
dp),
INTENT(in) :: conv
128 REAL(kind=
dp),
DIMENSION(nnodes),
INTENT(in), &
129 OPTIONAL :: tnodes_restart
131 CHARACTER(len=*),
PARAMETER :: routinen =
'simpsonrule_init'
133 INTEGER :: handle, icol, irow, ncols, nrows
134 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :), &
135 POINTER :: w_data, w_data_my
138 CALL timeset(routinen, handle)
143 nnodes = 4*((nnodes - 1)/4) + 1
145 sr_env%shape_id = shape_id
149 sr_env%error = huge(0.0_dp)
150 sr_env%error_conv = 0.0_dp
152 NULLIFY (sr_env%error_fm, sr_env%weights)
153 CALL cp_fm_get_info(weights, local_data=w_data, nrow_local=nrows, ncol_local=ncols, matrix_struct=fm_struct)
154 ALLOCATE (sr_env%error_fm, sr_env%weights)
162 w_data_my(irow, icol) = abs(w_data(irow, icol))/15.0_dp
166 NULLIFY (sr_env%integral, sr_env%integral_conv)
167 NULLIFY (sr_env%integral_abc, sr_env%integral_cde, sr_env%integral_ace)
169 ALLOCATE (sr_env%tnodes(nnodes))
171 IF (
PRESENT(tnodes_restart))
THEN
172 sr_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes)
178 CALL timestop(handle)
341 TYPE(
cp_cfm_type),
DIMENSION(:),
INTENT(inout) :: zdata_next
343 CHARACTER(len=*),
PARAMETER :: routinen =
'simpsonrule_refine_integral'
346 COMPLEX(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: zscale
347 COMPLEX(kind=dp),
CONTIGUOUS,
DIMENSION(:, :), &
348 POINTER :: error_zdata
349 INTEGER :: handle, interval, ipoint, jpoint, &
350 nintervals, nintervals_exist, npoints
351 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: inds
352 LOGICAL :: interval_converged, interval_exists
353 REAL(kind=
dp) :: my_bound, rscale
354 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: errors
355 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :), &
356 POINTER :: error_rdata
358 TYPE(simpsonrule_subinterval_type),
ALLOCATABLE, &
359 DIMENSION(:) :: subintervals
361 CALL timeset(routinen, handle)
363 npoints =
SIZE(zdata_next)
364 IF (
ASSOCIATED(sr_env%integral))
THEN
368 cpassert(npoints > 0 .AND. mod(npoints, 4) == 0)
374 cpassert(npoints > 1 .AND. mod(npoints, 4) == 1)
378 nintervals_exist =
SIZE(sr_env%tnodes)
379 cpassert(nintervals_exist >= npoints)
380 ALLOCATE (zscale(npoints))
383 sr_env%a, sr_env%b, sr_env%shape_id, weights=zscale)
386 DO ipoint = 1, npoints
393 nintervals = npoints/4
394 IF (
ASSOCIATED(sr_env%integral))
THEN
396 nintervals_exist =
SIZE(sr_env%subintervals)
397 cpassert(nintervals <= nintervals_exist)
399 ALLOCATE (subintervals(nintervals_exist + nintervals))
401 DO interval = 1, nintervals
402 subintervals(2*interval - 1)%lb = sr_env%subintervals(interval)%lb
403 subintervals(2*interval - 1)%ub = 0.5_dp*(sr_env%subintervals(interval)%lb + sr_env%subintervals(interval)%ub)
404 subintervals(2*interval - 1)%conv = 0.5_dp*sr_env%subintervals(interval)%conv
405 subintervals(2*interval - 1)%fa = sr_env%subintervals(interval)%fa
406 subintervals(2*interval - 1)%fb = zdata_next(4*interval - 3)
407 subintervals(2*interval - 1)%fc = sr_env%subintervals(interval)%fb
408 subintervals(2*interval - 1)%fd = zdata_next(4*interval - 2)
409 subintervals(2*interval - 1)%fe = sr_env%subintervals(interval)%fc
411 subintervals(2*interval)%lb = subintervals(2*interval - 1)%ub
412 subintervals(2*interval)%ub = sr_env%subintervals(interval)%ub
413 subintervals(2*interval)%conv = subintervals(2*interval - 1)%conv
414 subintervals(2*interval)%fa = sr_env%subintervals(interval)%fc
415 subintervals(2*interval)%fb = zdata_next(4*interval - 1)
416 subintervals(2*interval)%fc = sr_env%subintervals(interval)%fd
417 subintervals(2*interval)%fd = zdata_next(4*interval)
418 subintervals(2*interval)%fe = sr_env%subintervals(interval)%fe
420 zdata_next(4*interval - 3:4*interval) = cfm_null
423 DO interval = nintervals + 1, nintervals_exist
424 subintervals(interval + nintervals) = sr_env%subintervals(interval)
426 DEALLOCATE (sr_env%subintervals)
430 ALLOCATE (sr_env%integral, sr_env%integral_conv, &
431 sr_env%integral_abc, sr_env%integral_cde, sr_env%integral_ace)
440 ALLOCATE (subintervals(nintervals))
442 rscale = 1.0_dp/real(nintervals, kind=
dp)
444 DO interval = 1, nintervals
446 subintervals(interval)%lb = sr_env%tnodes(4*interval - 3)
447 subintervals(interval)%ub = sr_env%tnodes(4*interval + 1)
448 subintervals(interval)%conv = rscale*sr_env%conv
450 subintervals(interval)%fa = zdata_next(4*interval - 3)
451 subintervals(interval)%fb = zdata_next(4*interval - 2)
452 subintervals(interval)%fc = zdata_next(4*interval - 1)
453 subintervals(interval)%fd = zdata_next(4*interval)
454 subintervals(interval)%fe = zdata_next(4*interval + 1)
460 zdata_next(1:npoints) = cfm_null
467 sr_env%error = sr_env%error_conv
468 nintervals_exist =
SIZE(subintervals)
470 DO interval = 1, nintervals_exist
471 rscale = subintervals(interval)%ub - subintervals(interval)%lb
472 CALL do_simpson_rule(sr_env%integral_ace, &
473 subintervals(interval)%fa, &
474 subintervals(interval)%fc, &
475 subintervals(interval)%fe, &
477 CALL do_simpson_rule(sr_env%integral_abc, &
478 subintervals(interval)%fa, &
479 subintervals(interval)%fb, &
480 subintervals(interval)%fc, &
482 CALL do_simpson_rule(sr_env%integral_cde, &
483 subintervals(interval)%fc, &
484 subintervals(interval)%fd, &
485 subintervals(interval)%fe, &
492 CALL do_boole_rule(sr_env%integral_abc, &
493 subintervals(interval)%fa, &
494 subintervals(interval)%fb, &
495 subintervals(interval)%fc, &
496 subintervals(interval)%fd, &
497 subintervals(interval)%fe, &
498 0.5_dp*rscale, sr_env%integral_cde)
504 error_rdata(:, :) = abs(error_zdata(:, :))
505 CALL cp_fm_trace(sr_env%error_fm, sr_env%weights, subintervals(interval)%error)
507 sr_env%error = sr_env%error + subintervals(interval)%error
510 IF (subintervals(interval)%error <= subintervals(interval)%conv)
THEN
512 sr_env%error_conv = sr_env%error_conv + subintervals(interval)%error
516 IF (sr_env%error <= sr_env%conv)
THEN
527 DO interval = nintervals_exist, 1, -1
528 interval_exists = .false.
529 my_bound = subintervals(interval)%lb
530 DO jpoint = 1, nintervals_exist
531 IF (subintervals(jpoint)%ub == my_bound)
THEN
532 interval_exists = .true.
536 IF (.NOT. interval_exists)
THEN
539 ELSE IF (interval_converged)
THEN
549 ALLOCATE (errors(nintervals_exist), inds(nintervals_exist))
552 DO interval = 1, nintervals_exist
553 errors(interval) = subintervals(interval)%error
555 IF (subintervals(interval)%error > subintervals(interval)%conv) &
556 nintervals = nintervals + 1
559 CALL sort(errors, nintervals_exist, inds)
561 IF (nintervals > 0) &
562 ALLOCATE (sr_env%subintervals(nintervals))
565 DO ipoint = nintervals_exist, 1, -1
566 interval = inds(ipoint)
568 IF (subintervals(interval)%error > subintervals(interval)%conv)
THEN
569 nintervals = nintervals + 1
571 sr_env%subintervals(nintervals) = subintervals(interval)
575 interval_exists = .false.
576 my_bound = subintervals(interval)%lb
577 DO jpoint = 1, nintervals_exist
578 IF (subintervals(jpoint)%ub == my_bound)
THEN
579 interval_exists = .true.
583 IF (.NOT. interval_exists)
THEN
586 ELSE IF (interval_converged)
THEN
595 interval_exists = .false.
596 interval_converged = .false.
597 my_bound = subintervals(interval)%ub
598 DO jpoint = 1, nintervals_exist
599 IF (subintervals(jpoint)%lb == my_bound)
THEN
600 interval_exists = .true.
601 IF (subintervals(jpoint)%error <= subintervals(jpoint)%conv) interval_converged = .true.
605 IF (.NOT. interval_exists .OR. interval_converged)
THEN
611 DEALLOCATE (errors, inds)
614 DEALLOCATE (subintervals)
616 CALL timestop(handle)