Fortran Coder

标题: Lapack中几个难以理解的现象 [打印本页]

作者: cstg 时间: 2017-8-14 11:39
标题: Lapack中几个难以理解的现象
[attach]1557[/attach][attach]1557[/attach]最近在编写一个算法时发现了以下两个难以理解的问题：
第一个是使用了write语句输出后会影响后续计算结果。。：

[Fortran] 纯文本查看 复制代码

        write(*,*) resid

          

            call zgeqp3(n, d, resid, n, workint,

     &               workqr,workqr(1+min(n,d)), 

     &      worklen, rwork, info)



           write(*,*) resid

这里zgeqp3是一个计算QR分解的subroutine，是来自于lapack库的。residue作为存放要分解的矩阵，在经过zgeqp3后上三角部分存放得到的R矩阵。
可以看到里面有两处write语句。非常奇怪的是，如果把第一个write语句注释掉，则计算得到的结果residue会出现NAN；反之，如果不注释掉，则计算得到的residue显示正常。这实在难以理解。

第二个是两个数相乘得到了NAN。

[Fortran] 纯文本查看 复制代码

         do 15 i=1, d, 1

             call zcopy(min(d+1-i,nrv),resid(1+(d-i)*n),1,

     &      u(1,(-1+i)*(ncv+1)+1),1) 

             temp1=(1.0d0)/rnorm

         do 16 jj=1, min(d+1-i,nrv)

             write(*,*) u(jj,(i-1)*(ncv+1)+1)

             write(*,*) cmplx(temp1,0.0d0)

             write(*,*) u(jj,

     &      (i-1)*(ncv+1)+1)*cmplx(temp1,0.0d0)  

              u(jj,(i-1)*(ncv+1)+1)=u(jj,

     &     (i-1)*(ncv+1)+1)*cmplx(temp1,0.0d0)     

   16    continue        

   15  continue

这里是将一个矩阵residue复制到另一个矩阵u中，之后再对u中元素乘以一个数。cmplx是使得temp1从实数转为复数。三个write语句输出了每次乘法的过程。结果见图：

可以看到我相乘的两个数中并没有nan，且也不存在0与无穷相乘的情况。但结果却有nan。这对我完全无法理解。。。

我使用的是MinGW编译器。现在程序就耗在这，不止到如何解决。。。如果大家发现了什么问题，能提供相应的思路的话，那就太感谢了。

1.png (10.87 KB, 下载次数: 282)

作者: pasuka 时间: 2017-8-14 12:37
无法提供可以再现问题的最短fortran代码，就只能对照zgeqp3的函数说明看看有没有参数输入错误
http://www.netlib.org/lapack/exp ... b3f0d9468d99c28d36e

作者: chiangtp 时间: 2017-8-14 13:12
建議試試:
1. compile-time-opions儘量non-optimization
2. 用不同的compilers試試
"lapack库"是用別人的".lib", 還是自己建立的

作者: cstg 时间: 2017-8-14 13:53

pasuka 发表于 2017-8-14 12:37
无法提供可以再现问题的最短fortran代码，就只能对照zgeqp3的函数说明看看有没有参数输入错误
http://www.n ...

谢谢，我检查过了，应该是对的。而且我用其他例子简单测试过这个函数的用法，并没有出现相同问题。

作者: cstg 时间: 2017-8-14 13:54

pasuka 发表于 2017-8-14 12:37
无法提供可以再现问题的最短fortran代码，就只能对照zgeqp3的函数说明看看有没有参数输入错误
http://www.n ...

而且无论是否出现nan，zgeqp3返回的info都是0，就是说正常结束。。。

作者: cstg 时间: 2017-8-14 14:31

chiangtp 发表于 2017-8-14 13:12
建議試試:
1. compile-time-opions儘量non-optimization
2. 用不同的compilers試試

是用别人的lib，具体的说是用如下网址中的“Prebuilt dynamic libraries using Mingw”
：http://icl.cs.utk.edu/lapack-for-windows/lapack/#libraries

我也试着用其他的编译器，比如VS的编译器，但是使用时却提示我找不到库函数（无法解析的命令）。。。而我似乎是设置正确的：

6.png (31.11 KB, 下载次数: 268)

5.png (66.68 KB, 下载次数: 278)

4.png (63.23 KB, 下载次数: 279)

3.png (59.55 KB, 下载次数: 271)

2.png (32.47 KB, 下载次数: 266)

作者: cstg 时间: 2017-8-14 17:25
额！我大概找到原因了，很可能是库文件的问题。。。我刚才换了一个版本的库文件后，这些神奇的现象就都消失了。。具体原因不知道。我猜想的是原来的库文件是专门for windows用户的，而我在编译时是用MinGw的shell，会不会这里出现了问题（这方面我不懂，很可能是瞎想）；后来我换的新版本就不是给windows用户的了（且这次库文件是用MinGw编译的，而前面是现成的）。

最后，谢谢大家的关心！

作者: chiangtp 时间: 2017-8-14 19:42
恭喜, 解決了問題, 一樂也。

三個可以算是 "难以理解"的, 與君同樂

(a) 光看主程式, "WRITE(*,*) 5", write出來說不是"5", 不可能吧!

[Fortran] 纯文本查看 复制代码

PROGRAM test

  IMPLICIT NONE

  WRITE(*,*) 5

  CALL abc(5)

  WRITE(*,*) 5

  CALL abc(5) 

  WRITE(*,*) 5

CONTAINS

  SUBROUTINE abc(five)

    IMPLICIT NONE

    INTEGER :: five

    five = five + 1

  END SUBROUTINE abc

END PROGRAM test

compile-time給wraning或fatal error, 或是run-time出現不知所云error的, 是"嚴謹"的compilers
run-time告訴你"5"不是"5", 不是Fortran的錯, 也不是compiler的錯(compiler太老實了點), 是程式員的錯(logical error)
compile-time沒任何訊息, run-time告訴你"5"就是"5"的compilers, 太自作聰明了
用Fortrn 90"INTENT"attribute, 是一個好習慣, 但是, 請看 (b)

(b) 光看副程式, 兩個"WRITE(*,*) x", 沒理由值不同吧?

[Fortran] 纯文本查看 复制代码

PROGRAM test

  IMPLICIT NONE

  REAL :: x

  x = 5.0

  CALL abc(x, x) !---> OK, if: CALL ABC((x), x)

CONTAINS

  SUBROUTINE abc(x, y)

    IMPLICIT NONE

    REAL, INTENT(IN   ) :: x

    REAL, INTENT(  OUT) :: y

    WRITE(*,*) 'x=', x

    y = x + 1.0

    WRITE(*,*) 'x=', x

  END SUBROUTINE abc

END PROGRAM test

原因出在主程式 "CALL abc(x, x)"
主/副程式, "語法"都是正確的, 程式員的錯(logical error)
不要以為這是"無聊"的把戲, "x"帶訊息進去責任已了, 能省則省新的訊息由"x"帶出來, 挺合理的
多加一對小括號 "CALL abc((x), x)", 就完全沒問題了

(c) 說一套(syntax defined by language standrad), 做一套(compilation by compiler)

[Fortran] 纯文本查看 复制代码

PROGRAM test

  IMPLICIT NONE

  REAL :: a, b, c

  a=1.0; b=2.0; c=3.0

  CALL abc(a, b, c)

  WRITE(*,*) c

CONTAINS

  SUBROUTINE abc(a, b, c)

    IMPLICIT NONE

    REAL, INTENT(IN   ) :: a

    REAL, INTENT(INOUT) :: b

    REAL, INTENT(  OUT) :: c

    b = a + b

   !WRITE(*,*) c

  END SUBROUTINE abc

END PROGRAM test



! All INTENT(OUT) Dummy Arguments become Undefined on Entry to a Procedure.

! (the intension is that it be used only to pass information out, so it

! becomes undefined on entry to the procedure)

! You cannot use an INTENT(OUT) dummy until it has been defined

(a), (b), (c) compile-time/run-time 會怎樣, compile-time "options" dependent, "compiler" dependent,
甚至會是 compiler "version" dependent, 怎辦？ ["SAVE"也是一樣]

不怎辦, 都是程式員的錯 (logical errors)

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)