Lapack中几个难以理解的现象

cstg

最近在编写一个算法时发现了以下两个难以理解的问题：
第一个是使用了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编译器。现在程序就耗在这，不止到如何解决。。。如果大家发现了什么问题，能提供相应的思路的话，那就太感谢了。

pasuka

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

chiangtp

建議試試:
1. compile-time-opions儘量non-optimization
2. 用不同的compilers試試
"lapack库"是用別人的".lib", 還是自己建立的

cstg

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

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

cstg

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

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

cstg

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的编译器，但是使用时却提示我找不到库函数（无法解析的命令）。。。而我似乎是设置正确的：

cstg

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

最后，谢谢大家的关心！

chiangtp

恭喜, 解決了問題, 一樂也。

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

(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)