Запуск циклов только для указанных наборов целых чисел в Фортран - PullRequest
Новая
       25

Запуск циклов только для указанных наборов целых чисел в Фортран

0 голосов
/ 13 февраля 2020 
        do i=1,n
                s=0
                do l=1,n
                do m=1,n


                    s=s-a(i,l,m)*q0(l)*q0(m)

                end do
                end do


                f0(i)=s-g(i)*q0(i)
        end do

Это часть моего кода. Поскольку мне приходится одновременно запускать три цикла друг в друге, общее выполнение становится ужасно медленным.

Важным фактом является то, что здесь массив a (i, l, m) отличен от нуля только для набора значения a (l, m, n). Ниже приведен код для установки (i, l, m). 

do i=1,n
    do l=1,n
        do m=1,n

        if(i.eq.l+m .or. i.eq.-l+m .or. i.eq.l-m) then
        a1=1
        else 
        a1=0
        end if



        if(i+l+m.eq.n+n+2 .or. i-l+m.eq.n+n+2 .or. i+l-m.eq.n+n+2 .or. i-l-m.eq.n+n+2) then
        b1=1
        else
        b1=0
        end if

    a(i,l,m)=(a1-b1)    (!multiplied with some long function, erased for ease of understanding)

end do
    end do
        end do

Теперь, есть ли какой-нибудь способ в фортране запустить l oop только для значений (i, l, m) для которого a (i, l, m) отличен от нуля? (ненулевое значение устанавливает только значение только в вычислениях, как можно видеть). Это сэкономит огромное количество времени.

Ответы [ 2 ]

3 голосов
/ 13 февраля 2020

Вам не нужно хранить весь массив a, вы можете вычислять его элементы по мере необходимости, а также использовать защитные элементы (дополнительные элементы массива, чтобы избежать индексов за пределами границ), чтобы избежать условий if. Вот как вы можете это сделать, уменьшив проблему до O (n ^ 2) и используя намного меньше памяти. Также обратите внимание, что я предоставил полную тестовую программу, это намного облегчает ответы на запросы, НАМНОГО проще - пожалуйста, сделайте это сами в будущем! 

ijb@ianbushdesktop ~/work/stack $ cat o3.f90
Program o3

  Implicit None

  Integer, Parameter :: wp = Selected_real_kind( 12, 70 )

  Real( wp ), Dimension( :, :, : ), Allocatable :: a

  Real( wp ), Dimension( : ), Allocatable :: q0, g, f0, s2

  Real( wp ) :: a1, b1
  Real( wp ) :: s

  Integer :: n
  Integer :: start, finish, rate
  Integer :: i, l, m

  Write( *, * ) 'n ?'
  Read ( *, * ) n 

  Allocate( a( 1:n, 1:n, 1:n ) )
  Allocate( q0( 1:n ) )
  Allocate( f0( 1:n ) )
  Allocate(  g( 1:n ) )
  Allocate( s2( -4 * n - 2:4 * n + 2 ) ) ! guards to avoid out of bounds - haven't thought very carefully about
                                         ! what they should be!!

  Call Random_number( q0 )
  Call Random_number( g )

  Call system_clock( start , rate )
  b1 = 0.0_wp
  do i=1,n
     do l=1,n
        do m=1,n

           if(i.eq.l+m .or. i.eq.-l+m .or. i.eq.l-m) then
              a1=1.0_wp
           else 
              a1=0.0_wp
           end if

           if(i+l+m.eq.n+n+2 .or. i-l+m.eq.n+n+2 .or. i+l-m.eq.n+n+2 .or. i-l-m.eq.n+n+2) then
              b1=1.0_wp
           else
              b1=0.0_wp
           end if

           a(i,l,m)=(a1-b1)    !(multiplied with some long function, erased for ease of understanding)

        end do
     end do
  end do
  do i=1,n
     s=0.0_wp
     do l=1,n
        do m=1,n

           s=s-a(i,l,m)*q0(l)*q0(m)

        end do
     end do
     f0(i)=s-g(i)*q0(i)
  end do
  Call system_clock( finish, rate )
  Write( *, * ) 'Sum f0, time: ', Sum( f0 ), Real( finish - start ) / rate

  Call system_clock( start , rate )
  s2 = 0.0_wp
  Do l = 1, n
     Do m = 1, n

        ! First condition
        i = l + m
        a1 = 1.0_wp
        s2( i ) = s2( i ) - a1 * q0( l ) * q0( m )

        ! Second condition
        i = - l + m
        a1 = 1.0_wp
        s2( i ) = s2( i ) - a1 * q0( l ) * q0( m )

        ! Third condition
        i = l - m
        a1 = 1.0_wp
        s2( i ) = s2( i ) - a1 * q0( l ) * q0( m )

        ! Fourth Condition
        i = 2 * n + 2 - l - m
        b1 = 1.0_wp
        s2( i ) = s2( i ) - ( - b1 ) * q0( l ) * q0( m )

        ! Fifth Condition
        i = 2 * n + 2 + l - m
        b1 = 1.0_wp
        s2( i ) = s2( i ) - ( - b1 ) * q0( l ) * q0( m )


        ! Sixth Condition
        i = 2 * n + 2 - l + m
        b1 = 1.0_wp
        s2( i ) = s2( i ) - ( - b1 ) * q0( l ) * q0( m )

        ! Seventh Condition
        i = 2 * n + 2 + l + m
        b1 = 1.0_wp
        s2( i ) = s2( i ) - ( - b1 ) * q0( l ) * q0( m )


     End Do
  End Do
  Do i = 1, n
     f0( i ) = s2( i ) - g( i ) * q0( i )
  End Do
  Call system_clock( finish, rate )
  Write( *, * ) 'Sum f0, time: ', Sum( f0 ), Real( finish - start ) / rate
End Program o3
ijb@ianbushdesktop ~/work/stack $ gfortran --version
GNU Fortran (GCC) 7.4.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ijb@ianbushdesktop ~/work/stack $ gfortran -Wall -Wextra -std=f2008 -fcheck=all -O o3.f90
ijb@ianbushdesktop ~/work/stack $ ./a.out
 n ?
300
 Sum f0, time:   -23660.711846511185       0.446999997    
 Sum f0, time:   -23660.711846511185        1.00000005E-03
ijb@ianbushdesktop ~/work/stack $ gfortran -Wall -Wextra -std=f2008 -O3 o3.f90
ijb@ianbushdesktop ~/work/stack $ ./a.out
 n ?
300
 Sum f0, time:   -21932.467299817898       0.298999995    
 Sum f0, time:   -21932.467299817898        0.00000000    
ijb@ianbushdesktop ~/work/stack $ ./a.out
 n ?
1000
 Sum f0, time:   -238036.00437753636        52.4760017    
 Sum f0, time:   -238036.00437753636        2.00000009E-03
ijb@ianbushdesktop ~/work/stack $
2 голосов
/ 13 февраля 2020

Фортран, а также множество других языков программирования, имеют операторы jump, которые позволяют манипулировать итерацией l oop сверх стандартного элемента управления l oop. В Фортране это операторы CYCLE и EXIT:

оператор CYCLE: Выполнение итераций al oop может быть сокращено путем выполнения оператора CYCLE, принадлежащего конструкции Оператор EXIT: Оператор EXIT предоставляет один из способов завершения al oop или завершения выполнения другой конструкции.

Используя эти конструкции, теперь можно быстро циклически проходить циклы когда конкретный индекс не имеет отношения к вычислению. В случае OP можно сделать что-то вроде: 

do i=1,n
   s=0
   do l=1,n
      do m=1,n
         if (a(i,l,m) == 0) cycle
         s=s-a(i,l,m)*q0(l)*q0(m)
      end do
   end do
   f0(i)=s-g(i)*q0(i)
end do

Конечно, вы всегда должны учитывать, что это останется проблемой O (n ^ 3).

Там Тем не менее, больше информации о том, как вы строите свой 3d-массив a. Поскольку a(i,l,m) = a1 - b1 и a1 и b1 могут иметь значения 0 или 1 только в зависимости от условия, элемент a(i,l,m) отличается от 0, если выполняется только 1 из условий. Теперь очень легко проверить, что если первое условие выполнено: 

i == l+m .or. i == -l+m .or. i == l-m

, второе условие никогда выполнено: 

i+l+m == 2*n+2 .or. i-l+m == 2*n+2 .or. i+l-m == 2*n+2 .or. i-l-m == 2*n+2

Таким образом, только одно из условия могут быть выполнены одновременно. Это дает вам дополнительные рычаги для ускорения процесса и удаления внутреннего l oop, делая это O (n ^ 2): 

do i=1,n
  s=0
  do l=1,n
     m=i-l
     if (m > 0 .and. m <= n) s=s-a(i,l,m)*q0(l)*q0(m)
     m=i+l
     if (m > 0 .and. m <= n) s=s-a(i,l,m)*q0(l)*q0(m)
     m=l-i
     if (m > 0 .and. m <= n) s=s-a(i,l,m)*q0(l)*q0(m)
     m=2*n+2-i-l
     if (m > 0 .and. m <= n) s=s-a(i,l,m)*q0(l)*q0(m)
     m=2*n+2-i+l
     if (m > 0 .and. m <= n) s=s-a(i,l,m)*q0(l)*q0(m)
     m=-(2*n+2-i-l)
     if (m > 0 .and. m <= n) s=s-a(i,l,m)*q0(l)*q0(m)
     m=-(2*n+2-i+l)
     if (m > 0 .and. m <= n) s=s-a(i,l,m)*q0(l)*q0(m)
  end do
  f0(i)=s-g(i)*q0(i)
end do

Дальнейшие улучшения наверняка еще возможны.

...