# YHSCALAPACK

**Repository Path**: e-level-parallel-algorithm/yhscalapack

## Basic Information

- **Project Name**: YHSCALAPACK
- **Description**: 基于稠密线性系统开源软件SCALAPACK开发了异构融合HU-SCALAPACK，本算法库提供异构计算接口，支持GPU、MIC、MATRIX2000等加速卡。
- **Primary Language**: Unknown
- **License**: Apache-2.0
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2021-02-24
- **Last Updated**: 2021-06-22

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# YHSCALAPACK

#### 介绍
基于稠密线性系统开源软件SCALAPACK开发了异构融合HU-SCALAPACK，本算法库提供异构计算接口，支持GPU、MIC、MATRIX2000等加速卡。

#### 软件架构
近年来，异构系统硬件飞速发展，为了解决相应的编程和执行效率问题，异构并行编程已被广泛使用和研究。Scalapack是基于分布式存储的并行计算库，不支持异构并行计算，而HU-Scalapack是在Scalapack的基础上加入了对异构并行的支持，部分函数能够在具有异构加速器平台上进行异构并行。HU-Scalapack集成了MIC offload、OpenMP target、COI和SCIF的异构方式，能够根据不同类型的计算任务选择所需要的异构并行方法。其中MIC offload模式是在Intel Xeon Phi加速卡上计算的异构方式；OpenMP target模式是使用OpenMP的target方式来实现异构并行；COI模式指的是在Xeon Phi协处理器上使用COI的方式进行异构并行；SCIF模式是使用SCIF通信模式来让cpu和加速卡之间通信从而实现异构计算。

- **MIC offload的异构方式**

  MIC（Many Integrated Core）是Intel的集成众核架构协处理器，用作高性能计算的超级计算机或服务器加速卡。Offload是MIC语言扩展里最基本的关键字，是MIC程序的根本，其基本作用是在offload作用范围内的程序代码在MIC卡上运行。Scalapack中BLAS负责基本线性代数运算，例如矩阵-矩阵乘、矩阵向量乘、向量-向量乘等。Scalapack求解线性方程组最终需要调用BLAS库中的基本线性代数函数，而且这部分计算量在求解线性方程组中占了很大的一部分，所以Scalapack的MIC offload异构方式的整体思路是：把Scalapack中计算最耗时的部分（即BLAS库中的矩阵乘法计算部分）通过MIC offload的方式加载到MIC卡中计算，充分发挥MIC卡多核的优势加速运算。具体实现方式是：针对BLAS库的LEVEL3矩阵-矩阵乘法函数，因为这些矩阵-矩阵运算的函数是三重循环的结构，只需要在循环语句前添加编译指导语句即可以将这部分的计算加载到MIC中。为了发挥MIC卡的众核优势，在MIC上采用了OpenMP多线程的方法实现并行化。

- **OpenMP target的异构方式**

  MIC offload异构方式只能针对于Intel众核加速器Xeon Phi，具有比较大的局限性。OpenMP从4.0版本（2013年发出）之后开始支持异构编程(CPU+加速器的异构系统)，可以实现在不同的加速卡上异构并行，通用性强。OpenMP的异构并行方式与MIC的异构并行方式类似，通过添加编译指导语句方式将代码块加载到加速卡上运行，其基本语法结构如下：

  ```fortran
  !$omp target clause [[,]  clause],…]
  	structured-block
  !$omp end target
  device(scalar-integer-expression)
  map(alloc | to |  from | tofrom: list)  if(scalar-expr)
  ```

  OpenMP target异构的具体实现方式与MIC offload类似：在BLAS库中针对LEVEL3矩阵-矩阵乘法函数中的循环语句块前添加omp target编译指导语句即可以将这部分的计算加载到MIC中，并且使用OpenMP多线程并行的方式实现在加速卡上开启多线程并行计算。

- **SCIF异构并行方式**

  SCIF(Symmetric Communications Interface) ，对称通信接口，是基于主处理器与协处理器之间的通信工具。SCIF的相较于COI与OpenMP等实现更为底层，比较接近硬件，具有更低的通信延迟和更高的通信带宽，使用SCIF可以有效的减少主处理器与协处理器之间的通信开销。

  SCIF实现包含两个部分：用户态的接口函数（User mode library）和内核态的接口函数（Kernel Mode library），两者接口相似的函数一共有19个，主要实现部分为基于driver、library做封装。上层软件全部基于用户态的接口函数实现。因此实际上用户在利用SCIF进行通信优化的时候只需要调用用户态的接口函数。SCIF的通信过程是指信息的传递过程，是两个endpoint之间的信息传递，一个endpoint可以监听另一个与他建立对应连接的endpoint，而另一个endpoint则等待这个连接请求。两个endpoint之间要想进行数据传输，就必须先建立连接，建立endpoint的过程类似于socket网络编程。通信的过程大致如下：生成新的endpoint(open)、绑定地址(bind)、建立连接(connect)、监听(listen)、接受连接请求(accept)、数据传输。SCIF需要两套代码来实现数据传递，一套在Host端执行，另一套需拷贝到Device端执行。

  HU-Scalapack SCIF的异构并行方式可以利用SCIF的低通信延迟和高通信带宽的优势，提高Host端和Device端之间的数据传输效率，提高计算效率。HU-Scalapack SCIF异构实现方式是通过在Host端和Device端建立通信连接后，使用scif_send()、scif_recv()或者RMA函数进行数据传输，而在Device端使用OpenMP多线程的方式并行。

-  **COI异构并行方式**

  COI是MIC架构下的分载模式的一个库，MIC虽然提供了通过简单的编译指导语句（#pragma offload）的方式来将部分代码和数据加载到Xeon Phi上进行计算。但这样的方式太过简单，且自由度较小。但如果使用COI，用户可以获取更大的自由度，包括控制CPU和MIC的同步，控制MIC卡上程序的创建和退出；起止端的异步操作；起止端的数据缓冲。为开发更加灵活的MIC程序提供了便利。

  利用COI编写的程序在编译的时候就会生成两个可执行文件。一个在CPU上执行，一个在MIC上执行。启动程序只需要调用CPU端的可执行文件即可。系统就会在需要的时候把MIC端可执行文件和所需要的库发送到MIC端并执行MIC端的程序。COI的指令和数据传输通过COIPipeline函数实现，COIPipeline类似于RPC（Remote Procedure Call）的机制，可以将一系列指令序列插入到COIPipeline中，这些指令可以顺序地在sink端执行，这里的指令序列主要是要在远程调用的函数序列。在COIPipeline中插入的函数会在sink端按序执行，因为可以在远端调用完整函数，所以有了比单一offload更多的自由度。COIPipeline在插入函数时除了插入该函数需要的参数外，还可以传递一块buffer，实现从source端到sink端的数据传输。

  HU-Scalapack COI的异构并行是通过COIPipeline函数把BLAS的矩阵-矩阵乘法函数从source端传输到sink端并在sink端执行，另外sink端所需要的数据也是通过COIPipeline函数的buffer块传输，按照这种方式source端和sink端之间的指令和数据通信得以实现，与其他方式一样，在sink端的矩阵乘法并行计算是通过OpenMP多线程方式来实现的。


#### 安装教程

- **HU-Scalapack MIC offload和OpenMP targets异构版本的编译安装**

1. 进入HU-Scalapack/lapack/BLAS目录下，编辑common.inc文件，在#define后面选择MIC_OFFLOAD或者OPENMP_TARGET。

   ```shell
   vim common.inc
   ```

   <img src="D:\HU-scalapack\common.png" alt="image-20210622160014524" style="zoom:80%;" />

2. 进入HU-Scalapack下的lapack文件夹

   ```shell
   $cd ./HU-Scalapack/lapack
   ```

   复制make.inc.example文件并重命名为make.inc：

   ```shell
   $cp make.inc.example make.inc
   ```

   打开make.inc文件，查看配置文件以及最后所生成的库文件：

   ```shell
   $vim make.inc
   ```

   修改BLAS库和LAPACK库的输出路径和名称：

   ```shell
   BLASLIB			= ../../miclibblas.a
   CBLASLIB		= ../../miclibcblas.a
   LAPACKLIB       = micliblapack.a
   TMGLIB         = miclibtmglib.a
   ```

   对源码进行编译：

   ```shell
   $make
   ```

   最后生成以下几个库文件：

   > miclibblas.a：BLAS库安装成功的库文件。
   >
   > micliblapack.a：LAPACK库安装成功的库文件
   >
   > miclibtmglib.a

   返回上一层目录，进入HU-Scalapack目录：

   ```shell
   $cd ..
   ```

   复制SLmake.inc.example，重命名为SLmake.inc：

   ```shell
   $cp SLmake.inc.example SLmake.inc
   ```

   编辑SLmake.inc文件，修改BLAS库文件路径和LAPACK库文件路径，设置为安装LAPACK时生成的库文件路径：

   ```shell
   $vim SLmake.inc
   ```

   ```shell
   "BLASLIB       = /your_directory/HU-scalapack/lapack/miclibblas.a"
   "LAPACKLIB     = /your_directory/HU-scalapack/lapack/micliblapack.a"
   "LIBS          = $"("LAPACKLIB")" $"("BLASLIB") 
   ```

   再对源码进行make编译：

   ```shell
   $make
   ```

   最后生成Scalapack库文件libScalapack.a，安装完成。

- **HU-Scalapack COI安装**

  进入HU-Scalapack/lapack/BLAS目录，打开并编辑common.inc文件，在#define后面选择COI。

  ```shell
  $vim common.inc
  ```

  <img src="D:\HU-scalapack\3.2.png" style="zoom:80%;" />

  设置完毕后，进入HU-Scalapack下的COI目录。该目录下是COI的level3矩阵乘函数程序，其中包括了单精度sgemm、双精度dgemm、单精度复数cgemm和双精度复数zgemm的矩阵乘函数：

  ```shell
  $cd HU-Scalapack/COI
  ```

  编译COI源码：

  ```shell
  $make
  ```

  如果编译过程中出现了找不到COI的库文件，分别进入每个coi_*gemm目录，修改makefile下的库路径，其中包括COI和OpenBlas库。

  <img src="D:\HU-scalapack\coi.png" style="zoom:80%;" />

  编译完成后，./bin目录下会有相应的可执行文件，该文件是在加速卡上运行的矩阵乘计算程序，在CPU端不能运行：

  > coi_cgemm_mic
  >
  > coi_dgemm_mic
  >
  > coi_sgemm_mic
  >
  > coi_zgemm_mic

  编译完成后，./COI目录下生成了一个库文件：libcoigemm.a，当需要的时候可以调用这个库来编译程序。编译完成了COI的库文件后，依照3.1小节的编译方法编译安装Lapack和Scalapack即可完成该版本的安装。

- **HU-Scalapack SCIF安装**

  进入HU-Scalapack/lapack/BLAS目录下，编辑common.inc文件，在#define后面选择SCIF。

  ```shell
  $vim common.inc
  ```

  设置完毕后，进入HU-Scalapack下的SCIF_GEMM目录。该目录是COI的level3矩阵乘函数，其中包括了单精度sgemm、双精度dgemm、单精度复数cgemm和双精度复数zgemm的矩阵乘函数：

  > scif_sgemm
  >
  > scif_dgemm
  >
  > scif_cgemm
  >
  > scif_zgemm

  编译SCIF源码:

  ```shell
  $make
  ```

  如果编译过程中出现找不到scif的库文件，分别进入每个scif_*gemm文件夹，修改makefile下的库路径：

  <img src="D:\HU-scalapack\scif.png" style="zoom:80%;" />

  编译完成后，./bin目录下会有相应的SCIF可执行文件，该可执行文件是在加速卡上运行的程序，在CPU端不能执行：

  > ```
  > scif_cgemm_mic
  > 
  > scif_dgemm_mic
  > 
  > scif_sgemm_mic
  > 
  > scif_zgemm_mic
  > ```

  编译完成了SCIF的库文件后，依照3.1小节的编译方法编译安装即可完成。编译完成后，./ SCIF_GEMM目录下生成了一个库文件：libscifgemm.a，当需要的时候可以调用这个库来编译程序。

#### 使用说明

- **调用HU-Scalapack库**

  HU-Scalapack库的接口和调用方式与Scalapack一致，在用户自己的程序中调用Scalapack程序需要四个基本步骤：

  1.  初始化进程网格

     一台抽象的并行计算机的*P*个进程可以表示为一维数组，通常将这个一维数组映射到二维矩形网格可以更方便地表示算法，这个网格称为进程网格。在程序的开始用户需要初始化进程网格，得到默认的系统上下文。用户也可以查询进程网格以识别每个进程的坐标。

     调用Scalapack TOOLS程序SL_INIT初始化进程网格。这个程序使用进程行优先排序初始化一个`PrxPc`（在源代码中以`NPROWxNPCOL`表示）进程网格，得到默认的系统上下文。用户可以通过调用BLACS_GRIDINFO查询进程网格以识别每个进程的坐标(MYROW*,* MYCOL)。

     完成这项任务的典型代码如下：

     ```fortran
     CALL BLACS_GET( -1,0,ICTXT )
     CALL BLACS_GRIDINIT( ICTXT,'Row-major',NPROW,NPCOL )
     CALL BLACS_GRIDINFO( ICTXT,NPROW,NPCOL,MYROW,MYCOL )
     ```

  2.   将矩阵分布到进程网格上

     在调用Scalapack程序之前所有全局矩阵必须分布在进程网格上。执行数据分布是用户的责任。每个要分布在进程网格上的全局矩阵都必须分配一个数组描述符，数组描述符通过调用Scalapack TOOLS的程序DESCINIT可以很简单地初始化，它必须在调用Scalapack程序之前设置。

     矩阵的数组描述符用以下代码分配：

     ```fortran
     CALL DESCINIT( DESCA,M,N,MB,NB,RSRC,CSRC,ICTXT, MXLLDA,INFO )
     CALL DESCINIT( DESCB,N,NRHS,NB,NBRHS,RSRC,CSRC, ICTXT,MXLLDB,INFO )
     CALL DESCINIT( DESCX,N,NRHS,NB,NB,0,0,ICTXT, MAX(1,NP),INFO ) 
     ```

     之后可以使用如下的代码调用从数据文件中读入矩阵的数据并将其分布到进程网格上：

     ```fortran
     CALL PDLAREAD( 'SCAEXMAT.dat',MEM(IPA),DESCA,0,0,MEM(IPW)  )
     CALL PDLAREAD( 'SCAEXRHS.dat',MEM(IPB),DESCB,0,0,MEM(IPW)  )
     ```

  3. 调用Scalapack程序

     当数据正确分布到进程网格上之后，用户就可以调用相应的驱动、计算、辅助程序进行需要的计算。各个程序的调用参数的详细解释可以在Scalapack源代码的注释中找到。

     下面的例子调用求解线性方程组*AX*＝*B*的简单驱动程序PDGESV：

     ```fortran
     CALL PDGESV( N,NRHS,MEM(IPA),1,1,DESCA,MEM(IPPIV),MEM(IPB),1,1,DESCB,INFO )
     ```

  4.  释放进程网络

     在进程网格上执行完计算后，通过调用BLACS_GRIDEXIT释放进程网格。当所有计算完成后，程序通过调用BLACS_EXIT退出。

     完成这个步骤的典型代码是：

     ```fortran
     CALL BLACS_GRIDEXIT( ICTXT )
     CALL BLACS_EXIT( 0 )
     ```

     在编译链接程序时，需要指定libScalapack.a、micliblapack.a、miclibblas.a文件的位置，确保程序可以被正确链接。如果是COI异构模式，则编译的时候需要指定libcoigemm.a库位置；如果是SCIF异构模式，在编译的时候需要指定libscifgemm.a库文件位置。编译结束后，当需要执行可执行文件时，若是COI异构模式，则需要把COI/bin目录下的文件复制到可执行文件当前目录下，以供程序传输到加速器端；若是SCIF异构模式，则需要把SCIF_GEMM/bin目录下的文件复制到可执行文件当前目录下，否则程序运行时会报错，提示找不到文件。

- **调用HU-Scalapack下COI或者SCIF矩阵乘库**

  HU-Scalapack/COI/libcoigemm.a是COI异构并行的矩阵乘法库函数，HU-Scalapack/SCIF_GEMM/libscifgemm.a是SCIF异构并行的矩阵乘法库函数，这两个库函数可以分别提供COI和SCIF方式异构并行的矩阵乘运算，具体调用方法如下。

  1. libcoigemm.a库的调用

     libcoigemm.a库分别提供了4种不同精度的运算：coi_sgemm（单精度）、coi_dgemm（双精度）、coi_cgemm（单精度复数）、coi_zgemm（双精度复数）。其中API分别为：

     ```c
     void coi_sgemm_(int *Order,int *TransA,
                int *TransB,int *M,int *N,int *K,float *alpha,float *A,
                int *lda,float *B,int *ldb,float *beta,float *C,int *ldc,int *mic_num)
     void coi_dgemm_(int *Order,int *TransA,
                 int *TransB,int *M,int *N,int *K,double *alpha,double *A,
                 int *lda,double *B,int *ldb,double *beta,double *C,int *ldc,int *mic_num)
     void coi_cgemm_(int *Order,int *TransA,
                 int *TransB,int *M,int *N,int *K,float *alpha,float *A,
                 int *lda,float *B,int *ldb,float *beta,float *C,int *ldc,int *mic_num)
     void coi_zgemm_(int *Order,int *TransA,
                 int *TransB,int *M,int *N,int *K,double *alpha,double *A,
                 int *lda,double *B,int *ldb,double *beta,double *C,int *ldc,int *mic_num)
     
     ```

     > 其中：Order候选值有CblasRowMajor和CblasColMajor，这两个参数决定了一维数组怎样存储在内存中，一般选用CblasRowMajow。
     >
     > TransA和TransB：表示矩阵A和B是否进行转置，候选参数CblasTrasn，CblasNoTrans。
     >
     > 参数M：表示A或C的行数。
     >
     > 参数N：表示B或C的行数。
     >
     > 参数K：表示A的列数或者B的行数。
     >
     > 参数LDA：表示A的列数，与转置与否无关。
     >
     > 参数LDB：表示B的列数，与转置与否无关。
     >
     > 参数LDC：始终等于N。
     >
     > 参数alpha：表示矩阵A的系数$\alpha$。
     >
     > 参数beta：表示矩阵C的系数$\beta$。
     >
     > 参数mic_num：表示可执行文件运行在mic端的编号。

     以单精度coi_sgemm函数为例，下面的程序coi_test.c展示如何使用该库函数，该程序的源代码为：

     ```c
     #include <stdio.h>
     #include <cblas.h>"
     #include <intel-coi/source/COIProcess_source.h>"
     #include <intel-coi/source/COIEngine_source.h>"
     #include <intel-coi/source/COIPipeline_source.h>"
     #include <intel-coi/source/COIBuffer_source.h>"
     #include <intel-coi/source/COIEvent_source.h>"
     #include <intel-coi/common/COIMacros_common.h>"
     
     int main"()
     {
         int size = 10;
         float A[size][size];
         float B[size][size];
         float C[size][size];
         for (int i = 0; i < size; i++)
         {
             for (int j = 0; j < size; j++)
             {
                 A[i][j] = 2.0 ;
                 B[i][j] = 1.0 ;
                 C[i][j] = 0.0;
         } 
     
     float alpha = 1.0;"
         float beta = 0.0;"
     	int mic_num = 1;"
     	OPENBLAS_CONST enum CBLAS_ORDER LAYOUT = CblasColMajor;
     	OPENBLAS_CONST enum CBLAS_TRANSPOSE isTransA = CblasNoTrans;
     	OPENBLAS_CONST enum CBLAS_TRANSPOSE isTransB = CblasNoTrans;
     	printf("CblasColMajor = %d\n",CblasColMajor);
      	printf("CblasNoTrans = %d\n",CblasNoTrans);
     	printf("CblasRowMajor = %d\n",CblasRowMajor);	
         coi_sgemm_(&LAYOUT ,&isTransA,&isTransB,
                     &size,&size,&size,
                     &alpha,(float *)A,&size,
                     (float *)B,&size,&beta,
                     (float *)C,&size,&mic_num);  
     
     	for (int i = 0; i < size; i++)
         {
             for (int j = 0; j < size; j++)
             {            
                 printf("C[%d][%d] = %f\n",i,j,C[i][j]);
             }
         }
     } 
     
     ```

  2. libscifgemm.a库的调用

     ibscifgemm.a库分别提供了4种不同精度的运算：scif_sgemm（单精度）、scif_dgemm（双精度）、scif_cgemm（单精度复数）、scif_zgemm（双精度复数），其中API分别为：

     ```c
     void scif_sgemm_(int *Order,int *TransA,int *TransB,int *M,int *N,int *K,float *alpha,
       float *A,int *LDA,float *B,int *LDB,float *BETA,const float *C,int *LDC) 
     void scif_dgemm_(int *Order,int *TransA,int *TransB,int *M,int *N,int *K,double *alpha,
       double *A,int *lda,double *B,int *ldb,double *beta,double *C,int *ldc)
     void scif_cgemm_(int *Order,int *TransA,int *TransB,int *M,int *N,int *K,const void *alpha,
       const void *A,int *lda,const void *B,int *ldb,const void *beta,const void *C,int *ldc)
     void scif_zgemm_(int *Order,int *TransA,int *TransB,int *M,int *N,int *K,const void *alpha,
        const void *A,int *lda,const void *B,int *ldb,const void *beta,const void *C,int *ldc) 
     ```

     > 其中：Order候选值有CblasRowMajow和CblasColMajor，这两个参数决定了一维数组怎样存储在内存中，一般选用CblasRowMajow。
     >
     > TransA和TransB：表示矩阵A和B是否进行转置，候选参数CblasTrasn，CblasNoTrans。
     >
     > 参数M：表示A或C的行数。
     >
     > 参数N：表示B或C的行数。
     >
     > 参数K：表示A的列数或者B的行数。
     >
     > 参数LDA：表示A的列数，与转置与否无关。
     >
     > 参数LDB：表示B的列数，与转置与否无关。
     >
     > 参数LDC：始终等于N。
     >
     > 参数alpha：表示矩阵A的系数$\alpha$。
     >
     > 参数beta：表示矩阵C的系数$\beta$。
     
     以单精度scif_sgemm为例，下面的程序scif_test.c展示如何调用该函数，该程序的源代码为：
     
     ```c
     #include <stdio.h>
     #include <stdlib.h>
     #include <math.h>
     #include "cblas.h"
     #define SIZE 16"
     int main(void)
     {
     	int Order    = CblasColMajor; 
     	int isTransa = CblasNoTrans;
     	int isTransb = CblasNoTrans;
     	int M = SIZE;
     	int N = SIZE;
     	int K = SIZE;
     	int lda;
     	int ldb;
     	int ldc;  
       
      	ldb = max(K,N);
     	ldc = max(M,N);
     	float *a = (float *)malloc(M * K * sizeof(float));
     	float *b = (float *)malloc(K * N * sizeof(float));
     	float *c = (float *)malloc(M * N * sizeof(float));
     	float *c1= (float *)malloc(M * N * sizeof(float));
     	int len_i,len_j;
     	for( len_i = 0; len_i < M; len_i++)
     	{
     		for( len_j = 0; len_j < N; len_j++)
     	{
     		*(a  + len_i*N + len_j)= 2e0;
     		*(b  + len_i*N + len_j)= 1e0;
     		*(c  + len_i*N + len_j)= 1e0;
     		*(c1 + len_i*N + len_j)= 1e0;
     	}
     	} 
     	float alpha=1.0;
     	float beta =1.0;
     	printf("#  data has ready.  #\n");
     	scif_sgemm_(&Order,&isTransa,&isTransb,&M,&N,&K,&alpha,a,&lda,b,&ldb,&beta,c1,&ldc);
     	for(int i = 0; i < M; i++)
     	{
     		for(int j = 0; j < N; j++)
     		{	
     			printf("scif_c[%d][%d] = %lf\n",i,j,*(c1+i*N+j));
     		}
     	}
     	free(a);free(b);free(c);
     	return 0;
     }
     
     ```