UCGLE-master学习

0、简介

大致分成3部分，
1、对数据进行处理分析、判断是否适合经典的lsa算法。
2、如果适合就执行else{}部分（classicalGMRES(&v, &A)）。
3、不适合就执行if{}部分(项目作者写的）

1、数据分析

main()函数中定义一个结构体

com_lsa com;//一个mpi通信组

数据初始化函数

 non_lsa = mpi_lsa_init(argc,argv,&com);  // 初始化com，进程分组，给进程组属性区分GMRES_GR组和其他组

在初始化时候判断是否使用lsa

int argsHandleComponents(int argc, char ** argv, int * lsa_gmres, int * lsa_arnoldi)
//解析命令行argv，判断是否使用lsa 
//修改指针指向的值，不修改指针的值，把返回值传出去
//函数内部从到argc-1作比较(0代表路径)
//int 返回解析的状态
for(i=1;i<argc;i++){
if(strcmp(argv[i],"-lsa_gmres")==0){ //等于这个命令
  		cmpt1=1;
  		if(i+1<=argc-1){
  			i++;
  			*lsa_gmres=atoi(argv[i]);// 并对返回值进行修改，指针做返回值
  		}else{
  			if(rank==0)
  				printf("Warning: Number of nodes for GMRES isn't set, using 1 node\n");
  			*lsa_gmres=1;
  		}
  	}



if(cmpt1 && !cmpt2){ //匹配第一个命令，不匹配第二个，把第二个标记成1
  	*lsa_arnoldi=1;
  	cmpt2=1;
  	if(rank==0)  // 组内0号进程，打印提示信息
  		printf("Warning: Number of nodes for ARNOLDI isn't set, using 1 node\n");//

if(cmpt1 && cmpt2)  //正常返回0，因为前一个if-else把两个都变成了1
  	return 0;
  else 
  	return 1;  // 异常返回1

在初始化数据时创建组内通信组

int mpi_lsa_create_groups(com_lsa * com){ // 创建组内进程的通信组 

for(i=0;i<4;i++){  // 获取组内进程属性，获取组内第一个进程的id
		tmp_int=0;
		for(j=0;j<i;j++){
			tmp_int+=com->size.com[j];
		}
		com->master.com[i]=tmp_int;
	}


MPI_Comm_split(MPI_COMM_WORLD,com->color_group,com->rank_world,&com->com_group);  // 一个通信组内切分进程
	MPI_Comm_rank(com->com_group,&com->rank_group);
	MPI_Comm_size(com->com_group,&tmp_int);	// possibible to collect the size of com_group
	

创建进程组之间的通信组

int mpi_lsa_create_intercoms(com_lsa * com){ ///创建进程组之间的通信组
		com->master.com[i]=tmp_int;

//第一种方式：调用系统函数创建
MPI_Intercomm_create(com->com_group,0,
					MPI_COMM_WORLD,com->master.com[4-((com->color_group)+1)],
					com->rank_group,
					&(com->inter.com[4-((com->color_group)+1)]));
	}
//第二种方式
MPI_Intercomm_create(com->com_group,0,
					com->com_world,com->master.com[(4-((com->color_group)-1)%4)%4],
	 				com->rank_group,
	 				&(com->inter.com[(4-((com->color_group)-1)%4)%4]));

//设置输入输出通信组
com->out_com=com->inter.com[next];
	MPI_Comm_test_inter(com->inter.com[next],&flag);
		if(!flag){
			mpi_lsa_print("\n\n\n\nproblem with inter.[next]\n\n\n\n\n", com);
		}
	com->in_com=com->inter.com[prev];
	MPI_Comm_test_inter(com->inter.com[prev],&flag);
	if(!flag){
		mpi_lsa_print("\n\n\n\n\nproblem with inter.[prev]\n\n\n\n", com);
	}

//同步之后，根据进程属性对应的值输出执行的算法
MPI_Barrier(MPI_COMM_WORLD);

	if(com->color_group==0) 		 printf("GMRES :   ");
	else if(com->color_group==1) printf("MAIN :    ");
	else if(com->color_group==2) printf("ARNOLDI : ");
	else if(com->color_group==3) printf("LS :      ");
	
	

int mpi_lsa_print(char * s,com_lsa * com){  //  0号进程换行

2、else{}

执行经典的lsa算法

read_matrix_vector(&A, &v,PETSC_COMM_WORLD);//把文件A，文件v数据读进计算机的变量A(是一个稀疏矩阵形式)和v
//此函数内部：
// A
ierr=PetscOptionsGetString(NULL,PETSC_NULL,"-mfile",filea,PETSC_MAX_PATH_LEN-1,&flaga);CHKERRQ(ierr); // 根据命令行，把参数传给filea
PetscPrintf(PETSC_COMM_WORLD,"Loading Matrix : %s\n",filea);// 
	ierr=PetscViewerBinaryOpen(PETSC_COMM_WORLD,filea,FILE_MODE_READ,&fd);CHKERRQ(ierr);  // 从filea中获取文件路径
	ierr = MatCreate(PETSC_COMM_WORLD,A);CHKERRQ(ierr);  // 创建A
ierr=MatLoad(*A,fd);CHKERRQ(ierr); // 把fd文件加载到A
ierr=MatGetSize(*A,&size,&sizea);CHKERRQ(ierr);  // 获取A的行列信息
 
//b
ierr=PetscOptionsGetString(NULL,PETSC_NULL,"-vfile",fileb,PETSC_MAX_PATH_LEN-1,&flagb);CHKERRQ(ierr);  // -vfile"命令行传入b

if (!flagb) {  // 随机产生一个b向量
		/* the user did not provide a vector, so generate it*/
		generate_random_seed_vector(size, -10.0, 10.0, 0, v);
		VecNorm(*v, NORM_2, &vnorm);  // 2范数
		VecScale(*v, 0.01/vnorm);	  // 单位化	
		PetscPrintf(PETSC_COMM_WORLD,"Generated right hand side matrix b\n");		
	} else {  直接加载b传给v指向的内存
		PetscPrintf(PETSC_COMM_WORLD,"Loading Vector : %s\n",fileb);
		ierr=PetscViewerBinaryOpen(PETSC_COMM_WORLD,fileb,FILE_MODE_READ,&fd);CHKERRQ(ierr);
		VecCreate(PETSC_COMM_WORLD,v);
		ierr=VecLoad(*v,fd);CHKERRQ(ierr);
		ierr=PetscViewerDestroy(&fd);CHKERRQ(ierr);
		PetscPrintf(PETSC_COMM_WORLD,"Loaded Vector of size : %d\n",size);
		

3、if()

执行自己写的算法

PETSC_COMM_WORLD=com.com_group;  // 初始化一个通信组
PetscPrintf(com.com_world,"]> Initializing PETSc/SLEPc\n");
ierr=PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the GMRES/LS-Arnoldi solver", "KSP");CHKERRQ(ierr);// 输出提示信息
MPI_Barrier(MPI_COMM_WORLD);  // mpi同步

if(com.color_group==GMRES_GR || (com.color_group==ARNOLDI_GR && nols!=0)){  // 按照分好的GMRES_GR组进行读取
        	read_matrix_vector(&A, &v,(MPI_Comm*) &com.color_group);  //读取数据文件
		MatGetSize(A,NULL, &matsize);  // 获取数据
MPI_Bcast(&matsize, 1, MPIU_INT,0, MPI_COMM_WORLD);  // 广播到所有进程


if(com.color_group==GMRES_GR || (com.color_group==ARNOLDI_GR && nols!=0)){  // 对A操作
		Af=A;
		vf=v;

		PetscOptionsName("-Aclx","Convert complex matrix to real","PetscAclx",&flag);// // 解析命令行，赋值给flag
		if (flag) {  
			ierr= real2complexMat(&A,&Af);CHKERRQ(ierr);  //矩阵变形转成复数
			ierr= real2complexVec(&v, &vf);CHKERRQ(ierr);
		}

		A=Af;  // A阵用变形过的
		vf=v;  // 变过的vf还用原来的v

else if((com.color_group==FATHER_GR && nols!=0) || (com.color_group==LS_GR && nols!=0)){  // 预处理 解析命令行，初始化b向量
		VecCreate(PETSC_COMM_WORLD,&v);  // 创建v
		VecSetSizes(v,PETSC_DECIDE,matsize);  // 设置v的大小
		VecSetFromOptions(v);
	}
        }


mpi_lsa_com_vecsize_init(&com,&v);  // 初始化lsa
ierr=PetscOptionsGetInt(NULL,PETSC_NULL,"-ksp_ls_nopc",&nols,&flag);CHKERRQ(ierr);  // 解析命令行，初始化flag
if(!flag) nols=1;  // flag==0 时预处理
	else if(nols==0) PetscPrintf(com.com_world,"]> not using LSA preconditionning %d\n",nols);

PetscOptionsHasName(NULL,NULL,"-GMRES_FT",&gft_flg);  // 解析命令行，赋值到变量

if(com.color_group==GMRES_GR){
			PetscPrintf(com.com_group,"]> Launching GMRES\n");
			launchGMRES(&com, &v, &A);  // 启动GMRES
		}else if(com.color_group==ARNOLDI_GR && nols!=0){
			PetscPrintf(com.com_group,"]> Launching Arnoldi\n");
			Arnoldi(&com,&A,&v)	;
			}
		else if(com.color_group==FATHER_GR && nols!=0){  // 预处理的组
			PetscPrintf(com.com_group,"]> Launching Father\n");
			Father(&com,&v);
		} else if(com.color_group==LS_GR && nols!=0){  // 预处理的组
			VecGetSize(v,&vsize);
			PetscPrintf(com.com_group,"]> Launching LS\n");
			LSQR(&com,&vsize);