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编译器。现在程序就耗在这，不止到如何解决。。。如果大家发现了什么问题，能提供相应的思路的话，那就太感谢了。

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编译的，而前面是现成的）。

最后，谢谢大家的关心！