Fortran Coder

标题: 一个数组运算的问题 [打印本页]

作者: andy8496    时间: 2017-8-4 18:27
标题: 一个数组运算的问题
遇到一个数组片段的问题。

代码1:
do i=1,size(ia_index,2)
  ra_val(ia_index(1,i)) = ra_val(ia_index(1,i)) + ra_add(i)
  ra_val(ia_index(2,i)) = ra_val(ia_index(2,i)) + ra_add(i)
enddo

代码2:
  ra_val(ia_index(1,:)) = ra_val(ia_index(1,:)) + ra_add
  ra_val(ia_index(2,:)) = ra_val(ia_index(2,:)) + ra_add

我一直以为代码1和2得到的结果是一样的,今天遇到一个问题发现并不一样,代码1的结果是正确的,当ia_index(1,:)表示的某个下标索引出现1-2次的时候,两段代码结果一致,当出现3次的时候,代码1可以累加,代码2中ra_val索引指向的成员就只有最后一个值,没有累加。不知何故!

事实上,我最初使用代码2之前,由于不放心,是做过实验的:
[Fortran] 纯文本查看 复制代码
program main
implicit none
   integer :: ia(5),ib(9)
   ia = [1,5,5,5,2]
   ib = 0
   ib(ia) = ib(ia) + [1,1,1,1,1]
   write(*,*) ib
   stop
end program

结果见附件图片,无误,能够叠加我才这么写的,真是郁闷了...
还请大伙指点,先谢谢了!


aa.jpg (15.14 KB, 下载次数: 458)

aa.jpg

作者: fcode    时间: 2017-8-5 04:48
It's hard to check why ? May be we need more code !
作者: andy8496    时间: 2017-8-5 10:21
我编了一组数据:
[Fortran] 纯文本查看 复制代码

program main
implicit none
  real :: ra_val1(11),ra_val2(11),ra_add(9)
  integer :: ia_index(2,9)
  
  integer :: i
   
  ra_add = [0.014000,0.016500,0.016500,0.014000,0.014000,0.014750,0.014750,0.010750,0.010750]
  
  ia_index(:,1) = [ 1  , 2  ]
  ia_index(:,2) = [ 1  , 9  ]
  ia_index(:,3) = [ 1  , 11 ]
  ia_index(:,4) = [ 2  , 3  ]
  ia_index(:,5) = [ 3  , 4  ]
  ia_index(:,6) = [ 4  , 8  ]
  ia_index(:,7) = [ 4  , 10 ]
  ia_index(:,8) = [ 5  , 6  ]
  ia_index(:,9) = [ 6  , 7  ]
  

  
  ra_val1 = 0.
  ra_val2 = 0.
  
  ! 写法1
  do i=1,9
    ra_val1(ia_index(1,i)) = ra_val1(ia_index(1,i)) + ra_add(i)
    ra_val1(ia_index(2,i)) = ra_val1(ia_index(2,i)) + ra_add(i)
  enddo
  
  ! 写法2
  ra_val2(ia_index(1,:)) = ra_val2(ia_index(1,:)) + ra_add
  ra_val2(ia_index(2,:)) = ra_val2(ia_index(2,:)) + ra_add
  
  
  write(*,"(A6,2A12)") "No.","Val1","Val2"
  do i=1,11
    write(*,"(I6,2F12.6)") i,ra_val1(i),ra_val2(i)
  enddo
  
  stop
end program

aa.jpg (29.29 KB, 下载次数: 430)

aa.jpg

作者: andy8496    时间: 2017-8-5 10:23
本帖最后由 andy8496 于 2017-8-5 10:27 编辑
fcode 发表于 2017-8-5 04:48
It's hard to check why ? May be we need more code !

上面的代码为什么写法1和2结果不同?理解不了,我认为结果应该是一样的,都是写法1的结果。还请老大赐教!
作者: 维尼猴    时间: 2017-8-5 15:52
本帖最后由 维尼猴 于 2017-8-5 15:58 编辑

哇塞,竟然能这样操作,但感觉可能和缓存的处理有关系吧

我试了一下,比如加个 中间变量  temp

[Fortran] 纯文本查看 复制代码
temp = ia_index(1,:)
ra_val3(temp ) = ra_val3(temp ) + ra_add


这样处理之后是可以累加的,可能是你数组声明成  2,9 的问题,内存不够连续??
来我试试  9,2 的

诶………………我试了一下,失败了,把  ia_index 交换了维数也没有得到对应效果

作者: 维尼猴    时间: 2017-8-5 16:00
我补充测试了一下,似乎下标数组不能取切片

[Fortran] 纯文本查看 复制代码
program main
implicit none
  integer :: ia(5),ib(9)
  ia = [1,5,5,5,2]
  ib = 0
  ib(ia) = ib(    ia  (:)   ) + [1,1,1,1,1]
  write(*,*) ib
  stop
end program


只是在你原来例子上给ia加了个切片操作,就不能累加了

作者: andy8496    时间: 2017-8-5 21:29
唉,我都忘了我哪些地方用了这种写法,但是肯定不止目前发现的这个位置。真是郁闷了!
要是能通过什么编译选项能设置解决就好了!
谢了!
作者: 维尼猴    时间: 2017-8-6 16:07
andy8496 发表于 2017-8-5 21:29
唉,我都忘了我哪些地方用了这种写法,但是肯定不止目前发现的这个位置。真是郁闷了!
要是能通过什么编译 ...

我感觉这种有些“讨巧”的语法,还是少用心里更安一些呢
作者: pasuka    时间: 2017-8-6 17:22
ia_index的第1行数据为啥非要设置重复元素?

作者: andy8496    时间: 2017-8-6 18:30
pasuka 发表于 2017-8-6 17:22
ia_index的第1行数据为啥非要设置重复元素?

这个是实际应用的需要。
作者: pasuka    时间: 2017-8-6 21:28
andy8496 发表于 2017-8-6 18:30
这个是实际应用的需要。

症结即在于此呀,没有看出来?
作者: andy8496    时间: 2017-8-6 21:40
本帖最后由 andy8496 于 2017-8-6 21:41 编辑
pasuka 发表于 2017-8-6 21:28
症结即在于此呀,没有看出来?

真看不出来。我试验用的那段代码,咋就能按照我的理解运行呢?索引数组也是有重复元素的:
program main
implicit none
   integer :: ia(5),ib(9)
   ia = [1,5,5,5,2]
   ib = 0
   ib(ia) = ib(ia) + [1,1,1,1,1]
   write(*,*) ib
   stop
end program

aaa.jpg (15.14 KB, 下载次数: 237)

aaa.jpg

作者: chiangtp    时间: 2017-8-6 23:20
a many-one array section must NOT Assign To NOR Read In
Compiler may not Check at Runtime owing to the Performance Cost

many-one-array.pdf

74.9 KB, 下载次数: 13


作者: andy8496    时间: 2017-8-7 09:15
好吧,总算看到一点依据了。虽然心中不服,比如:五楼的大哥的方法为什么又可行?
还有两点疑问,还要继续麻烦大家指点:
chiangtp附件中提到的,
①a many-one array section must NOT Assign To NOR Read In Compiler may not Check at Runtime owing to the Performance Cost
VS+IVF中,Debug中如何打开这个检查?
②the use of vector subscripting is very inefficient and should not be used unless absolutely necessary (although easy to coding)
这句话该如何理解?数据片段只是写起来方便,但是运行效率低?
作者: pasuka    时间: 2017-8-7 11:04
本帖最后由 pasuka 于 2017-8-7 11:17 编辑
andy8496 发表于 2017-8-6 21:40
真看不出来。我试验用的那段代码,咋就能按照我的理解运行呢?索引数组也是有重复元素的:
program main
  ...

左端赋值包含重复元素的切片这样的写法本身就是不正确的,起码IBM XL和Intel Fortran的用户手册都明确强调这点
传送门1: http://w3.pppl.gov/~hammett/comp/f90tut/f90.tut5.html
Note that a vector subscript with duplicate values cannot appear on the left-hand side of an assignment as it would be ambiguous. Thus,
       b( (/ 1, 7, 3, 7 /) ) = (/ 1, 2, 3, 4 /)
is illegal.
传送门2: https://www.ibm.com/support/know ... ectorsubscript.html
An array section with a vector subscript in which two or more elements of the vector subscript have the same value is called a many-one section. Such a section must not:
Appear on the left side of the equal sign in an assignment statement传送门3: https://software.intel.com/en-us/node/678553
An array section with a vector subscript that has two or more elements with the same value is called a many-one array section. For example:
  REAL A(3, 3), B(4)  INTEGER K(4)  ! Vector K has repeated values  K = (/3, 1, 1, 2/)  ! Sets all elements of A to 5.0  A = 5.0  B = A(3, K)
The array section A(3,K) consists of the elements:
  A(3, 3)  A(3, 1)  A(3, 1)  A(3, 2)
A many-one section must not appear on the left of the equal sign in an assignment statement, or as an input item in a READ statement.



作者: chiangtp    时间: 2017-8-7 13:03
1. 15F, posuka見多識廣, 佩服
2. whole array的assignment(多個要做的規則), 定義上是in-any-order+parallel(如果能), runtime時 "五楼的大哥的方法", "Debug打开检查" 都是compiler/compile-time-options/... dependent
3. "vector subscripting is very inefficient "應可理解(?), 我沒有理想的實證

作者: chiangtp    时间: 2017-8-7 13:05
抱歉, pasuka誤為posuka
作者: pasuka    时间: 2017-8-7 13:52
本帖最后由 pasuka 于 2017-8-7 14:01 编辑

再补充一些内容
Fortran 2003 handbook的第6章,第184页如是说:
A many-one section must not appear on the left of the equal sign in an assignment statement or as an input item in a READ statement. The reason is that the result will depend on the order of evaluation of the subscripts, which is not specified by the language. The results would not be predictable and the program containing such a statement would not be portable.

也即这个属于Fortran规范的盲点,lz有心的话,不妨邮件咨询RAL的J K Reid爵士或者I Duff
btw,无需见多识广,反倒是科学上网很重要

作者: andy8496    时间: 2017-8-7 14:52
多谢各位了!
现在补救的办法就是依赖编译器检查了,不知道可行否?在哪里设置?


2. whole array的assignment(多個要做的規則), 定義上是in-any-order+parallel(如果能), runtime時 "五楼的大哥的方法", "Debug打开检查" 都是compiler/compile-time-options/... dependent

这个设置我没找到在哪儿……

setting.jpg (75.14 KB, 下载次数: 220)

setting.jpg

作者: pasuka    时间: 2017-8-7 15:34
本帖最后由 pasuka 于 2017-8-7 15:37 编辑
andy8496 发表于 2017-8-7 14:52
多谢各位了!
现在补救的办法就是依赖编译器检查了,不知道可行否?在哪里设置?

直接改回方法1不行吗?
在ia index数组包含重复元素且不能剔除的前提下,IBM XL和Intel Fortran编译器手册均明确对方法2亮红牌,这下总得遵守规则被罚下场吧?
仍然不服气的话,建议研读g95或gfortran前端代码,自行实现需要的功能。g95的代码值得看看,既能加深Fortran的了解,也能助益C语言进阶
作者: andy8496    时间: 2017-8-7 15:55
惭愧!惭愧!
已经是心服口服了!
目前遇到的困难是:工程已经相当大,目前不知道那些位置有这个问题;现在想改……
作者: chiangtp    时间: 2017-8-7 16:16
你ˋ可能還需要考慮另一個問題:

array-expression-assignment.pdf

117.25 KB, 下载次数: 3


作者: 维尼猴    时间: 2017-8-7 19:57
chiangtp 发表于 2017-8-7 16:16
你ˋ可能還需要考慮另一個問題:

哇,学到好多~~



能问一下这个资料从哪里可以找到么?我可能找法不对,没有搜到


作者: chiangtp    时间: 2017-8-8 00:02
维尼猴 发表于 2017-8-7 19:57
哇,学到好多~~

http://pan.baidu.com/s/1i5l8Ypf
作者: 维尼猴    时间: 2017-8-8 16:50
chiangtp 发表于 2017-8-8 00:02
http://pan.baidu.com/s/1i5l8Ypf

O(∩_∩)O谢谢!!
作者: mltx    时间: 2017-8-9 21:51
[Fortran] 纯文本查看 复制代码
program main
implicit none
   integer :: ia(5),ib(9)
   ia = [1,5,5,5,2]
   ib = 0
   ib(ia) = ib(ia) + [1,1,1,1,1]
   write(*,*) ib
   stop
end program

没特别仔细看,试着说两句:
我用CVF和ELF90试了一下,并不累加(我还是用老古董的编译器)。
感觉写法2的语法就不应该累加。因为从语法上看,看不出来右边的ib一定要用更新的值,或对ib(5)的三次操作要累加起来。
而写法1就非常明确,等号右边总是用更新的值,就有累加效果。

写法2的累加效果可能是编译器dependent,还是用写法1比较稳妥。
作者: chiangtp    时间: 2017-8-10 13:18
1. [代碼2] assign to many-one array section, 語法的 logical (非compiler責任) error, 語意就遑論了
2. [Squential DO] 能否coding為 [ARRAY assignment], 端看要做得多個能否 in any order
3 [代碼1 loop], DO I=1,5 與 DO I=5,1,-1 是不同的, 是否為andy8496的本意???
作者: mltx    时间: 2017-8-10 19:46
找到了依据(by Michael Metcalf / CERN CN-AS):

Note that a vector subscript with duplicate values cannot appear on the left-hand side of an assignment as it would be ambiguous. Thus,

       b( (/ 1, 7, 3, 7 /) ) = (/ 1, 2, 3, 4 /)

is illegal.

作者: andy8496    时间: 2017-8-10 21:25
本帖最后由 andy8496 于 2017-8-10 21:27 编辑

当时就是想以最少的代码实现一个累加,也略微表示怀疑,但是试验后来发现成功了,就采用了,但后来才发现了这样一个大bug。
等号左边的数组片段中,表示下表的数组不能有重复元素。后来在各位的帮助下找到依据了,如15楼pasuka的资料就非常详细。

当时做这样的语法规定,个人理解是为了防止在数组初始化的时候,某个重复下表指向的成员得到了不确定的值,因为数组片段可能采用了并行的方法来赋值(纯猜测),这样就不知道最终这些重复下标的成员最终得到的是何值。
个人认为:
①在单纯的数组赋值时,编译器对于这种错误程序运行时应该报错,而不应该任由错误发生。如:
b( (/ 1, 7, 3, 7 /) ) = (/ 1, 2, 3, 4 /) ! 运行时应该有报错提示
②对于累加的需求,应该从语法层面予以支持,如:
b( (/ 1, 7, 3, 7 /) ) =b( (/ 1, 7, 3, 7 /) ) +1 !等号左右有相同的数组片段应该实现循环相同的效果。

最后,再次感谢各位的热心帮助!





欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2