《并行算法在電阻抗斷層成像中的應(yīng)用研究》15000字（論文）

并行算法在電阻抗斷層成像中的應(yīng)用研究目錄TOC\o"1-3"\f\h\z\u5897400引言 、簡要介紹了電阻抗層析成像技術(shù)的背景，國內(nèi)外研究的重要性以及當(dāng)前的研究現(xiàn)狀，并簡要介紹了當(dāng)前的EIT算法.對EIT場域進(jìn)行了數(shù)學(xué)描述并推導(dǎo)了邊界微分方程，并且描述了EIT正問題常規(guī)的求解方法，有限元方法(FEM)、有限體元法(FVM)、邊界元法(BEM)，并著重介紹了有限元法，推導(dǎo)了其求解正問題的方程的過程.同時也介紹了EIT的逆問題.逆問題是一個高度病態(tài)的非線性問題，直接求解難度較大，本文介紹了一些常見的求解算法，并著重介紹了牛頓-拉夫遜算法，但是由于推導(dǎo)過程中對高階導(dǎo)變量的舍棄，使得最終解會影響最終成像質(zhì)量，于是就介紹了正則化技術(shù)對牛頓-拉夫遜算法的修正，著重介紹并推導(dǎo)了Tikhonov正則化算法對其的優(yōu)化.在有限元法求解EIT的基礎(chǔ)上，由于其程序的核心算法是對矩陣迭代計算，傳統(tǒng)的算法模型時間長，冗余度高，于是提出了熱門的并行計算模型，并且介紹了基于并行計算，以及一些先進(jìn)的數(shù)學(xué)問題求解的編程庫以及工具，如DUMEC++庫，PETSC工具等實現(xiàn)的并行機器和集群上運行的求解器(PEIT)，PEIT在CPU上的執(zhí)行效率已經(jīng)大大優(yōu)于了傳統(tǒng)EIT求解程序，然而這個求解器雖然已經(jīng)優(yōu)化了大部分的計算模型，但是Jacobian矩陣的計算還是基于預(yù)處理，計算速度不是很理想.于是著重介紹了NVIDIA公司的GPU架構(gòu)體系，對最新架構(gòu)進(jìn)行了介紹和分析，并介紹以及搭建了基于GPU的CUDA編程模型.對基于DUNE-FEM，PETSC的Jacobian矩陣求解算法，使用基于GPU的架構(gòu)，以及多個GPU架構(gòu)的對Jacobian矩陣求解進(jìn)行加速.

參考文獻(xiàn)王雅娟.腔內(nèi)電阻抗成像正問題邊界元法求解的并行計算研究[D].天津:河北工業(yè)大學(xué),2014.代月霞.基于深度學(xué)習(xí)的EIT圖像重建算法研究[D].天津:天津工業(yè)大學(xué),2017.丁紅軍,邢克禮.醫(yī)學(xué)成像技術(shù)的進(jìn)展[J].醫(yī)療衛(wèi)生裝備,2006(11):22-23.常曉敏.提高電阻抗重建圖像質(zhì)量的方法研究[D].天津:天津科技大學(xué),2018.向勝昭.基于三維磁場聚焦技術(shù)的磁感應(yīng)成像系統(tǒng)中激勵線圈的設(shè)計[D].成都:四川大學(xué),2005.張夏婉.多頻電阻抗成像算法研究[D].南京:南京郵電大學(xué),2018.唐磊.電阻抗斷層成像算法研究及系統(tǒng)軟件設(shè)計[D].天津:天津大學(xué),2006.田明武.帶裂縫水泥漿體阻抗特性研究[D].成都:四川大學(xué),2006.周舟.差分電阻抗斷層成像算法研究[D].長沙:國防科學(xué)技術(shù)大學(xué),2015.楊尚琴.多層次并行算法與MPI-2新特性的研究及應(yīng)用[D].成都:成都理工大學(xué),2009.李克臣.波動方程的速度估計和疊前深度偏移辛算法的并行實現(xiàn)[D].合肥:中國科學(xué)技術(shù)大學(xué),2002.王宏斌.模擬胸腔斷層的二維橢圓域阻抗成像研究[D].天津:河北工業(yè)大學(xué),2007.付峰,秦明新.電阻抗斷層圖像重構(gòu)算法[J].國外醫(yī)學(xué).生物醫(yī)學(xué)工程分冊,1995(05):255-262.曹海燕,王化祥.電容層析成像系統(tǒng)傳感器仿真研究[J].計算機工程與應(yīng)用,2012,48(15):207-211.徐桂芝,王明時,李有余,等.醫(yī)用電阻抗成像系統(tǒng)的模塊化設(shè)計[J].天津:天津大學(xué)學(xué)報,2006(06):133-137.任燕芝.改進(jìn)的粒子群算法及其工程應(yīng)用研究[D].西安:西安電子科技大學(xué),2018.戎立鋒.電阻抗成像技術(shù)的研究與系統(tǒng)設(shè)計[D].南京:南京理工大學(xué),2006.張建國.基于阻抗層析成像的冰—水兩相流檢測技術(shù)與系統(tǒng)研究[D].太原:太原理工大學(xué),2013.邵興龍.基于FPGA的圖像輪廓提取系統(tǒng)設(shè)計與實現(xiàn)[D].無錫:江南大學(xué),2014.蘇麗麗.基于CPU-GPU集群的分子動力學(xué)并行計算研究[D].大連:大連理工大學(xué),2009.王貞東.增強現(xiàn)實中虛實融合光照一致性研究[D].蘇州:蘇州大學(xué),2010.趙芳斌.基于FPGA和GPU的外輻射源雷達(dá)信號處理平臺的研究[D].成都:電子科技大學(xué),2013.陸開誠.LAPW基組第一性原理計算的GPU加速方法及其應(yīng)用[D].杭州:浙江大學(xué),2019.高躍清,張焱,劉偉光.基于CUDA的SAR成像CS算法研究[J].計算機與網(wǎng)絡(luò),2012,38(07):55-57.JehlMarkus,DednerAndreas,BetckeTimo,etal.Afastparallelsolverfortheforwardprobleminelectricalimpedancetomography[J].IEEEtransactionsonbio-medicalengineering,2015,62(1):1-14.吳玫華.在GPU上實現(xiàn)Jacobi迭代法的分析與設(shè)計[J].電子設(shè)計工程,2012,20(10):28-30.董秀珍,秦明新,湯孟興,等.電阻抗斷層成像系統(tǒng)及重構(gòu)算法[J].第四軍醫(yī)大學(xué)學(xué)報,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]);

}