Let's break down what's going wrong with your code and fix it step by step—your trapezoidal integration logic is actually correct, so the issue isn't with the method itself, nor is it related to parallel computing (your code doesn't include any parallel directives anyway). The problem lies in several critical bugs in variable handling and initialization.

Key Bugs in Your Original Code #

• Integer/Real Type Mismatch: You declared j as integer(8) but assigned it j=0. (a real value). This causes compiler warnings at best, and worse, leads to invalid array indexing when you set sum(j)=0.—since sum is dimensioned (1317), index 0 is out of bounds, triggering undefined behavior.
• Uninitialized Arrays: sum and factor are never properly initialized. You only set the invalid index 0 to 0, leaving all other elements with garbage values that get included in your integral calculation.
• Complex Type Mismatch: You use cmplx(..., kind=16) to create quadruple-precision complex numbers, but your array R is declared as complex(8) (double-precision). This type mismatch leads to unintended data conversion or truncation.
• Potential Naming Conflict: sum is a built-in Fortran function—using it as a variable name can cause confusion and unexpected behavior.

Corrected Fortran Code #

Here's a fixed version with comments explaining the changes:

program array_integration
    implicit none
    integer(8) :: i, j
    real(8) :: dt
    ! Rename 'sum' to avoid conflict with Fortran's built-in sum() function
    real(8) :: RR(1317,26123), t(26123), IR(1317,26123), f(8,26123)
    complex(8) :: integral_sum(1317), factor(1317), fnl(1317)
    complex(8), dimension(:,:), allocatable :: R
    
    ! Initialize all complex arrays to zero to avoid garbage values
    integral_sum = (0.0d0, 0.0d0)
    factor = (0.0d0, 0.0d0)
    
    ! Allocate memory for the complex array
    allocate(R(1317,26123))
    
    ! Use meaningful unit numbers (avoid 0, which is often tied to stdin)
    Open(unit=10, File='/Users/clyrionchen/Desktop/result1-4.log', status='old')
    Open(unit=11, File='/Users/clyrionchen/Desktop/ReR', status='old')
    Open(unit=22, File='/Users/clyrionchen/Desktop/ImR', status='old')
    
    ! Read input data (ensure file formats match array dimensions)
    do i = 1, 26123
        read(11, *) RR(:, i)
        read(10, *) f(:, i)
        read(22, *) IR(:, i)
    end do
    
    ! Extract time points from the f array
    t(:) = f(1, :)
    
    ! Create complex array with matching kind (8 for double-precision)
    R(:,:) = cmplx(RR(:,:), IR(:,:), kind=8)
    
    ! Trapezoidal integration logic (your original logic was correct!)
    do j = 1, 1317
        do i = 2, 26123
            ! Calculate area of current trapezoid
            factor(j) = (R(j,i) + R(j,i-1)) / 2.0d0 * (t(i) - t(i-1))
            ! Accumulate the integral
            integral_sum(j) = integral_sum(j) + factor(j)
        end do
    end do
    
    ! Print formatted results for readability
    do j = 1, 1317
        write(*, '(A, I5, A, 2F16.8)') 'Integral for group ', j, ': ', integral_sum(j)
    end do
    
    ! Clean up: close files and deallocate memory
    close(10)
    close(11)
    close(22)
    deallocate(R)
    
end program array_integration

Additional Notes #

• If you later want to add parallelization (to speed up the 1317 groups), you can use OpenMP on the outer j loop. Add this line right before the do j=1,1317 loop:

  !$OMP PARALLEL DO PRIVATE(i,factor) REDUCTION(+:integral_sum)
  

  Compile with OpenMP support (e.g., gfortran -fopenmp your_code.f90).
• Always verify your input files: make sure ReR and ImR contain exactly 26123 lines, each with 1317 values, to match your array dimensions.

内容的提问来源于stack exchange，提问作者ClyrionChen