《并行算法在电阻抗断层成像中的应用研究》15000字（论文）

并行算法在电阻抗断层成像中的应用研究目录TOC\o"1-3"\f\h\z\u5897400引言 、简要介绍了电阻抗层析成像技术的背景，国内外研究的重要性以及当前的研究现状，并简要介绍了当前的EIT算法.对EIT场域进行了数学描述并推导了边界微分方程，并且描述了EIT正问题常规的求解方法，有限元方法(FEM)、有限体元法(FVM)、边界元法(BEM)，并着重介绍了有限元法，推导了其求解正问题的方程的过程.同时也介绍了EIT的逆问题.逆问题是一个高度病态的非线性问题，直接求解难度较大，本文介绍了一些常见的求解算法，并着重介绍了牛顿-拉夫逊算法，但是由于推导过程中对高阶导变量的舍弃，使得最终解会影响最终成像质量，于是就介绍了正则化技术对牛顿-拉夫逊算法的修正，着重介绍并推导了Tikhonov正则化算法对其的优化.在有限元法求解EIT的基础上，由于其程序的核心算法是对矩阵迭代计算，传统的算法模型时间长，冗余度高，于是提出了热门的并行计算模型，并且介绍了基于并行计算，以及一些先进的数学问题求解的编程库以及工具，如DUMEC++库，PETSC工具等实现的并行机器和集群上运行的求解器(PEIT)，PEIT在CPU上的执行效率已经大大优于了传统EIT求解程序，然而这个求解器虽然已经优化了大部分的计算模型，但是Jacobian矩阵的计算还是基于预处理，计算速度不是很理想.于是着重介绍了NVIDIA公司的GPU架构体系，对最新架构进行了介绍和分析，并介绍以及搭建了基于GPU的CUDA编程模型.对基于DUNE-FEM，PETSC的Jacobian矩阵求解算法，使用基于GPU的架构，以及多个GPU架构的对Jacobian矩阵求解进行加速.

参考文献王雅娟.腔内电阻抗成像正问题边界元法求解的并行计算研究[D].天津:河北工业大学,2014.代月霞.基于深度学习的EIT图像重建算法研究[D].天津:天津工业大学,2017.丁红军,邢克礼.医学成像技术的进展[J].医疗卫生装备,2006(11):22-23.常晓敏.提高电阻抗重建图像质量的方法研究[D].天津:天津科技大学,2018.向胜昭.基于三维磁场聚焦技术的磁感应成像系统中激励线圈的设计[D].成都:四川大学,2005.张夏婉.多频电阻抗成像算法研究[D].南京:南京邮电大学,2018.唐磊.电阻抗断层成像算法研究及系统软件设计[D].天津:天津大学,2006.田明武.带裂缝水泥浆体阻抗特性研究[D].成都:四川大学,2006.周舟.差分电阻抗断层成像算法研究[D].长沙:国防科学技术大学,2015.杨尚琴.多层次并行算法与MPI-2新特性的研究及应用[D].成都:成都理工大学,2009.李克臣.波动方程的速度估计和叠前深度偏移辛算法的并行实现[D].合肥:中国科学技术大学,2002.王宏斌.模拟胸腔断层的二维椭圆域阻抗成像研究[D].天津:河北工业大学,2007.付峰,秦明新.电阻抗断层图像重构算法[J].国外医学.生物医学工程分册,1995(05):255-262.曹海燕,王化祥.电容层析成像系统传感器仿真研究[J].计算机工程与应用,2012,48(15):207-211.徐桂芝,王明时,李有余,等.医用电阻抗成像系统的模块化设计[J].天津:天津大学学报,2006(06):133-137.任燕芝.改进的粒子群算法及其工程应用研究[D].西安:西安电子科技大学,2018.戎立锋.电阻抗成像技术的研究与系统设计[D].南京:南京理工大学,2006.张建国.基于阻抗层析成像的冰—水两相流检测技术与系统研究[D].太原:太原理工大学,2013.邵兴龙.基于FPGA的图像轮廓提取系统设计与实现[D].无锡:江南大学,2014.苏丽丽.基于CPU-GPU集群的分子动力学并行计算研究[D].大连:大连理工大学,2009.王贞东.增强现实中虚实融合光照一致性研究[D].苏州:苏州大学,2010.赵芳斌.基于FPGA和GPU的外辐射源雷达信号处理平台的研究[D].成都:电子科技大学,2013.陆开诚.LAPW基组第一性原理计算的GPU加速方法及其应用[D].杭州:浙江大学,2019.高跃清,张焱,刘伟光.基于CUDA的SAR成像CS算法研究[J].计算机与网络,2012,38(07):55-57.JehlMarkus,DednerAndreas,BetckeTimo,etal.Afastparallelsolverfortheforwardprobleminelectricalimpedancetomography[J].IEEEtransactionsonbio-medicalengineering,2015,62(1):1-14.吴玫华.在GPU上实现Jacobi迭代法的分析与设计[J].电子设计工程,2012,20(10):28-30.董秀珍,秦明新,汤孟兴,等.电阻抗断层成像系统及重构算法[J].第四军医大学学报,1999(03):218-219.ABorsic,EAAttardo,RJHalter.Multi-GPUJacobianacceleratedcomputingforsoft-fieldtomography[J].PhysiologicalMeasurement,2012,33(10):1703-1715.

附录intmain(intargc,char**argv)

{if(argc>1){printf("Thisprogramacceptsnoarguments\n");exit(EXIT_FAILURE);}matrix_t

A;

matrix_t

B;

matrix_treference_x;

matrix_tgpu_naive_solution_x;

matrix_tgpu_opt_solution_x;

srand(time(NULL));

printf("\nGenerating%dx%dsystem\n",MATRIX_SIZE,MATRIX_SIZE);A=create_diagonally_dominant_matrix(MATRIX_SIZE,MATRIX_SIZE);if(A.elements==NULL){

printf("Errorcreatingmatrix\n");

exit(EXIT_FAILURE);}

B=allocate_matrix_on_host(MATRIX_SIZE,1,1);reference_x=allocate_matrix_on_host(MATRIX_SIZE,1,0);gpu_naive_solution_x=allocate_matrix_on_host(MATRIX_SIZE,1,0);

gpu_opt_solution_x=allocate_matrix_on_host(MATRIX_SIZE,1,0);#ifdefDEBUGprint_matrix(A);print_matrix(B);print_matrix(reference_x);#endif

printf("\nPerformingJacobiiterationontheCPU\n");

structtimevalstart,stop;

gettimeofday(&start,NULL);

compute_gold(A,reference_x,B);

gettimeofday(&stop,NULL);

fprintf(stderr,"CPUruntime=%fs\n",(float)(stop.tv_sec-start.tv_sec+(stop.tv_usec-start.tv_usec)/(float)1000000));

display_jacobi_solution(A,reference_x,B);/*Displaystatistics*/

printf("\nPerformingJacobiiterationondevice\n");

compute_on_device(A,gpu_naive_solution_x,gpu_opt_solution_x,B);

display_jacobi_solution(A,gpu_naive_solution_x,B);/*Displaystatistics*/

display_jacobi_solution(A,gpu_opt_solution_x,B);

free(A.elements);

free(B.elements);

free(reference_x.elements);

free(gpu_naive_solution_x.elements);

free(gpu_opt_solution_x.elements);

exit(EXIT_SUCCESS);}voidcompute_on_device(constmatrix_tA,matrix_tgpu_naive_sol_x,

matrix_tgpu_opt_sol_x,constmatrix_tB){

structtimevalstart,stop;

matrix_tAd=allocate_matrix_on_device(A);

copy_matrix_to_device(Ad,A);

matrix_tBd=allocate_matrix_on_device(B);

copy_matrix_to_device(Bd,B);

matrix_tXd=allocate_matrix_on_device(B);

copy_matrix_to_device(Xd,B);

gettimeofday(&start,NULL);

jacobi_iteration_naive(Ad,Bd,Xd);

copy_matrix_from_device(gpu_naive_sol_x,Xd);

gettimeofday(&stop,NULL);

fprintf(stderr,"GPUruntimefornaive=%fs\n",(float)(stop.tv_sec-start.tv_sec+(stop.tv_usec-start.tv_usec)/(float)1000000));

copy_matrix_to_device(Xd,B);

constmatrix_tnewA=transpose_matrix(A);

copy_matrix_to_device(Ad,newA);

gettimeofday(&start,NULL);

jacobi_iteration_optimized(Ad,Bd,Xd);

copy_matrix_from_device(gpu_opt_sol_x,Xd);

gettimeofday(&stop,NULL);

fprintf(stderr,"CPUruntimeforoptimized=%fs\n",(float)(stop.tv_sec-start.tv_sec+(stop.tv_usec-start.tv_usec)/(float)1000000));

free((void*)newA.elements);

free_matrix_on_device(&Ad);

free_matrix_on_device(&Bd);

free_matrix_on_device(&Xd);

return;}voidjacobi_iteration_naive(matrix_tAd,matrix_tBd,matrix_tXd){

dim3threads(THREAD_BLOCK_SIZE,1,1);

dim3grid((Xd.num_rows+THREAD_BLOCK_SIZE-1)/THREAD_BLOCK_SIZE,1);

double*ssdDevice;

cudaMalloc((void**)&ssdDevice,sizeof(double));

unsignedintdone=0;

doublessd,mse;

unsignedintnum_iter=0;

while(!done){

cudaMemset((void*)ssdDevice,0,sizeof(double));

jacobi_iteration_kernel_naive<<<grid,threads>>>(Ad,Bd,Xd,ssdDevice);

cudaDeviceSynchronize();/*ForceCPUtowaitforGPUtocomplete*/

cudaMemcpy((void*)&ssd,(void*)ssdDevice,sizeof(double),cudaMemcpyDeviceToHost);

num_iter++;

mse=sqrt(ssd);/*Meansquarederror.*/

if(mse<=THRESHOLD)

done=1;

}

printf("\nConvergenceachievedafter%diterations\n",num_iter);

cudaFree(ssdDevice);}voidjacobi_iteration_optimized(matrix_tAd,matrix_tBd,matrix_tXd){

dim3threads(THREAD_BLOCK_SIZE,1,1);

dim3grid((Xd.num_rows+THREAD_BLOCK_SIZE-1)/THREAD_BLOCK_SIZE,1);

double*ssdDevice;

cudaMalloc((void**)&ssdDevice,sizeof(double));

unsignedintdone=0;

doublessd,mse;

unsignedintnum_iter=0;

while(!done){

cudaMemset((void*)ssdDevice,0,sizeof(double));

jacobi_iteration_kernel_optimized<<<grid,threads>>>(Ad,Bd,Xd,ssdDevice);

cudaDeviceSynchronize();/*ForceCPUtowaitforGPUtocomplete*/

cudaMemcpy((void*)&ssd,(void*)ssdDevice,sizeof(double),cudaMemcpyDeviceToHost);

num_iter++;

mse=sqrt(ssd);/*Meansquarederror.*/

if(mse<=THRESHOLD)

done=1;

}

printf("\nConvergenceachievedafter%diterations\n",num_iter);

cudaFree(ssdDevice);}voidfree_matrix_on_device(matrix_t

*M)

{cudaFree(M->elements);M->elements=NULL;}voidfree_matrix_on_host(matrix_t*M){free(M->elements);M->elements=NULL;}matrix_tallocate_matrix_on_device(constmatrix_tM){

matrix_tMdevice=M;

intsize=M.num_rows*M.num_columns*sizeof(float);

cudaMalloc((void**)&Mdevice.elements,size);

returnMdevice;}matrix_ttranspose_matrix(constmatrix_tM){

matrix_tnewM;

newM.num_columns=M.num_columns;

newM.num_rows=M.num_rows;

intsize=M.num_rows*M.num_columns;

newM.elements=(float*)malloc(size*sizeof(float));

introw,cols;

for(inti=0;i<size;i++){

row=i/M.num_rows;

cols=i%M.num_rows;

newM.elements[cols*M.num_rows+row]=M.elements[i];

}

returnnewM;}matrix_tallocate_matrix_on_host(intnum_rows,intnum_columns,intinit){

matrix_tM;

M.num_columns=num_columns;

M.num_rows=num_rows;

intsize=M.num_rows*M.num_columns;M.elements=(float*)malloc(size*sizeof(float));for(unsignedinti=0;i<size;i++){if(init==0)

M.elements[i]=0;

else

M.elements[i]=get_random_number(MIN_NUMBER,MAX_NUMBER);}

returnM;}voidcopy_matrix_to_device(matrix_tMdevice,constmatrix_tMhost){

intsize=Mhost.num_rows*Mhost.num_columns*sizeof(float);

Mdevice.num_rows=Mhost.num_rows;

Mdevice.num_columns=Mhost.num_columns;

cudaMemcpy(Mdevice.elements,Mhost.elements,size,cudaMemcpyHostToDevice);

return;}voidcopy_matrix_from_device(matrix_tMhost,constmatrix_tMdevice){

intsize=Mdevice.num_rows*Mdevice.num_columns*sizeof(float);

cudaMemcpy(Mhost.elements,Mdevice.elements,size,cudaMemcpyDeviceToHost);

return;}voidprint_matrix(constmatrix_tM){for(unsignedinti=0;i<M.num_rows;i++){

for(unsignedintj=0;j<M.num_columns;j++){printf("%f",M.elements[i*M.num_rows+j]);

}