作业2-1：矩阵协方差

备注：待亲手实验

讲评：

（1）编译器参数：

　　改makefile

　　1.降精度 -O3：精度误差累加可能会出错

　　2.xHost：intel编译优化，告诉编译器生成编译主机处理器上可用的最高指令集的指令。无脑加就行。

　　3.开向量模块（intel编译器）

修改过的makefile

FC=ifort
LD= ifort
FFLAGS= -qopenmp -mcmodel=large -O3 -xHost //增加 -O3 -xHost

main.b: main.o Covariance1.o Covariance.o check.o
    ${FC} ${FFLAGS} main.o Covariance1.o Covariance.o check.o -o main.b
main.o: main.f90
    ${FC} ${FFLAGS} -c main.f90
Covariance1.o: Covariance1.f90 
    ${FC} ${FFLAGS} -c Covariance1.f90
Covariance.o: Covariance.f90
    ${FC} ${FFLAGS} -c Covariance.f90
check.o: check.f90
    ${FC} ${FFLAGS} -c check.f90

clean:
    rm -rf *.o *.i main

（2）改代码

　　1.函数inline可以跨文件，直接拷过来

　　2.fortran：列优先，若代码为行逻辑则考虑先矩阵转置（可并行可分块，按cache层）；

　　　　　　　fortran语言：先声明，后申请空间（数组），然后才是语句

修改后converiance代码SUBROUTINE FindCOV(NVAR, NROW, NV, M1Type, COV)

        IMPLICIT NONE

        INTEGER, INTENT(IN)    :: NVAR              !Number of variables
        INTEGER, INTENT(IN)    :: NROW              !Number of rows in MType
        INTEGER, INTENT(IN)    :: NV                !Number of cals in MType
        REAL*8,    INTENT(IN)    :: M1Type(NVAR, NROW) !Data used to calculate covariance
        REAL*8,    INTENT(INOUT) :: COV(NVAR, NV)   !Covariance matrix, output

        INTEGER :: IVar1
        INTEGER :: IVar2
        INTEGER :: I2 
        REAL*8, ALLOCATABLE::MType(:, :)//数组声明
        REAL*8,EXTERNAL :: Covariance
        ALLOCATE(MType(NROW, NVAR))//数组申空间
        !$OMP PARALLEL DO private(IVar1,IVar2) shared(M1Type, MType) schedule(guided) //并行，每个并行块都要写
        DO IVar1 = 1, NVAR//数组转置
        DO IVar2 = 1, NROW
                MType(IVar1, IVar2) = M1Type(IVar2, IVar1) 
        ENDDO
        ENDDO
        !$OMP END PARALLEL DO
        !$OMP PARALLEL DO private(IVar1,IVar2,I2) shared(NVAR,COV,MType) schedule(guided) //并行 加调度参数
        DO IVar1 = 1, NVAR
        DO IVar2 = 1, NV
           I2 = mod(IVar1 + IVar2 - 1, NVAR) + 1
           COV(IVar1, IVar2) = Covariance(NROW, MType(:,IVar1), MType(:,I2))//正确遍历：列优先 
//fortran冒号：隐含遍历
　　     ENDDO
        ENDDO
        !$OMP END PARALLEL DO
        RETURN

    END SUBROUTINE FindCOV

    !DEC$ ATTRIBUTES FORCEINLINE :: Covariance
    REAL*8 FUNCTION Covariance(NObs, X1, X2)

        IMPLICIT NONE
        
        INTEGER, INTENT(IN) :: NObs
        REAL*8,    INTENT(IN) :: X1(NObs)
        REAL*8,    INTENT(IN) :: X2(NObs)

        INTEGER :: IObs, NumMissing
        REAL*8 :: MeanX1, MeanX2
        LOGICAL :: IsMissingPheno//将外面的函数拷进来，删掉external，一开始把external的函数忽略了。。
        NumMissing=0
        MeanX1=0.0
        MeanX2=0.0
        Covariance=0.0 
        !!$OMP PARALLEL DO private(Iobs) shared(Covariance, X1, X2, NObs) schedule(guided) 
        DO IObs=1, NObs
            IF(ISNAN(X1(IObs)) .OR. & //IsMissngPheno的展开
                ISNAN(X2(IObs)) .OR. &
                ABS(X1(IObs)-65.0) .LE. 1.e-03 .OR. ABS(X2(IObs)-65.0) .LE. 1.e-03) THEN
                NumMissing=NumMissing+1
       
            ELSE
                MeanX1=MeanX1+X1(IObs)
                MeanX2=MeanX2+X2(IObs)

                Covariance=Covariance+X1(IObs)*X2(IObs)
     
            ENDIF
            
        ENDDO
        !!$OMP END PARALLEL DO
        MeanX1=MeanX1/(NObs-NumMissing)
        MeanX2=MeanX2/(NObs-NumMissing)

        IF(NObs-NumMissing .LE. 0) THEN
            Covariance=0.0
        ELSE
            Covariance=Covariance/(NObs-NumMissing)-MeanX1*MeanX2
   
        ENDIF


        RETURN

    END FUNCTION Covariance

经过：

（1）将本地文件传到服务器：使用pscp显示unreachable。下载xshell 使用rz传：

https://blog.csdn.net/weixin_37909391/article/details/80530575

下载完xshell觉得很好看，我要舍弃putty了。

（2）解压zip：unzip+名字

https://www.cnblogs.com/zdz8207/p/3765604.html

（3）运行，需要编译。f90是fortran，需要ifort编译。

　　按ifort：服务器上有

　　使用命令： module load intel/2017.6.064 （module：环境模块 http://blog.sciencenet.cn/blog-468005-920774.html）

　　如何编译：http://blog.sciencenet.cn/blog-1982759-890635.html ：

1.直接生成可执行文件

假如我们手头上有一个intel的fortran编译器进行编译，就要采用如下命令：

ifort global.f90 main.f90 function1.f90 subrouine1.f90

这样会自动生成文件a.out，这是一个可执行文件。


 如果想让生成的可执行程序有一个名字，将命令改成

 ifort -o execname global.f90 main.f90 function1.f90 subrouine1.f90

即可。

输入命令

 ./execname

即可执行该程序。

　　这个题的readme写了直接make就行。一开始make报错是因为没有load ifort，load完直接make就可以了。运行要等好久，我还以为卡了。。

（4）安装VTune：

教程：https://blog.csdn.net/WY_stutdy/article/details/79106501

注册试用版最多需要两天，还要等。