00001 #include "schurcompl.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <time.h>
00005 #include <limits.h>
00006
00007 #include "precond.h"
00008 #include "iotools.h"
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 schurcompl::schurcompl(int np,int mr,long nd)
00020 {
00021
00022 nproc=np;
00023
00024 myrank=mr;
00025
00026 ndom=nd;
00027
00028
00029 ndof=0;
00030
00031 nidof=0;
00032
00033 nbdof=0;
00034
00035 maxnbdof=0;
00036
00037
00038
00039
00040
00041 trssol = (redsystsolver) 0;
00042
00043 rsmstor = (storagetype) 0;
00044
00045 nicgsch=0;
00046
00047 errcgsch=0.0;
00048
00049
00050
00051
00052
00053 anicgsch=0;
00054
00055 aerrcgsch=0.0;
00056
00057 zero=1.0e-15;
00058
00059
00060 nbdofmas=NULL;
00061
00062
00063 gtop=NULL;
00064 arr=NULL;
00065
00066
00067 ssle = new slesolv ();
00068
00069 destr_nbdofmas = 0;
00070 }
00071
00072 schurcompl::~schurcompl()
00073 {
00074 if (myrank==0){
00075 if (destr_nbdofmas)
00076 delete [] nbdofmas;
00077
00078 delete gtop;
00079 delete arr;
00080 }
00081 delete ssle;
00082 }
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366 void schurcompl::initiate (partop *ptop,selnodes *selnodschur)
00367 {
00368
00369 ndof = ptop->ndof;
00370
00371 nidof = ptop->nidof;
00372
00373 nbdof = ptop->nbdof;
00374
00375
00376 maxnbdof = selnodschur->maxsnndof;
00377
00378 if (myrank==0){
00379
00380
00381 ndofcp = selnodschur->tndofsn;
00382
00383
00384 nbdofmas = selnodschur->snndofmas;
00385
00386 gtop = new gtopology;
00387 gtop->initiate (selnodschur->cnmas,selnodschur->snndofmas,nproc);
00388 gtop->cnstate=1;
00389 }
00390 }
00391
00392 void schurcompl::initiate (partopjk *ptop,FILE *out)
00393 {
00394 long i;
00395
00396
00397 ndof = ptop->ndof;
00398
00399 nidof = ptop->nidof;
00400
00401 nbdof = ptop->nbdof;
00402
00403
00404 maxnbdof = ptop->maxnbdofd;
00405
00406 fprintf (out,"\n\n jouda \n\n");
00407 fprintf (out,"\n ndof %ld",ndof);
00408 fprintf (out,"\n nidof %ld",nidof);
00409 fprintf (out,"\n nbdof %ld",nbdof);
00410 fprintf (out,"\n maxnbdof %ld",maxnbdof);
00411
00412 if (myrank==0){
00413
00414
00415 ndofcp = ptop->ndofc;
00416
00417 fprintf (out,"\n ndofcp %ld",ndofcp);
00418
00419
00420 nbdofmas = ptop->nbdofd;
00421
00422 for (i=0;i<nproc;i++){
00423 fprintf (out,"\n nbdofmas %ld %ld",i,nbdofmas[i]);
00424 }
00425
00426 gtop = new gtopology;
00427 gtop->initiate (ptop->bdofd,ptop->nbdofd,nproc);
00428 gtop->cnstate=1;
00429 }
00430
00431 }
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599 void schurcompl::solve_red_sys_iter (long *domproc,double *condmat,double *condvect,FILE *out)
00600 {
00601 int ind;
00602 long i,j,k,sizebuff,gndofe,*cn;
00603 double alpha,beta,nom,denom,norrhs;
00604 double *lhs,*rhs,*buff,*p,*d,*r;
00605 MPI_Status stat;
00606 rhs=NULL;
00607 lhs=NULL;
00608 sizebuff = maxnbdof+1;
00609 buff = new double [sizebuff];
00610 memset (buff,0,sizebuff*sizeof(double));
00611
00612
00613 if (myrank==0){
00614 lhs = new double [ndofcp];
00615 memset (lhs,0,ndofcp*sizeof(double));
00616 rhs = new double [ndofcp];
00617 memset (rhs,0,ndofcp*sizeof(double));
00618
00619 gndofe = gtop->give_ndofe (domproc[0]);
00620 cn = new long [gndofe];
00621 gtop->give_code_numbers (domproc[0],cn);
00622 locglob (rhs,condvect,cn,gndofe);
00623 delete [] cn;
00624
00625 for (i=1;i<nproc;i++){
00626 MPI_Recv(buff,sizebuff,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00627 j=domproc[stat.MPI_TAG]; ind=0;
00628
00629 gndofe = gtop->give_ndofe (j);
00630 cn = new long [gndofe];
00631 gtop->give_code_numbers (j,cn);
00632 locglob (rhs,buff,cn,gndofe);
00633 delete [] cn;
00634
00635 }
00636 }
00637 else{
00638 copyv (condvect,buff,ndof-nidof);
00639 MPI_Send(buff,sizebuff,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00640 }
00641
00642
00643 if (myrank==0){
00644 r = new double [ndofcp];
00645 d = new double [ndofcp];
00646 p = new double [ndofcp];
00647
00648
00649 norrhs=0.0; nom=0.0;
00650 for (i=0;i<ndofcp;i++){
00651 lhs[i]=0.0;
00652 r[i]=rhs[i];
00653 d[i]=rhs[i];
00654 norrhs+=rhs[i]*rhs[i];
00655 nom+=r[i]*r[i];
00656 }
00657
00658
00659 for (i=0;i<nicgsch;i++){
00660
00661
00662 for (j=1;j<nproc;j++){
00663 memset (buff,0,sizebuff*sizeof(double));
00664 gndofe = gtop->give_ndofe (domproc[j]);
00665 cn = new long [gndofe];
00666 gtop->give_code_numbers (domproc[j],cn);
00667 globloc (d,buff,cn,gndofe);
00668 delete [] cn;
00669
00670 MPI_Send(buff,sizebuff,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
00671 }
00672
00673 memset (buff,0,sizebuff*sizeof(double));
00674 gndofe = gtop->give_ndofe (domproc[0]);
00675 cn = new long [gndofe];
00676 gtop->give_code_numbers (domproc[0],cn);
00677 globloc (d,buff,cn,gndofe);
00678 delete [] cn;
00679
00680
00681 mxv (condmat,buff,condvect,nbdof,nbdof);
00682
00683
00684 memset (p,0,ndofcp*sizeof(double));
00685
00686 gndofe = gtop->give_ndofe (domproc[0]);
00687 cn = new long [gndofe];
00688 gtop->give_code_numbers (domproc[0],cn);
00689 locglob (p,condvect,cn,gndofe);
00690 delete [] cn;
00691
00692 for (j=1;j<nproc;j++){
00693 MPI_Recv(buff,sizebuff,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00694 k=domproc[stat.MPI_TAG];
00695
00696 gndofe = gtop->give_ndofe (k);
00697 cn = new long [gndofe];
00698 gtop->give_code_numbers (k,cn);
00699 locglob (p,buff,cn,gndofe);
00700 delete [] cn;
00701
00702 }
00703
00704
00705 denom=ss(d,p,ndofcp);
00706
00707 if (fabs(denom)<zero){
00708 if (i!=nicgsch-1){
00709 buff[maxnbdof]=1.0;
00710 for (j=1;j<nproc;j++){
00711 MPI_Send(buff,sizebuff,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
00712 }
00713 }
00714 fprintf (stderr,"\n\n denominator in alpha expression in parallel conjugate gradient method is equal to zero.\n");
00715 break;
00716 }
00717
00718
00719 alpha=nom/denom;
00720
00721
00722 for (j=0;j<ndofcp;j++){
00723 r[j]-=alpha*p[j];
00724 lhs[j]+=alpha*d[j];
00725 }
00726
00727
00728 denom=nom;
00729 nom=ss(r,r,ndofcp);
00730
00731
00732 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,nom/norrhs);
00733 if (nom/norrhs<errcgsch){
00734 fprintf (stdout,"\n last iteration %ld norres/norrhs %e",i,nom/norrhs);
00735 if (i!=nicgsch-1){
00736 buff[maxnbdof]=1.0;
00737 for (j=1;j<nproc;j++){
00738 MPI_Send(buff,sizebuff,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
00739 }
00740 }
00741 break;
00742 }
00743
00744
00745 beta=nom/denom;
00746
00747
00748 for (j=0;j<ndofcp;j++){
00749 d[j]=r[j]+beta*d[j];
00750 }
00751
00752 }
00753 anicgsch=i;
00754 aerrcgsch=nom/norrhs;
00755
00756
00757
00758
00759
00760
00761
00762
00763 }
00764 else{
00765
00766 for (i=0;i<nicgsch;i++){
00767
00768 MPI_Recv(buff,sizebuff,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00769
00770 if (buff[maxnbdof]>0.5) break;
00771
00772
00773 mxv (condmat,buff,condvect,nbdof,nbdof);
00774
00775
00776 copyv (condvect,buff,nbdof);
00777 MPI_Send(buff,sizebuff,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00778
00779 }
00780 anicgsch=i;
00781 }
00782
00783
00784 if (myrank==0){
00785 for (i=1;i<nproc;i++){
00786
00787 memset (buff,0,sizebuff*sizeof(double));
00788 gndofe = gtop->give_ndofe (domproc[i]);
00789 cn = new long [gndofe];
00790 gtop->give_code_numbers (domproc[i],cn);
00791 globloc (lhs,buff,cn,gndofe);
00792 delete [] cn;
00793
00794 MPI_Send(buff,sizebuff,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
00795 }
00796
00797 nullv (condvect,nbdofmas[domproc[0]]);
00798
00799 gndofe = gtop->give_ndofe (domproc[0]);
00800 cn = new long [gndofe];
00801 gtop->give_code_numbers (domproc[0],cn);
00802 globloc (lhs,buff,cn,gndofe);
00803 delete [] cn;
00804
00805 }
00806 else{
00807 MPI_Recv(buff,sizebuff,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00808 }
00809 copyv (buff,condvect,nbdof);
00810
00811 delete [] buff;
00812
00813 }
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828 void schurcompl::solve_red_sys_fin (long *domproc,double *condmat,double *condvect,long decomp,FILE *out)
00829 {
00830 long i,j,k,l,ii,gndofe,*cn;
00831 double *lhs,*rhs;
00832 double *buff;
00833 MPI_Status stat;
00834 precond prec;
00835
00836 if (decomp==0)
00837 buff = new double [maxnbdof*maxnbdof+maxnbdof];
00838 else
00839 buff = new double [maxnbdof];
00840
00841 ii=0;
00842 if (decomp==0)
00843 {
00844 for (i=0;i<nbdof;i++){
00845 for (j=0;j<nbdof;j++){
00846 buff[ii]=condmat[i*nbdof+j]; ii++;
00847 }
00848 }
00849 }
00850 for (i=0;i<nbdof;i++){
00851 buff[ii]=condvect[i]; ii++;
00852 }
00853
00854 if (myrank==0){
00855 lhs = new double [ndofcp];
00856 memset (lhs,0,ndofcp*sizeof(double));
00857 rhs = new double [ndofcp];
00858 memset (rhs,0,ndofcp*sizeof(double));
00859 if (decomp==0){
00860 if (arr){
00861
00862 delete arr;
00863 }
00864 arr = new gmatrix;
00865 arr->ts = rsmstor;
00866 arr->setval (ssle);
00867 arr->alloc ();
00868 arr->initiate (gtop,ndofcp,rsmstor,1,out);
00869
00870 gndofe = gtop->give_ndofe (domproc[0]);
00871 cn = new long [gndofe];
00872 gtop->give_code_numbers (domproc[0],cn);
00873 arr->localized (condmat,cn,gndofe,domproc[0]);
00874 locglob (rhs,condvect,cn,gndofe);
00875 delete [] cn;
00876
00877
00878 for (i=1;i<nproc;i++){
00879 MPI_Recv(buff,maxnbdof*maxnbdof+maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00880 l=domproc[stat.MPI_TAG];
00881
00882 ii=0;
00883 for (j=0;j<nbdofmas[l];j++){
00884 for (k=0;k<nbdofmas[l];k++){
00885 condmat[j*nbdofmas[l]+k]=buff[ii]; ii++;
00886 }
00887 }
00888 for (j=0;j<nbdofmas[l];j++){
00889 condvect[j]=buff[ii]; ii++;
00890 }
00891
00892 gndofe = gtop->give_ndofe (l);
00893 cn = new long [gndofe];
00894 gtop->give_code_numbers (l,cn);
00895 arr->localized (condmat,cn,gndofe,l);
00896 locglob (rhs,condvect,cn,gndofe);
00897 delete [] cn;
00898
00899 }
00900
00901 arr->prepmat (0.0,1);
00902 }
00903 else
00904 {
00905 gndofe = gtop->give_ndofe (domproc[0]);
00906 cn = new long [gndofe];
00907 gtop->give_code_numbers (domproc[0],cn);
00908 locglob (rhs,condvect,cn,gndofe);
00909 delete [] cn;
00910 for (i=1;i<nproc;i++){
00911 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00912 l=domproc[stat.MPI_TAG];
00913
00914 ii=0;
00915 for (j=0;j<nbdofmas[l];j++){
00916 condvect[j]=buff[ii]; ii++;
00917 }
00918
00919 gndofe = gtop->give_ndofe (l);
00920 cn = new long [gndofe];
00921 gtop->give_code_numbers (l,cn);
00922 locglob (rhs,condvect,cn,gndofe);
00923 delete [] cn;
00924 }
00925 }
00926
00927
00928
00929
00930
00931
00932
00933 ssle->solve_system (gtop,arr,lhs,rhs,out);
00934
00935
00936
00937
00938
00939
00940
00941
00942 for (i=1;i<nproc;i++){
00943 nullv (condvect,nbdofmas[domproc[i]]);
00944
00945 gndofe = gtop->give_ndofe (domproc[i]);
00946 cn = new long [gndofe];
00947 gtop->give_code_numbers (domproc[i],cn);
00948 globloc (lhs,condvect,cn,gndofe);
00949 delete [] cn;
00950
00951 for (long ij=0;ij<nbdofmas[i];ij++){
00952 buff[ij]=condvect[ij];
00953 }
00954 MPI_Send(buff,maxnbdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
00955
00956 }
00957
00958 nullv (condvect,nbdofmas[ndom]);
00959
00960 gndofe = gtop->give_ndofe (domproc[0]);
00961 cn = new long [gndofe];
00962 gtop->give_code_numbers (domproc[0],cn);
00963 globloc (lhs,condvect,cn,gndofe);
00964 delete [] cn;
00965
00966 delete [] lhs; delete [] rhs;
00967
00968 }
00969 else{
00970 if (decomp == 0)
00971 MPI_Send(buff,maxnbdof*maxnbdof+maxnbdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00972 else
00973 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
00974
00975 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00976
00977 for (long ij=0;ij<nbdof;ij++){
00978 condvect[ij]=buff[ij];
00979 }
00980 }
00981
00982 delete [] buff;
00983
00984 }
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998 void schurcompl::solve_red_sys (long *domproc,double *condmat,double *condvect,long decomp,FILE *out)
00999 {
01000 switch (trssol){
01001 case master_sol:{
01002 solve_red_sys_fin (domproc,condmat,condvect,decomp,out);
01003 break;
01004 }
01005 case paral_sol:{
01006 solve_red_sys_iter (domproc,condmat,condvect,out);
01007 break;
01008 }
01009 default:{
01010 par_print_err(myrank,"unknown type of solver of reduced problem is required",__FILE__,__LINE__,__func__);
01011 }
01012 }
01013 }
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029 void schurcompl::solve_system (gtopology *top,gmatrix *gm,
01030 long *domproc,double *lhs,double *rhs,FILE *out)
01031 {
01032 long i,j;
01033 double *condmat,*condvect;
01034 time_t t1,t2,t3,t4;
01035 long decomp;
01036
01037
01038 if (myrank==0){
01039 condmat = new double [maxnbdof*maxnbdof];
01040 memset (condmat,0,maxnbdof*maxnbdof*sizeof(double));
01041 }
01042 else{
01043 condmat = new double [nbdof*nbdof];
01044 memset (condmat,0,nbdof*nbdof*sizeof(double));
01045 }
01046
01047
01048 if (myrank==0){
01049 condvect = new double [maxnbdof];
01050 memset (condvect,0,maxnbdof*sizeof(double));
01051 }
01052 else{
01053 condvect = new double [nbdof];
01054 memset (condvect,0,nbdof*sizeof(double));
01055 }
01056
01057 t1 = time (NULL);
01058
01059 decomp = gm->decomp();
01060
01061
01062 gm->condense (top,condmat,condvect,lhs,rhs,nbdof,1,out);
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105 t2 = time (NULL);
01106
01107
01108
01109 solve_red_sys (domproc,condmat,condvect,decomp,out);
01110
01111
01112 fprintf (out,"\n\n\n kontrola po hrubem problemu \n");
01113
01114 for (i=0;i<nbdof;i++){
01115
01116 fprintf (out,"\n condvect %le ",condvect[i]);
01117 }
01118
01119
01120
01121
01122 t3 = time (NULL);
01123
01124
01125 gm->condense (top,condmat,condvect,lhs,rhs,nbdof,2,out);
01126
01127
01128
01129
01130
01131 t4 = time (NULL);
01132
01133
01134
01135
01136 MPI_Status stat;
01137 long *buff;
01138 buff = new long[3];
01139 buff[0]=t2-t1;
01140 buff[1]=t3-t2;
01141 buff[2]=t4-t3;
01142
01143
01144 if (myrank==0){
01145
01146 long *fact,*cp,*tottime;
01147 fact=new long [nproc];
01148 cp=new long [nproc];
01149 tottime = new long [nproc];
01150
01151 fprintf (out,"\n\n\n\n\n");
01152 j=domproc[myrank];
01153 fact[j]=buff[0];
01154 cp[j]=buff[1];
01155 tottime[j]=buff[0]+buff[1]+buff[2];
01156
01157
01158
01159
01160
01161
01162
01163
01164 j=domproc[myrank];
01165
01166
01167
01168
01169
01170 for (i=1;i<nproc;i++){
01171 MPI_Recv (buff,3,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01172 j=domproc[stat.MPI_TAG];
01173 fact[j]=buff[0];
01174 cp[j]=buff[1];
01175 tottime[j]=buff[0]+buff[1]+buff[2];
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187 }
01188
01189 long maxfact,mincp,mintot,maxtot;
01190 maxfact=0;
01191 maxtot=0;
01192 mincp=LONG_MAX;
01193 mintot=LONG_MAX;
01194 for (i=0;i<nproc;i++){
01195 if (maxfact<fact[i])
01196 maxfact=fact[i];
01197
01198 if (mincp>cp[i])
01199 mincp=cp[i];
01200
01201 if (mintot>tottime[i])
01202 mintot=tottime[i];
01203 if (maxtot<tottime[i])
01204 maxtot=tottime[i];
01205 }
01206
01207
01208
01209
01210
01211
01212 delete [] fact;
01213 delete [] cp;
01214 delete [] tottime;
01215 }
01216 else{
01217 MPI_Send(buff,3,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01218 }
01219 MPI_Barrier (MPI_COMM_WORLD);
01220
01221 delete [] buff;
01222
01223
01224 delete [] condmat;
01225 delete [] condvect;
01226
01227 }
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242 void schurcompl::gather_bound_vect (double *lv,double *gv,long *domproc)
01243 {
01244 long i,j,l,gndofe,*cn;
01245 double *buff;
01246 MPI_Status stat;
01247
01248 buff = new double [maxnbdof];
01249
01250
01251 j=0;
01252 for (i=nidof;i<ndof;i++){
01253 buff[j]=lv[i]; j++;
01254 }
01255
01256 if (myrank==0){
01257
01258
01259
01260 l=domproc[0];
01261
01262 gndofe = gtop->give_ndofe (l);
01263 cn = new long [gndofe];
01264
01265 gtop->give_code_numbers (l,cn);
01266
01267 locglob (gv,buff,cn,gndofe);
01268 delete [] cn;
01269
01270
01271
01272
01273 for (i=1;i<nproc;i++){
01274 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01275 l=domproc[stat.MPI_TAG];
01276
01277
01278 gndofe = gtop->give_ndofe (l);
01279 cn = new long [gndofe];
01280
01281 gtop->give_code_numbers (l,cn);
01282
01283 locglob (gv,buff,cn,gndofe);
01284 delete [] cn;
01285 }
01286 }
01287 else{
01288 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01289 }
01290 MPI_Barrier (MPI_COMM_WORLD);
01291
01292 delete [] buff;
01293 }
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307 double schurcompl::pss_gather_bound_vect (double *lv1,double *gv1,double *lv2,double *gv2,long *domproc)
01308 {
01309 long i,j,l,gndofe,*cn,buffsize;
01310 double dp,idp,*buff;
01311 MPI_Status stat;
01312
01313
01314
01315 buffsize=2*maxnbdof+1;
01316 buff = new double [buffsize];
01317
01318 idp=0.0;
01319
01320
01321 dp=0.0;
01322 for (i=0;i<nidof;i++){
01323 dp+=lv1[i]*lv2[i];
01324 }
01325 buff[buffsize-1]=dp;
01326
01327
01328 j=0;
01329 for (i=nidof;i<ndof;i++){
01330 buff[j]=lv1[i];
01331 buff[j+maxnbdof]=lv2[i]; j++;
01332 }
01333
01334 if (myrank==0){
01335
01336 idp=buff[buffsize-1];
01337
01338
01339
01340
01341 l=domproc[0];
01342
01343 gndofe = gtop->give_ndofe (l);
01344 cn = new long [gndofe];
01345
01346 gtop->give_code_numbers (l,cn);
01347
01348 locglob (gv1,buff,cn,gndofe);
01349 locglob (gv2,buff+maxnbdof,cn,gndofe);
01350 delete [] cn;
01351
01352
01353
01354
01355 for (i=1;i<nproc;i++){
01356 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01357 l=domproc[stat.MPI_TAG];
01358
01359
01360 gndofe = gtop->give_ndofe (l);
01361 cn = new long [gndofe];
01362
01363 gtop->give_code_numbers (l,cn);
01364
01365 locglob (gv1,buff,cn,gndofe);
01366 locglob (gv2,buff+maxnbdof,cn,gndofe);
01367 delete [] cn;
01368
01369
01370 idp+=buff[buffsize-1];
01371 }
01372 }
01373 else{
01374 MPI_Send(buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01375 }
01376 MPI_Barrier (MPI_COMM_WORLD);
01377
01378 delete [] buff;
01379
01380 return idp;
01381 }
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405 long schurcompl::selected_norm_calculation (long cid,double err,double thresh,partopjk *ptop,gtopology *top,double *lv1,double *gv1,double *rhs,double *grhs,long *domproc,long n,FILE *out)
01406 {
01407 long i,j,k,l,nid,dof,gndofe,*cn,buffsize,stop,ndofn;
01408 double dp,totdp,norrhs,totnorrhs,*buff;
01409 MPI_Status stat;
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423 buffsize=2*maxnbdof+2;
01424 buff = new double [buffsize];
01425 nullv (buff,buffsize);
01426
01427
01428 dp=0.0;
01429
01430 norrhs=0.0;
01431
01432
01433 for (i=0;i<ptop->nin;i++){
01434
01435 nid=ptop->inid[i];
01436
01437
01438 dof = top->give_dof (nid,cid);
01439
01440 if (dof>0){
01441 dof--;
01442 dp+=lv1[dof]*lv1[dof];
01443 norrhs+=rhs[dof]*rhs[dof];
01444 }
01445
01446 }
01447
01448 fprintf (stdout,"\n dp %le",dp);
01449
01450 buff[buffsize-2]=dp;
01451 buff[buffsize-1]=norrhs;
01452
01453
01454 k=0;
01455
01456 for (i=0;i<ptop->nbn;i++){
01457
01458 nid=ptop->bnid[i];
01459
01460
01461
01462 ndofn = top->give_ndofn (nid);
01463 for (j=0;j<ndofn;j++){
01464 dof = top->give_dof (nid,j);
01465 if (dof>0){
01466 if (cid==j){
01467 dof--;
01468 buff[k]=lv1[dof];
01469 buff[maxnbdof+k]=rhs[dof];
01470 }
01471 k++;
01472 }
01473 }
01474 }
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484 if (myrank==0){
01485
01486 totdp=buff[buffsize-2];
01487 totnorrhs=buff[buffsize-1];
01488
01489
01490
01491
01492 l=domproc[0];
01493
01494 gndofe = gtop->give_ndofe (l);
01495 cn = new long [gndofe];
01496
01497 gtop->give_code_numbers (l,cn);
01498
01499 fprintf (out,"\n\n domain %ld\n kontrola kodovych cisel na hranici\n",myrank);
01500 for (long ijk=0;ijk<gndofe;ijk++){
01501 fprintf (out," %ld",cn[ijk]);
01502 }
01503 fprintf (out,"\n\n");
01504
01505
01506
01507 locglob (gv1,buff,cn,gndofe);
01508 locglob (grhs,buff+maxnbdof,cn,gndofe);
01509 delete [] cn;
01510
01511
01512
01513
01514
01515
01516
01517
01518 for (i=1;i<nproc;i++){
01519 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01520
01521 l=domproc[stat.MPI_TAG];
01522
01523
01524 gndofe = gtop->give_ndofe (l);
01525 cn = new long [gndofe];
01526
01527 gtop->give_code_numbers (l,cn);
01528
01529 fprintf (out,"\n\n domain %ld\n kontrola kodovych cisel na hranici\n",l);
01530 for (long ijk=0;ijk<gndofe;ijk++){
01531 fprintf (out," %ld",cn[ijk]);
01532 }
01533 fprintf (out,"\n\n");
01534
01535
01536 locglob (gv1,buff,cn,gndofe);
01537 locglob (grhs,buff+maxnbdof,cn,gndofe);
01538 delete [] cn;
01539
01540
01541
01542
01543
01544
01545 totdp+=buff[buffsize-2];
01546 totnorrhs+=buff[buffsize-1];
01547 }
01548
01549 }
01550 else{
01551 MPI_Send(buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01552 }
01553 MPI_Barrier (MPI_COMM_WORLD);
01554
01555
01556 delete [] buff;
01557
01558
01559 stop=0;
01560 if (myrank==0){
01561
01562 fprintf (out,"\n Kontrola gv na masterovi");
01563 for (i=0;i<ndofcp;i++){
01564 fprintf (out,"\n gv1 %6ld % 12.9le grhs % 12.9le",i,gv1[i],grhs[i]);
01565 totdp+=gv1[i]*gv1[i];
01566 totnorrhs+=grhs[i]*grhs[i];
01567 }
01568
01569
01570 totdp=sqrt(totdp);
01571 totnorrhs=sqrt(totnorrhs);
01572
01573 fprintf (stdout,"\n\n\n\n component ID %ld",cid);
01574 fprintf (stdout,"\n totdp %le",totdp);
01575 fprintf (stdout,"\n totnorrhs %le",totnorrhs);
01576 fprintf (stdout,"\n thresh %le",thresh);
01577 fprintf (stdout,"\n err %le",err);
01578 fprintf (stdout,"\n\n\n\n\n");
01579
01580
01581 if (totnorrhs<thresh){
01582 if (totdp<err){
01583
01584 stop=1;
01585 }
01586 }
01587 else{
01588 if (totdp/totnorrhs<err){
01589
01590 stop=1;
01591 }
01592 }
01593
01594 for (i=1;i<nproc;i++){
01595 MPI_Send(&stop,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01596 }
01597
01598 }else{
01599 MPI_Recv(&stop,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01600 }
01601 MPI_Barrier (MPI_COMM_WORLD);
01602
01603 return stop;
01604 }
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617 double schurcompl::pss_gather_bound_vect_old (double *lv,double *gv,long *domproc)
01618 {
01619 long i,j,l,gndofe,*cn,buffsize;
01620 double dp,idp,*buff;
01621 MPI_Status stat;
01622
01623 buffsize=maxnbdof+1;
01624 buff = new double [buffsize];
01625
01626 idp=0.0;
01627
01628
01629 dp=0.0;
01630 for (i=0;i<nidof;i++){
01631 dp+=lv[i]*lv[i];
01632 }
01633 buff[maxnbdof]=dp;
01634
01635
01636 j=0;
01637 for (i=nidof;i<ndof;i++){
01638 buff[j]=lv[i]; j++;
01639 }
01640
01641 if (myrank==0){
01642
01643 idp=buff[maxnbdof];
01644
01645
01646
01647
01648 l=domproc[0];
01649
01650 gndofe = gtop->give_ndofe (l);
01651 cn = new long [gndofe];
01652
01653 gtop->give_code_numbers (l,cn);
01654
01655 locglob (gv,buff,cn,gndofe);
01656 delete [] cn;
01657
01658
01659
01660
01661 for (i=1;i<nproc;i++){
01662 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01663 l=domproc[stat.MPI_TAG];
01664
01665
01666 gndofe = gtop->give_ndofe (l);
01667 cn = new long [gndofe];
01668
01669 gtop->give_code_numbers (l,cn);
01670
01671 locglob (gv,buff,cn,gndofe);
01672 delete [] cn;
01673
01674
01675 idp+=buff[maxnbdof];
01676 }
01677 }
01678 else{
01679 MPI_Send(buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01680 }
01681 MPI_Barrier (MPI_COMM_WORLD);
01682
01683 delete [] buff;
01684
01685 return idp;
01686 }
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699 void schurcompl::scatter_bound_vect (double *lv,double *gv,long *domproc)
01700 {
01701 long i,j,gndofe,*cn;
01702 double *buff;
01703 MPI_Status stat;
01704
01705 buff = new double [maxnbdof];
01706
01707 if (myrank==0){
01708 for (i=1;i<nproc;i++){
01709 nullv (buff,nbdofmas[domproc[i]]);
01710
01711 gndofe = gtop->give_ndofe (domproc[i]);
01712 cn = new long [gndofe];
01713 gtop->give_code_numbers (domproc[i],cn);
01714 globloc (gv,buff,cn,gndofe);
01715 delete [] cn;
01716
01717 MPI_Send(buff,maxnbdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01718 }
01719
01720 nullv (buff,nbdofmas[domproc[0]]);
01721
01722 gndofe = gtop->give_ndofe (domproc[0]);
01723 cn = new long [gndofe];
01724 gtop->give_code_numbers (domproc[0],cn);
01725 globloc (gv,buff,cn,gndofe);
01726 delete [] cn;
01727
01728 }
01729 else{
01730 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01731 }
01732 MPI_Barrier (MPI_COMM_WORLD);
01733
01734 j=0;
01735 for (i=nidof;i<ndof;i++){
01736 lv[i]=buff[j]/2.0; j++;
01737 }
01738
01739 delete [] buff;
01740
01741 }
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755 double schurcompl::unbalanced_values (double *lv,long *domproc,FILE *out)
01756 {
01757 long i;
01758 double aux,norm,parlocnorm;
01759 double *gv=NULL;
01760 MPI_Status stat;
01761
01762 if (myrank==0){
01763 gv = new double [ndofcp];
01764 nullv (gv,ndofcp);
01765 }
01766
01767
01768 parlocnorm =0;
01769 for (i=0;i<nidof;i++){
01770 parlocnorm+=lv[i]*lv[i];
01771 }
01772
01773
01774 if (myrank==0){
01775 norm=parlocnorm;
01776 for (i=1;i<nproc;i++){
01777 MPI_Recv (&aux,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01778 norm+=aux;
01779 }
01780 }
01781 else{
01782 MPI_Send(&parlocnorm,1,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01783 }
01784 MPI_Barrier (MPI_COMM_WORLD);
01785
01786
01787 gather_bound_vect (lv,gv,domproc);
01788
01789
01790 if (myrank==0){
01791 for (i=0;i<ndofcp;i++){
01792 norm+=gv[i]*gv[i];
01793
01794 }
01795 }
01796
01797
01798 scatter_bound_vect (lv,gv,domproc);
01799
01800
01801 if (myrank==0){
01802 for (i=1;i<nproc;i++){
01803 MPI_Send(&norm,1,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01804 }
01805 }
01806 else{
01807 MPI_Recv (&norm,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01808 }
01809 MPI_Barrier (MPI_COMM_WORLD);
01810
01811 if (myrank==0)
01812 delete [] gv;
01813
01814 return norm;
01815 }
01816