Trilinos 简介 Trilinos 简介卢朓 Trilinos 简介卢朓 Trilinos 简介 Trilinos 简介 Trilinos 的安装和使用 Trilinos 简介 Trilinos 简介 Trilinos 的安装和使用 Trilinos 简介 Trilinos Epetra

Size: px

Start display at page:

Download "Trilinos 简介 Trilinos 简介卢朓 Trilinos 简介卢朓 Trilinos 简介 Trilinos 简介 Trilinos 的安装和使用 Trilinos 简介 Trilinos 简介 Trilinos 的安装和使用 Trilinos 简介 Trilinos Epetra"

尤娴凌
5 years ago
Views:

1 Department of Science and Engineering Computing School of Mathematics School Peking University Department of Science and Engineering Computing School of Mathematics School Peking University August 14, 2007 August 14, 2007 Trilinos Packages Sandia 国家实验室的产品把 MPI 的 C 接口改成了 C++ 接口提供了 blas, lapack 的 C++ 接口各个软件包之间无缝的协作利用已有的软件提高效率屏蔽内部细节, 便于灵活使用三大模块 : ANA, LAL, APP 关键部分采用 Fortran, 提高计算性能 Trilinos 是一个大的软件项目, 包含了几十个软件包 Each package is Focus on important, state-of-the-art algorithms regime. Developed by a small team of domain experts. Self-contained: No explicit dependencies on any other software packages (with some special exceptions). Congurable/buildable/documented on its own. Sample packages: NOX, AztecOO, ML, IFPACK. Special package collection: Petra (, Tpetra, Jpetra): 具体数据对象 Thyra: 抽象概念界面 Teuchos: 一般工具 New_package: 帮助构造新软件包

2 Trilinos 的软件包 Basic Linear Algebra Libraries : Current Project C++ Library Tpetra: Next Generation C++ Library Jpetra: Java Library 简单的 congure, make, make install 不太好 Linear Solver: Amesos: OO interface to 3rd party direct solvers AztecOO: Preconditioned Krylov Packaage based on Aztec Komplex: Complex solver via equivalent real formulations Belos: Next generation Krylov and block Krylov solvers 选择一下自己要用的包遇到不能找到 blas 或者 lapack 库的报错, 这可能是你指明库位置的 congure 语句需要改变一下使用时需要注意链接库的次序 Preconditioners ML: Multi-level preconditioners Meros: Segregated/Block Preconditioners IFPACK: Algebraic preconditioner 参考我的博客 Common Services Teuchos: parameter Lists, Blas Interface, etc Trilinos 的简单例子前面例子的 Makele 的一个简单的例子 int main(int argc, char * argv []){ MPICC = mpicxx CFLAGS = -DHAVE_MPI -DMPICH_IGNORE_CXX_SEEK LDFLAGS = -lepetr -lteuchos -llapack -lblas -lg2c int NumProc; int MyPID; OBJS = epetra_hello.o all: epetra_hello.exe MPI_Init(&argc, &argv); _MpiComm Comm(MPI_COMM_WORLD); %.o:%.cpp $(MPICC) -c $(CFLAGS) $< -o $@ NumProc = Comm.NumProc(); MyPID = Comm.MyPID(); epetra_hello.exe: $(OBJS) $(MPICC) $(OBJS) $(LDFLAGS) -o $@ cout << Comm <<endl; epetra_hello.o:epetra_hello.cpp MPI_Finalize(); return EXIT_SUCCESS; clean: rm -f *.o rm -f *.exe

3 _Import and _Export provides 向量是基本的数据结构,trilinos 软件主要通过构造向量两种类型的数据 :double 和 int serial Vector or distributed ( 分布式 ) Vector Serial : usually small distrubted vectors: signicantly larger 如何分布 : 使用通信器和 _Map A map is basically a partitioning of a list of global IDs. The fundamental communicator object, _Comm The _Map object The creation and assembly of vectors. A tool to redistributing vectors across processes. _Import and _Export provides 纯虚类 _Comm 两个子类, 一个是提供串行通信, 一个提供并行通信使用 _Comm 的例子 I 使用 _Comm 的例子 II _Import and _Export provides #include ``_ConfigDefs.h'' #ifdef HAVE_MPI // 在编译的时候加选项 -DHAVE_MPI #include ``mpi.h'' #include ``_MpiComm.h'' #else #include ``_SerialComm.h'' #endif // 加上其他的头文件 int main(int argc, char *argv[]) #ifdef HAVE_MPI MPI_Init(&argc, &argv); _MpiComm Comm(MPI_COMM_WORLD); #else _SerialComm Comm; #endif _Import and _Export provides // 其他代码 #ifdef HAVE_MPI MPI_Finalize(); #endif

4 _Import and _Export provides 串行和并行的代码基本相同 _Comm 的方法和 MPI 函数 _Comm 的方法有 MyPID(), NumProc( ), Barrier( ), Broadcast( ), SumAll( ), GatherAll( ), MaxAll( ), MinAll( ), ScanSum( ) 例 int NumProc = Comm.NumProc(); // 处理器的个数 int MyPID = Comm.MyPID(); // 处理器的序号 (rank) 阅读类的头文件 _Import and _Export provides 多个 CPU 分布的一组编号称作一个 Map, (_Map 类 ) 定义了 NumGlobalElements( ): The total number of elements across all processes. NumMyElements(): The number of elements on the calling process. MinAllGID( )/MaxAllGID( ): The minimum/maximum global index value across all processes. MinMyGID( )/MaxGID( ): The minimum/maximum global index value on the calling process. MinLID( )/ MaxLID( ): The minimum/maximum local index value on the calling process. LinearMap( ): Return true if the elements are distributed linearly across processes, i.e., process 0 gets the frist n/p elemetns, process 1 gets the next n/p elements, etc, n is the number of elements and p is the number of processes. I II #include ``mpi.h'' #include ``_MpiComm.h'' #include ``_Map.h'' # int main(int argc, char * argv[]){ MPI_Init(&argc, &argv); _MpiComm comm(mpi_comm_world); if(map1.linearmap()) cout << ``Map1 is a linear Map''<<en else cout <<``Map1 is not a linear Map''<<endl; cout << Map1; // 构造 Map 的第二种方法 : 给定每个 process 上的单元数 // 构造一个线性 Map _Import and _Export provides // 构造 Map 的第一种方法 : 给定总的单元数, 构造一个线性 Map int NumGlobalPoint=5; _Map Map1(NumGlobalPoint,0,comm); //comm 前面的 0 是全局点编号开始的位置 // 如果用 3 个处理器, 每个处理器分布几个 Global IDs _Import and _Export provides int NumMyPoints = comm.mypid(); _Map Map2(-1,NumMyPoints,0,comm); //Map1 是 linear Map 吗?

5 构造 Map 的第三种方法 _Import and _Export provides _Import and _Export provides 每个进程上都定义了一个长度为 Length 的 double 向量 #include ``_SerialDenseVector.h'' _SerialDenseVector DoubleVector(Length); 每个进程上都定义了一个长度为 Length 的 integer 向量 #include ``_SerialDenseVector.h'' _SerialIntDenseVector IntVector(Length); 不用手动删除内存 [],() 都是可以取值和赋值的方法 DoubleVector[0] = 10.0; // 赋值 DoubleVector(Length-1) = 20.0; // 赋值 cout << IntVector[Length-1]<<endl; // 取值 ( ), 强制越界检查 ; [] 无强制越界检查 Figure: 指定每个 process 上的单元个数和 GID, 构造非线性 Map 把 _Vector 转换成指针使用 _Import and _Export provides A distributed object is an entity whose elements are partitioned across more than one process. 's distributed objects are created from a Map. 使用 _Map 的构造函数例 : _Vector x(map) Copy constructor. 例 : _Vector y(x); 例 : _Vector x(copy,map,localvalues); Localvalues is a pointer to an array of double precision values. 两种模式 : Copy: 分配空间, copy data View: 不分配空间, 使用 LocalValues 所指空间 _Import and _Export provides double * x_values; x.extractview( & x_values); for(int i=0;i<x.mylength();++i)x_values[i]*=10; x_values 指针指向 _Vector x 的数据存储空间对 x_values[i] 的操作等同于对 x[i] 的操作 x.extraccopy(&x_values); 和 ExtractView 有什么区别?

_Vector 的其它方法 _Import and _Export int Norm1(double * Result) const 返回 1-norm: P n i=1 jxi j int Norm2(double * Result) const 返回 2-norm: p P ni=1 jxi j2 _Import and _Export apply o-processor

6 _Vector 的其它方法 _Import and _Export int Norm1(double * Result) const 返回 1-norm: P n i=1 jxi j int Norm2(double * Result) const 返回 2-norm: p P ni=1 jxi j2 _Import and _Export apply o-processor communication. They construct a communication plan that can be called repeatedly _Import and _Export provides _Import and _Export provides 构造函数的参数 :SourceMap, TargetMap _Import 中 TargetMap 的 GID 是唯一的, 不会同一个 GID 出现在两个 Processor 上 _Export 中 SourceMap 的 GID 是唯一的, 不会同一个 GID 出现在两个 Processor 上 Export 的例子 I Export 的例子 II int main(int argc, char *argv[]) { //_Export 类 // 这个程序要用两个节点运行 #include "_ConfigDefs.h" #ifdef HAVE_MPI #ifdef HAVE_MPI MPI_Init(&argc, &argv); _MpiComm Comm(MPI_COMM_WORLD); #else _SerialComm Comm; _Import and _Export provides #include "mpi.h" #include "_MpiComm.h" #else #include "_SerialComm.h" #endif #include "_Map.h" #include "_Vector.h" #include "_Export.h" _Import and _Export provides #endif // 总的单元数 int NumGlobalElements = 4; // 局部的单元数 int NumMyElements; _IntSerialDenseVector MyGlobalElements;

7 Export 的例子 III Export 的例子 IV _Import and _Export provides if( Comm.MyPID() == 0 ) { NumMyElements = 3; MyGlobalElements.Size(NumMyElements); MyGlobalElements[0] = 0; MyGlobalElements[1] = 1; MyGlobalElements[2] = 2; else { NumMyElements = 3; MyGlobalElements.Size(NumMyElements); MyGlobalElements[0] = 1; MyGlobalElements[1] = 2; MyGlobalElements[2] = 3; //Map 构造函数的第一个参数 -1, // 表示第二个参数给出的是局部向量元素的个数, _Import and _Export provides // 第三个参数是一个指针类型的, 给出 // 局部序号对应的整体序号根据这些信息, // 总的向量元素的个数就可以计算出来 // 下面的 0 表示向量元素的序号从 0 开始计算 // 这些整体序号可能是重复, // 也就是该元素在多个处理器重复保存 _Map Map(-1,MyGlobalElements.Length(), MyGlobalElements.Values(),0, Comm); // create a vector based on map _Vector x(map); for( int i=0 ; i<nummyelements ; ++i ) x[i] = 10*( Comm.MyPID()+1 ); cout << x; // 构造一个目标 Map Export 的例子 V Export 的例子 VI int NumMyElements_target; if( Comm.MyPID() == 0 ) NumMyElements_target = NumGlobalElements; else NumMyElements_target = 0; _Map TargetMap(-1,NumMyElements_target,0,Comm); _Export Exporter(Map,TargetMap); // 元所在的处理器, 并且把数值加起来作为 // 结果存在目标向量里 y.export(x,exporter,add); cout << y; #ifdef HAVE_MPI MPI_Finalize(); #endif _Import and _Export provides // work on vectors _Vector y(targetmap); // 这个操作就会把在各个处理器上的和 TargetMap // 的 GID( 整体序号 ) 相同的元素传输到到目标单 _Import and _Export provides return(exit_success);

8 MPI_Barrier and MPI_Bcast 归约操作 : MPI_Reduce _Import and _Export MPI_Barrier 用于函数同步所有的进程都调用了此函数后才返回 MPI_Barrier(MPI_COMM_WORLD); comm.barrier() MPI_Bcast 将自己 buer 中的内容发送给通信器中所有其他进程 double dvalue; int root = 0; int MPI_Bcast(&dvalue, 1, MPI_DOUBLE, root, MPI_COMM_WORLD); comm.broadcast(&dvalue, 1, root); _Import and _Export int MPI_Reduce(void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm) Input Parameters sendbuf address of send buffer (choice) count number of elements in send buffer datatype data type of elements of send buffer op reduce operation (handle) root rank of root process (integer) comm communicator (handle) Output Parameter provides provides recvbuf address of receive buffer Algorithm: This implementation currently uses a simple tree algorithm. MPI_Reduce 的例子 :sum.cpp I MPI_Reduce 的例子 :sum.cpp II /* FileName: sum.cpp 计算 sum = 编译命令 : mpicxx -DMPICH_IGNORE_CXX_SEEK -o sum.exe sum.cpp 这个程序用的进程个数必须能整除 1000 */ int sum,totalsum; int imsg[4]; double a[9]; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&NumProc); MPI_Comm_rank(MPI_COMM_WORLD,&MyPID); _Import and _Export #include <iostream> #include "mpi.h" using namespace std; _Import and _Export startval = 1000*MyPID/NumProc + 1; endval = 1000*(MyPID+1)/NumProc; provides int main(int argc,char *argv[]) { int NumProc, MyPID, startval, endval,i; provides sum =0; for (i=startval; i< endval; i++) sum=sum + i;

$cpp III MPI_AllGather //recvbuf 和 sendbuf 必须不同 MPI_Reduce(&sum,&totalsum,1,MPI_INT,MPI_SUM, 0,MPI_COMM_WORLD); if( MyPID == 0) cout<<"the sum of i from 1 to 1000 \ is "<< totalsum <<endl;$

9 MPI_Reduce 的例子 :sum.cpp III MPI_AllGather //recvbuf 和 sendbuf 必须不同 MPI_Reduce(&sum,&totalsum,1,MPI_INT,MPI_SUM, 0,MPI_COMM_WORLD); if( MyPID == 0) cout<<"the sum of i from 1 to 1000 \ is "<< totalsum <<endl; MPI_Finalize(); double *dvalues, *dvalues = new double[numproc]; MPI_Allgather(dvalues, 1, MPI_DOUBLE, dvalues2, NumProc, MPI_DOUBLE,MPI_COMM_WORLD); Comm.GatherAll(dvalues, dvalues2, 1); _Import and _Export return 0; _Import and _Export provides provides Allgather 示意图 _Import and _Export provides _Import and _Export provides Calling FORTRAN from C++ varies from platform to platform _ handles the details of this issue

10 不使用的程序 I int main(int argc, char * argv[]) { 计算向量 ~a ~ b 分配向量 MPI_Init(&argc, &argv); int rank, size; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &size); 计算每个处理器上的局部内积使用 MPI 函数 MPI_Allreduce 把局部内积相加得到 Global inner product. if(size!= 2){ cout <<"Please use two processors"<<endl; MPI_Abort(MPI_COMM_WORLD,EXIT_FAILURE); int length; if(rank == 0 ) length = 3; if(rank == 1 ) length = 2; 不使用的程序 II 不使用的程序 III double * a = new double[length]; double * b = new double[length]; cout<<"cpu NO. "<<rank<< ", b["<<i<<"] = "<<b[i]<<endl; int i; /// 赋值 for(i = 0; i < length; i ++){ a[i] = (rank+1)*10.0; b[i] = (rank+2)*10.0; /**** 计算内积的部分的头 *************/ double * c = new double; /// 储存局部向量的内积 c[0] = 0.0; for(i = 0; i < length; i ++) c[0] += a[i]*b[i]; double result; //int MPI_Allreduce(void*, void*, int, // MPI_Datatype, MPI_Op, MPI_Comm); for(i = 0; i < length; i ++) cout<<"cpu NO. "<<rank<< ", a["<<i<<"] = "<<a[i]<<endl; for(i = 0; i < length; i ++) MPI_Allreduce(c,&result,1,MPI_DOUBLE, MPI_SUM,MPI_COMM_WORLD); /**** 计算内积的部分的尾 *************/

11 不使用的程序 IV 使用计算内积的程序 I cout <<"on the cpu No. "<<rank<<", a dot b = "<<result<<endl; delete[]a; delete[]b; delete c; MPI_Finalize(); #include "_ConfigDefs.h" #include "mpi.h" #include "_MpiComm.h" #include "_Map.h" #include "_Vector.h" int main(int argc, char *argv[]) { return EXIT_SUCCESS; MPI_Init(&argc, &argv); _MpiComm Comm(MPI_COMM_WORLD); int rank, size; rank = Comm.MyPID(); size = Comm.NumProc(); 使用计算内积的程序 II 使用计算内积的程序 III if(size!= 2){ cout <<"Please use two processors"<<endl; MPI_Abort(MPI_COMM_WORLD,EXIT_FAILURE); int length; if(rank == 0) length = 3; if(rank == 1) length = 2; // 构造一个 Map /// 赋值 for(int i = 0; i < length; i ++){ a[i] = (rank+1)*10.0; b[i] = (rank+2)*10.0; cout << a; cout << b; _Map Map(-1,length,0,Comm); double result; // 构造两个向量 _Vector a(map); /// 计算 a 和 b 的内积 a.dot(b, &result); _Vector b(map); cout <<"a dot b = "<<result<<endl;

12 使用计算内积的程序 IV 的程序是怎么写的? MPI_Finalize(); return EXIT_SUCCESS; 还是利用 MPI_Allreduce 函数利用计算局部内积利用 C++ 的函数重载, 使得不用指明数据类型重载了输出算符实现 Dot 的部分源程序 I 实现 Dot 的部分源程序 II //========================================================================= double **A_Pointers = A.Pointers(); // 程序是对 _Multivector 写的, _Vector // 是其子类 int _MultiVector::Dot(const _MultiVector& A, double *Result) const { // Dot product of two Multi int i; if (Num_!= A.Num()) EPETRA_CHK_ERR(-1); if (MyLength_!= A.MyLength()) EPETRA_CHK_ERR(-2); UpdateDoubleTemp();// 分配一个临时 double 变量 //DOT 是对库函数的调用的接口 for (i=0; i < Num_; i++) DoubleTemp_[i] = DOT(MyLength_, Pointers_[i], A_Pointers[i]); Comm_->SumAll(DoubleTemp_, Result, Num_); UpdateFlops(2*GlobalLength_*Num_); return(0); //======================================================= //===SumAll 函数的重载 (double 型 )===========================

13 实现 Dot 的部分源程序 III 实现 Dot 的部分源程序 IV // 就是调用 MPI_Allreduce( MPI_SUM ) 的接口 int _MpiComm::SumAll(double * PartialSums, double * GlobalSums, int Count) const { EPETRA_CHK_ERR(MPI_Allreduce(PartialSums, GlobalSums,Count, MPI_DOUBLE, MPI_SUM, MpiCommData_->Comm_)); return(0); //===SumAll 函数的重载 (int 型 )================================= int _MpiComm::SumAll(int * PartialSums, int * GlobalSums, int Count) const { EPETRA_CHK_ERR(MPI_Allreduce(PartialSums, GlobalSums,Count, MPI_INT, MPI_SUM, MpiCommData_->Comm_)); return(0); //======================================================= 1 参考我博客上的安装方法或者 Trilinos 的说明文档安装 Trilinos, 不要使用 enable-ml 安装 ml 包 2 仿照 lecture 上例子练习 _Map, 分别使用三种构造函数构造 _Map 的实例, 回答课件上提出的问题 3 仿照利用 _Vector 计算向量内积 ( 方法 Dot) 的例子, 学习使用其他 _Vector 的方法, 如 Norm2, Scale 等可以参考的帮助文档 : 4 课件使用 MPI_Reduce 计算 1 到 1000 的和, 请换成 MPI_Allreduce 函数计算, 并说明 MPI_Reduce 和 MPI_Allreduce 的区别

Department of Science and Engineering Computing School of Mathematics School Peking University August 14, 2007

Department of Science and Engineering Computing School of Mathematics School Peking University August 14, 2007 Department of Science and Engineering Computing School of Mathematics School Peking University

Trilinos 简介 Trilinos 简介 卢朓 Trilinos 简介 卢朓 Trilinos 简介 Trilinos 简介 Trilinos 的安装和使用 Trilinos 简介 Trilinos 简介 Trilinos 的安装和使用 Trilinos 简介 Trilinos Epetra

Trilinos 简介 Trilinos 简介卢朓 Trilinos 简介卢朓 Trilinos 简介 Trilinos 简介 Trilinos 的安装和使用 Trilinos 简介 Trilinos 简介 Trilinos 的安装和使用 Trilinos 简介 Trilinos Epetra