00001 #include "boundparcongrad.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <time.h>
00005 #include <vector.h>
00006 #include "petsc_boundcg.h"
00007
00008
00009
00010
00011
00012 boundparcongrad::boundparcongrad(int np,int mr,long nd,long mes)
00013 {
00014 mespr = mes;
00015
00016 nproc=np;
00017
00018 myrank=mr;
00019
00020 ndom=nd;
00021
00022 ndof=0;
00023
00024 nbdof=0;
00025
00026 nidof = 0;
00027
00028 maxnbdof=0;
00029
00030
00031 ni = 0;
00032
00033 err = 0.0;
00034
00035 prec = -1;
00036
00037 lidof= NULL;
00038
00039 lbdof = NULL;
00040
00041 bcn = NULL;
00042
00043 nbdofdom = NULL;
00044
00045 petscpoint = new PETSC_CONT;
00046
00047 }
00048
00049
00050
00051
00052
00053
00054 boundparcongrad::~boundparcongrad()
00055 {
00056 if(prec == petscilu){
00057 MatDestroy(petscpoint->fact);
00058 VecDestroy(petscpoint->r_rhs);
00059 VecDestroy(petscpoint->z_lhs);
00060
00061 PetscFinalize();
00062
00063 delete []petscpoint->val;
00064 delete []petscpoint->ix;
00065 }
00066 delete petscpoint;
00067
00068 }
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 void boundparcongrad::initiate (selnodes *selnodschur,partop *ptop,gtopology *top,FILE *out)
00082 {
00083 long i,j;
00084
00085
00086
00087 if(prec == petscilu){
00088 static char help[] = " ";
00089 PetscInitialize(&Argc,&Argv,(char *)0,help);
00090 }
00091
00092 if (myrank==0){
00093
00094 nbdofdom = selnodschur -> snndofmas;
00095
00096 bcn = selnodschur -> cnmas;
00097
00098 ndofcp = selnodschur -> tndofsn;
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 }
00109 maxnbdof = selnodschur->maxsnndof;
00110
00111 ndof=ptop->ndof;
00112 nidof=ptop->nidof;
00113 nbdof=ptop->nbdof;
00114
00115 }
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127 double boundparcongrad::v_ixu_i(double *u,double *v)
00128 {
00129 long i;
00130 double sp;
00131 sp = 0.0;
00132 for(i = 0; i < nidof; i++){
00133
00134
00135 sp+=u[i]*v[i];
00136 }
00137
00138 return sp;
00139 }
00140
00141
00142
00143
00144
00145
00146
00147
00148 void boundparcongrad::select_bound(double *buff,double *v)
00149 {
00150 long i;
00151
00152 for(i = nidof; i < ndof; i++){
00153 buff[i-nidof] = v[i];
00154 }
00155 }
00156
00157
00158
00159
00160
00161
00162
00163 void boundparcongrad::add_bound(double *buff,double *v)
00164 {
00165 long i;
00166 for(i = 0; i < nbdof; i++){
00167 v[i+nidof] = buff[i];
00168 }
00169 }
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180 void boundparcongrad::add_master_buff (double *buff,double *cv,long nd)
00181 {
00182 long i,j;
00183 for (i=0;i<nbdofdom[nd];i++){
00184
00185 j=bcn[nd][i]-1;
00186 cv[j]+=buff[i];
00187 }
00188
00189 }
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 void boundparcongrad::select_master_buff(double *cv,double *buff,long nd)
00200 {
00201 long i,j;
00202 for (i=0;i<nbdofdom[nd];i++){
00203
00204 j=bcn[nd][i]-1;
00205 buff[i] = cv[j];
00206 }
00207
00208 }
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219 void boundparcongrad::ilu_mat_petsc(gmatrix *gm,FILE *out)
00220 {
00221
00222 long adr1,adr2;
00223
00224
00225 PetscInt i,j,jj,nz,m;
00226 PetscInt *nnz;
00227 IS row,col;
00228
00229
00230
00231 PetscScalar v;
00232
00233
00234
00235
00236
00237 switch(gm->ts){
00238 case compressed_rows:{
00239
00240
00241 nz = 0;
00242 nnz = new PetscInt[gm->cr->n];
00243
00244
00245
00246 for(i = 0; i < gm->cr->n; i++){
00247
00248 adr1=gm->cr->adr[i];
00249 adr2=gm->cr->adr[i+1];
00250 nnz[i]=adr2-adr1;
00251 }
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261 MatCreateSeqAIJ(PETSC_COMM_SELF,gm->cr->n,gm->cr->n,nz,nnz,&petscpoint->fact);
00262
00263
00264 for(i = 0; i < gm->cr->n; i++){
00265 adr1=gm->cr->adr[i];
00266 adr2=gm->cr->adr[i+1];
00267 for(j = adr1; j < adr2; j++){
00268 v = gm->cr->a[j];
00269 jj=gm->cr->ci[j];
00270
00271 MatSetValues(petscpoint->fact,1,&i,1,&jj,&v,INSERT_VALUES);
00272 }
00273 }
00274 ISCreateStride(PETSC_COMM_WORLD,gm->cr->n,0,1,&row);
00275 ISCreateStride(PETSC_COMM_WORLD,ndof,0,1,&col);
00276 break;
00277 }
00278 case symm_comp_rows:{
00279
00280 nz = 0;
00281 nnz = new PetscInt[gm->scr->n];
00282 for(i = 0; i < gm->scr->n; i++){
00283 nnz[i] = 0;
00284 }
00285
00286
00287
00288
00289 for(i = 0; i < gm->scr->n; i++){
00290
00291 adr1=gm->scr->adr[i];
00292 adr2=gm->scr->adr[i+1];
00293 nnz[i] = adr2-adr1;
00294 for(j = i+1; j < gm->scr->n; j++){
00295 for(m = adr1=gm->scr->adr[j]; m < gm->scr->adr[j+1]; m++){
00296 if(gm->scr->ci[m] == i){
00297 nnz[i]++;
00298 }
00299 }
00300 }
00301 }
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313 MatCreateSeqAIJ(PETSC_COMM_SELF,gm->scr->n,gm->scr->n,nz,nnz,&petscpoint->fact);
00314
00315
00316 for(i = 0; i < gm->scr->n; i++){
00317 adr1=gm->scr->adr[i];
00318 adr2=gm->scr->adr[i+1];
00319 for(j = adr1; j < adr2; j++){
00320 v = gm->scr->a[j];
00321 jj=gm->scr->ci[j];
00322
00323 MatSetValues(petscpoint->fact,1,&i,1,&jj,&v,INSERT_VALUES);
00324 MatSetValues(petscpoint->fact,1,&jj,1,&i,&v,INSERT_VALUES);
00325 }
00326 }
00327 ISCreateStride(PETSC_COMM_WORLD,gm->scr->n,0,1,&row);
00328 ISCreateStride(PETSC_COMM_WORLD,ndof,0,1,&col);
00329 break;
00330 }
00331 default:{
00332
00333 par_print_err(myrank,"Unknown storage type\n", __FILE__, __LINE__, __func__);
00334 break;
00335 }
00336 }
00337
00338
00339 MatAssemblyBegin(petscpoint->fact,MAT_FINAL_ASSEMBLY);
00340 MatAssemblyEnd(petscpoint->fact,MAT_FINAL_ASSEMBLY);
00341
00342
00343
00344 MatFactorInfo info;
00345 info.fill = 1.0;
00346
00347
00348
00349
00350 info.levels=0;
00351
00352
00353 MatILUFactor(petscpoint->fact,row,col,&info);
00354
00355
00356 i = ndof;
00357 VecCreateSeq(PETSC_COMM_SELF,i,&petscpoint->r_rhs);
00358 VecAssemblyBegin(petscpoint->r_rhs);
00359 VecAssemblyEnd(petscpoint->r_rhs);
00360 i = ndof;
00361 VecCreateSeq(PETSC_COMM_SELF,i,&petscpoint->z_lhs);
00362 VecAssemblyBegin(petscpoint->r_rhs);
00363 VecAssemblyEnd(petscpoint->r_rhs);
00364
00365
00366 petscpoint->ix = new PetscInt[ndof];
00367 petscpoint->val = new PetscScalar[ndof];
00368
00369 delete []nnz;
00370 }
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382 void boundparcongrad::solve_fact_mat_petsc(double *r,double *z,long ndof,FILE *out)
00383 {
00384
00385 PetscInt i,ngv;
00386 PetscScalar v;
00387
00388
00389
00390
00391
00392
00393
00394
00395 for(i = 0; i < ndof; i++){
00396 v = r[i];
00397 VecSetValue(petscpoint->r_rhs,i,v,INSERT_VALUES);
00398 }
00399
00400
00401 MatSolve(petscpoint->fact,petscpoint->r_rhs,petscpoint->z_lhs);
00402
00403
00404 nullv(z,ndof);
00405 i = ndof;
00406
00407
00408 for(i = 0; i < ndof; i++){
00409 petscpoint->ix[i] = i;
00410 }
00411 ngv = ndof;
00412 VecGetValues(petscpoint->z_lhs,ngv,petscpoint->ix,petscpoint->val);
00413
00414 for(i = 0; i < ndof; i++){
00415
00416 z[petscpoint->ix[i]] =petscpoint->val[petscpoint->ix[i]];
00417 }
00418
00419
00420 }
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 void boundparcongrad::initialize_diag_prec(gmatrix *gm,double *precvec,FILE *out)
00432 {
00433 long i,j,jj,adr1,adr2;
00434
00435 switch(gm->ts){
00436 case compressed_rows:{
00437 for(i = 0; i < gm->cr->n; i++){
00438 adr1=gm->cr->adr[i];
00439 adr2=gm->cr->adr[i+1];
00440 for(j = adr1; j < adr2; j++){
00441 jj=gm->cr->ci[j];
00442 if(jj == i){
00443 precvec[i] = gm->cr->a[j];
00444 precvec[i]=1.00/precvec[i];
00445 break;
00446 }
00447
00448 }
00449
00450 }
00451 break;
00452 }
00453 case symm_comp_rows:{
00454 for(i = 0; i < gm->scr->n; i++){
00455 adr1=gm->scr->adr[i];
00456 adr2=gm->scr->adr[i+1];
00457 for(j = adr1; j < adr2; j++){
00458 jj=gm->scr->ci[j];
00459 if(jj == i){
00460 precvec[i] = gm->scr->a[j];
00461 precvec[i] = 1.00/precvec[i];
00462 }
00463 }
00464
00465 if(mespr == 1) fprintf(out,"prec[%ld]= %le\n",i,precvec[i]);
00466
00467 }
00468 break;
00469 }
00470 default:{
00471 printf("Unknown storage type of matrix");
00472 break;
00473 }
00474 }
00475
00476 }
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488 void boundparcongrad::jacobi_precondition(double *precvec,double *invec,double *outvec,long ndof,FILE *out)
00489 {
00490 long i;
00491 for(i = 0; i < ndof; i++){
00492 outvec[i] = precvec[i]*invec[i];
00493 if(mespr == 1) fprintf(out,"precvec[%ld]= %le incvec[%ld] = %le outvec[%ld] = %le\n",i,precvec[i],i,invec[i],i,outvec[i]);
00494 }
00495 }
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517 void boundparcongrad::solve_system (partop *ptop,gtopology *top,gmatrix *gm,long *domproc,double *lhs,double *rhs,FILE *out,long iv)
00518 {
00519 long i,j,k,iter;
00520 long buffsize;
00521 double nom_i,denom_i,alpha,beta,norrhs;
00522 double *p_dom, *r_dom, *d_dom,*x_dom,*buff,*z_dom,*aux;
00523 double *precvec;
00524 double tbacksub;
00525 double nom_master_i,nom_master_b,nom,denom;
00526 double *r_bound, *d_bound, *x_bound, *p_bound,*z_bound;
00527 double denom_master_i,denom_master_b;
00528 double tfacts,tfacte,tsols,tsole,titers,titere,ts,te;
00529 double nory;
00530
00531 MPI_Status stat;
00532
00533 tbacksub=0.0;
00534 tsols = clock();
00535 switch(prec){
00536 case noprec:{
00537 break;
00538 }
00539 case petscilu:{
00540 tfacts = clock();
00541 ilu_mat_petsc(gm,out);
00542 tfacte = clock();
00543 if(mespr == 1) fprintf (out,"Time of ILU factorization - ilu_mat_petsc function %le s\n",(tfacte-tfacts)/(double)CLOCKS_PER_SEC);
00544 break;
00545 }
00546 case pardiagprec:{
00547 precvec = new double[ndof];
00548 initialize_diag_prec(gm,precvec,out);
00549 break;
00550 }
00551 }
00552
00553
00554
00555
00556
00557 p_dom = new double [ndof];
00558 nullv (p_dom,ndof);
00559
00560 r_dom = new double [ndof];
00561 nullv (r_dom,ndof);
00562
00563 z_dom = new double [ndof];
00564 nullv (z_dom,ndof);
00565
00566 d_dom = new double [ndof];
00567 nullv (d_dom,ndof);
00568
00569 x_dom = new double [ndof];
00570 nullv (x_dom,ndof);
00571
00572 aux = new double[maxnbdof];
00573
00574 buffsize=2*maxnbdof+2;
00575
00576
00577
00578
00579 buff = new double [buffsize];
00580
00581
00582 if (iv==0){
00583 for (i=0;i<ndof;i++){
00584 lhs[i]=0.0;
00585 }
00586 }
00587
00588 nory = 0.0;
00589 for(i = 0; i < ndof; i++){
00590 nory +=rhs[i];
00591 }
00592 if(mespr == 1) fprintf (stdout,"\n nory %le",nory);
00593 if(myrank == 0){
00594 norrhs = nory;
00595 for (i=1;i<nproc;i++){
00596 MPI_Recv (&nory,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00597 norrhs +=nory;
00598 }
00599 if(mespr == 1) fprintf (stdout,"\n norrhs %le",nory);
00600 }
00601 else{
00602 MPI_Send(&nory,1,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00603 }
00604 MPI_Barrier (MPI_COMM_WORLD);
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615 switch(prec){
00616 case noprec:{
00617 gm->gmxv (x_dom,p_dom);
00618 break;
00619 }
00620 case petscilu:{
00621 gm->gmxv (x_dom,p_dom);
00622
00623 break;
00624 }
00625 case pardiagprec:{
00626 gm->gmxv (x_dom,p_dom);
00627 break;
00628 }
00629 }
00630
00631
00632
00633
00634 subv (rhs,p_dom,r_dom,ndof);
00635
00636
00637
00638
00639
00640 switch(prec){
00641 case noprec:{
00642 for(i = 0; i < ndof; i++){
00643 z_dom[i] = r_dom[i];
00644 }
00645 break;
00646 }
00647 case petscilu:{
00648 solve_fact_mat_petsc(r_dom,z_dom,ndof,out);
00649
00650
00651
00652 break;
00653 }
00654 case pardiagprec:{
00655 jacobi_precondition(precvec,r_dom,z_dom,ndof,out);
00656 break;
00657 }
00658 }
00659
00660
00661
00662 copyv (z_dom,d_dom,ndof);
00663
00664
00665
00666
00667
00668
00669
00670 nom_i = v_ixu_i(r_dom,z_dom);
00671
00672
00673
00674
00675
00676 nullv(aux,maxnbdof);
00677 select_bound(aux,r_dom);
00678 for(i = 0; i < nbdof; i++){
00679 buff[i] = aux[i];
00680 }
00681
00682 nullv(aux,maxnbdof);
00683 select_bound(aux,z_dom);
00684 k = 0;
00685 for(i = maxnbdof; i < maxnbdof+nbdof; i++){
00686 buff[i] = aux[k];
00687 k++;
00688 }
00689
00690 buff[2*maxnbdof] = nom_i;
00691
00692
00693 if(myrank == 0){
00694
00695
00696
00697 r_bound = new double [ndofcp];
00698 nullv (r_bound,ndofcp);
00699
00700 z_bound = new double [ndofcp];
00701 nullv (z_bound,ndofcp);
00702
00703 d_bound = new double [ndofcp];
00704 nullv (d_bound,ndofcp);
00705
00706 x_bound = new double [ndofcp];
00707 nullv (x_bound,ndofcp);
00708
00709 p_bound = new double [ndofcp];
00710 nullv (p_bound,ndofcp);
00711
00712
00713 nom_master_i = 0.0;
00714 nom_master_i += buff[2*maxnbdof];
00715 k=domproc[0];
00716
00717 nullv (aux,maxnbdof);
00718 for(j = 0; j < nbdofdom[k]; j++){
00719 aux[j] = buff[j];
00720 }
00721 add_master_buff (aux,r_bound,k);
00722
00723 nullv (aux,maxnbdof);
00724 k = 0;
00725 for(j = maxnbdof; j < maxnbdof+nbdofdom[k]; j++){
00726 aux[j-maxnbdof] = buff[j];
00727 }
00728 add_master_buff (aux,z_bound,k);
00729
00730
00731 for (i=1;i<nproc;i++){
00732 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00733 k=domproc[stat.MPI_TAG];
00734
00735 nom_master_i += buff[2*maxnbdof];
00736 nullv (aux,maxnbdof);
00737 for(j = 0; j < nbdofdom[k]; j++){
00738 aux[j] = buff[j];
00739 }
00740 add_master_buff (aux,r_bound,k);
00741
00742 nullv (aux,maxnbdof);
00743 for(j = maxnbdof; j < maxnbdof+nbdofdom[k]; j++){
00744 aux[j-maxnbdof] = buff[j];
00745 }
00746 add_master_buff (aux,z_bound,k);
00747 }
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758 nom_master_b = ss(r_bound,z_bound,ndofcp);
00759
00760
00761
00762
00763
00764 nom = nom_master_i+nom_master_b;
00765
00766
00767
00768
00769
00770 copyv (z_bound,d_bound,ndofcp);
00771
00772
00773
00774
00775
00776
00777
00778
00779 for (i=1;i<nproc;i++){
00780 k=domproc[i];
00781 nullv (buff,maxnbdof);
00782 select_master_buff(d_bound,buff,k);
00783 MPI_Send (buff,buffsize,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
00784 }
00785
00786
00787 k=domproc[0];
00788 nullv (buff,buffsize);
00789 select_master_buff(d_bound,buff,k);
00790
00791 }
00792 else{
00793
00794 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00795 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00796
00797 }
00798 MPI_Barrier (MPI_COMM_WORLD);
00799
00800
00801 add_bound(buff,d_dom);
00802
00803
00804
00805
00806
00807 delete []buff;
00808 delete []aux;
00809
00810 buffsize=maxnbdof+2;
00811
00812
00813 buff = new double [buffsize];
00814
00815
00816
00817 titers = clock();
00818 for(i = 0; i < ni; i++){
00819
00820
00821 gm->gmxv (d_dom,p_dom);
00822
00823
00824 denom_i = v_ixu_i(d_dom,p_dom);
00825
00826
00827
00828
00829 nullv (buff,maxnbdof);
00830
00831 select_bound(buff,p_dom);
00832
00833 buff[maxnbdof] = denom_i;
00834
00835
00836 if(myrank == 0){
00837
00838 nullv (p_bound,ndofcp);
00839
00840
00841 k=domproc[0];
00842 denom_master_i = 0.0;
00843 denom_master_i += buff[maxnbdof];
00844
00845 add_master_buff (buff,p_bound,k);
00846
00847
00848 for (j=1;j<nproc;j++){
00849 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00850 k=domproc[stat.MPI_TAG];
00851
00852 denom_master_i += buff[maxnbdof];
00853
00854 add_master_buff (buff,p_bound,k);
00855 }
00856
00857
00858
00859
00860
00861
00862
00863 denom_master_b = ss(d_bound,p_bound,ndofcp);
00864
00865
00866 denom = denom_master_i+denom_master_b;
00867
00868
00869
00870
00871
00872 alpha = nom/denom;
00873 denom = nom;
00874
00875
00876 buff[maxnbdof] = alpha;
00877
00878
00879
00880
00881
00882
00883 for (j=0;j<ndofcp;j++){
00884 x_bound[j]+=alpha*d_bound[j];
00885 r_bound[j]-=alpha*p_bound[j];
00886 }
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897 for (j=1;j<nproc;j++){
00898
00899 switch(prec){
00900 case noprec:{
00901 break;
00902 }
00903 case petscilu:{
00904 k=domproc[j];
00905 nullv (buff,maxnbdof);
00906 select_master_buff(r_bound,buff,k);
00907
00908
00909
00910 break;
00911 }
00912 case pardiagprec:{
00913 k=domproc[j];
00914 nullv (buff,maxnbdof);
00915 select_master_buff(r_bound,buff,k);
00916 break;
00917 }
00918 }
00919 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
00920 }
00921
00922 switch(prec){
00923 case noprec:{
00924 break;
00925 }
00926 case petscilu:{
00927 k=domproc[0];
00928 select_master_buff(r_bound,buff,k);
00929
00930
00931
00932 break;
00933 }
00934 case pardiagprec:{
00935 k=domproc[0];
00936 nullv (buff,maxnbdof);
00937 select_master_buff(r_bound,buff,k);
00938 break;
00939 }
00940 }
00941
00942 }
00943 else{
00944 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00945 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00946 }
00947 MPI_Barrier (MPI_COMM_WORLD);
00948
00949 alpha = buff[maxnbdof];
00950
00951
00952 for (j=0;j<ndof;j++){
00953 x_dom[j]+=alpha*d_dom[j];
00954 r_dom[j]-=alpha*p_dom[j];
00955 }
00956
00957
00958
00959
00960 switch(prec){
00961 case noprec:{
00962 for(j = 0; j < ndof; j++){
00963 z_dom[j] = r_dom[j];
00964 }
00965 break;
00966 }
00967 case petscilu:{
00968
00969 add_bound(buff,r_dom);
00970
00971 ts = clock();
00972 solve_fact_mat_petsc(r_dom,z_dom,ndof,out);
00973 te = clock();
00974 if(mespr == 1) fprintf (out,"Time of back substitution %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
00975 tbacksub += ((te-ts)/(double)CLOCKS_PER_SEC);
00976
00977
00978
00979
00980 nullv (buff,maxnbdof);
00981
00982 select_bound(buff,z_dom);
00983 break;
00984 }
00985 case pardiagprec:{
00986
00987 add_bound(buff,r_dom);
00988
00989 ts = clock();
00990 jacobi_precondition(precvec,r_dom,z_dom,ndof,out);
00991 te = clock();
00992 if(mespr == 1) fprintf (out,"Time of precondition %le s\n",(te-ts)/(double)CLOCKS_PER_SEC);
00993 tbacksub += ((te-ts)/(double)CLOCKS_PER_SEC);
00994
00995
00996
00997
00998 nullv (buff,maxnbdof);
00999
01000 select_bound(buff,z_dom);
01001 break;
01002 }
01003 }
01004
01005
01006
01007 nom_i = v_ixu_i(r_dom,z_dom);
01008 buff[maxnbdof] = nom_i;
01009
01010
01011 if(myrank == 0){
01012
01013 nom_master_i = buff[maxnbdof];
01014 k=domproc[0];
01015 nullv(z_bound,ndofcp);
01016 switch(prec){
01017 case noprec:{
01018 for(j = 0; j < ndofcp; j++){
01019 z_bound[j] = r_bound[j];
01020 }
01021 break;
01022 }
01023 case petscilu:{
01024
01025 add_master_buff (buff,z_bound,k);
01026 break;
01027 }
01028 case pardiagprec:{
01029
01030 add_master_buff (buff,z_bound,k);
01031 break;
01032 }
01033
01034 }
01035
01036
01037 for (j=1;j<nproc;j++){
01038 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01039 k=domproc[stat.MPI_TAG];
01040 nom_master_i += buff[maxnbdof];
01041 switch(prec){
01042 case noprec:{
01043 break;
01044 }
01045 case petscilu:{
01046
01047 add_master_buff (buff,z_bound,k);
01048 break;
01049 }
01050 case pardiagprec:{
01051
01052 add_master_buff (buff,z_bound,k);
01053 break;
01054 }
01055 }
01056
01057 }
01058
01059
01060
01061
01062
01063
01064
01065 nom_master_b = ss(r_bound,z_bound,ndofcp);
01066
01067
01068 nom = nom_master_i+nom_master_b;
01069
01070
01071
01072 if(sqrt(nom/norrhs) < err || i == ni-1){
01073
01074 buff[maxnbdof+1]=10000.0;
01075 iter = i;
01076 }
01077
01078
01079 if(mespr == 1) fprintf (stdout,"\n iteration %ld norres %le norrhs %le norres/norrhs %e",i,sqrt(nom),sqrt(norrhs),sqrt(nom/norrhs));
01080 beta = nom/denom;
01081
01082
01083
01084
01085
01086 for (j=0;j<ndofcp;j++){
01087 d_bound[j]=beta*d_bound[j]+z_bound[j];
01088 }
01089
01090
01091
01092
01093
01094
01095 for (j=1;j<nproc;j++){
01096 k=domproc[j];
01097 nullv (buff,maxnbdof);
01098
01099 select_master_buff(d_bound,buff,k);
01100
01101 buff[maxnbdof]=beta;
01102
01103 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
01104 }
01105
01106 k=domproc[0];
01107 nullv (buff,maxnbdof);
01108
01109 select_master_buff(d_bound,buff,k);
01110
01111 buff[maxnbdof] = beta;
01112
01113 }
01114 else{
01115 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01116 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01117 }
01118
01119
01120 add_bound(buff,d_dom);
01121
01122 beta = buff[maxnbdof];
01123
01124
01125
01126
01127 for (j=0;j<nidof;j++){
01128 d_dom[j]=beta*d_dom[j]+z_dom[j];
01129 }
01130
01131 if (buff[maxnbdof+1]==10000.0){
01132 if (myrank==0){
01133
01134 for (j=1;j<nproc;j++){
01135 k=domproc[j];
01136 nullv (buff,maxnbdof);
01137 select_master_buff(x_bound,buff,k);
01138 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
01139 }
01140
01141 k=domproc[0];
01142 nullv (buff,buffsize-1);
01143 select_master_buff(x_bound,buff,k);
01144
01145
01146 delete []r_bound;
01147 delete []d_bound;
01148 delete []x_bound;
01149 delete []p_bound;
01150 delete []z_bound;
01151
01152 }
01153 else{
01154 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01155 }
01156
01157 add_bound(buff,x_dom);
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171 copyv (x_dom,lhs,ndof);
01172
01173
01174 delete []r_dom;
01175 delete []d_dom;
01176 delete []x_dom;
01177 delete []p_dom;
01178 delete []z_dom;
01179
01180 break;
01181 }
01182
01183 }
01184 titere = clock();
01185 MPI_Barrier (MPI_COMM_WORLD);
01186
01187
01188 tsole = clock();
01189 if(myrank == 0){
01190 printf ("\n\nTime of computation\n\n\n");
01191 switch(prec){
01192 case noprec:{
01193 if(mespr == 1) fprintf (out,"Time of iteration %le s\n",(titere-titers)/(double)CLOCKS_PER_SEC);
01194 if(mespr == 1) fprintf (out,"Time of solution of system %le s\n",(tsole-tsols)/(double)CLOCKS_PER_SEC);
01195 printf ("Time of iteration %le s\n",(titere-titers)/(double)CLOCKS_PER_SEC);
01196 printf ("Time of solution of system %le s\n",(tsole-tsols)/(double)CLOCKS_PER_SEC);
01197 break;
01198 }
01199 case petscilu:{
01200 if(mespr == 1) fprintf (out,"Time of ILU factorisation by PETSC %le s\n",(tfacte-tfacts)/(double)CLOCKS_PER_SEC);
01201 if(mespr == 1) fprintf (out,"Time of iteration %le s\n",(titere-titers)/(double)CLOCKS_PER_SEC);
01202 if(mespr == 1) fprintf (out,"Time of solution of system %le s\n",(tsole-tsols)/(double)CLOCKS_PER_SEC);
01203 if(mespr == 1) fprintf (out,"Domain 0:\n");
01204 if(mespr == 1) fprintf (out,"Time of all back substitution %le s\n",tbacksub);
01205 if(mespr == 1) fprintf (out,"Time of back substitution per iteration %lf s\n",tbacksub/double(iter));
01206
01207 printf ("Time of ILU factorisation by PETSC %le s\n",(tfacte-tfacts)/(double)CLOCKS_PER_SEC);
01208 printf ("Time of iteration %le s\n",(titere-titers)/(double)CLOCKS_PER_SEC);
01209 printf ("Time of solution of system %le s\n",(tsole-tsols)/(double)CLOCKS_PER_SEC);
01210 printf ("Domain 0:\n");
01211 printf ("Time of all back substitution %le s\n",tbacksub);
01212 printf ("Time of back substitution per iteration %le s\n",tbacksub/double(iter));
01213
01214
01215
01216 for (i=1;i<nproc;i++){
01217 MPI_Recv (&tbacksub,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01218 j=domproc[stat.MPI_TAG];
01219 if(mespr == 1) fprintf (out,"Domain %ld:\n",j);
01220 if(mespr == 1) fprintf (out,"Time of all back substitution %le s\n",tbacksub);
01221 if(mespr == 1) fprintf (out,"Time of back substitution per iteration %lf s\n",tbacksub/double(iter));
01222 printf ("Domain %ld:\n",j);
01223 printf ("Time of all back substitution %le s\n",tbacksub);
01224 printf ("Time of back substitution per iteration %le s\n",tbacksub/double(iter));
01225 }
01226 break;
01227 }
01228 }
01229 }
01230 else{
01231 MPI_Send(&tbacksub,1,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01232 }
01233 MPI_Barrier (MPI_COMM_WORLD);
01234
01235 delete []buff;
01236 switch(prec){
01237 case pardiagprec:{
01238 delete []precvec;
01239 break;
01240 }
01241 }
01242 }
01243