00001 #include "feti1.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <time.h>
00005
00006 feti1::feti1(int np,int mr,long nd)
00007 {
00008 nproc=np; myrank=mr; ndom=nd;
00009
00010 ndof=0; maxlggl=0; nclggl=0; maxnrbm=0;
00011 enrbm=0; lithr=0.0; nrbm=0;
00012 nicg=0; anicg=0; errcg=0.0; aerrcg=0.0;
00013 zero=1.0e-15;
00014
00015 ndofcp=0; hsize=0;
00016
00017 lcn=NULL;
00018
00019 lggl=NULL;
00020 glcor=NULL;
00021
00022 nrbmdom=NULL;
00023 rbmadr=NULL;
00024
00025
00026
00027
00028 edofs = NULL;
00029
00030 ncnd = NULL;
00031
00032
00033
00034
00035 ndofncn = NULL;
00036 cncn = NULL;
00037
00038
00039 wscalmat = NULL;
00040 lwscalmat = NULL;
00041
00042 smsky = NULL;
00043
00044 fetiprecond = -1;
00045 }
00046
00047 feti1::~feti1()
00048 {
00049 long i;
00050
00051 delete [] lcn;
00052
00053
00054 for (i=0;i<nproc;i++){
00055 delete [] glcor[i];
00056 }
00057 delete [] glcor;
00058
00059 delete [] lggl;
00060 delete [] nrbmdom;
00061 delete [] rbmadr;
00062
00063
00064
00065
00066 delete [] edofs;
00067
00068 if (myrank==0){
00069 delete [] ncnd;
00070 }
00071
00072
00073
00074
00075 }
00076
00077
00078
00079
00080 void feti1::initiate (selnodes *selnodfeti,FILE *out)
00081 {
00082 long i;
00083
00084 maxlggl = selnodfeti->maxndof;
00085
00086 nclggl = selnodfeti->ndof;
00087
00088
00089 maxncdofd = selnodfeti->maxndof;
00090
00091
00092 ncdof = selnodfeti->ndof;
00093
00094
00095 lcn = selnodfeti->ldof;
00096
00097
00098 edofs = selnodfeti->ldof;
00099
00100 fprintf (out,"\n\n\n osel \n");
00101 for (long i=0;i<ncdof;i++){
00102 fprintf (out,"\n edofs %ld",edofs[i]);
00103 }
00104 fprintf (out,"\n\n\n osel \n");
00105
00106
00107 if (myrank==0){
00108
00109 ndofcp =selnodfeti->tndof;
00110
00111 ncdofd = selnodfeti->ndofdom;
00112
00113
00114
00115 ccn = selnodfeti->cndom;
00116
00117 wscalmat = new double [ndofcp];
00118 for (i = 0; i < ndofcp; i++){
00119 wscalmat[i] = double(selnodfeti->dofmultip[i]);
00120 }
00121
00122 }
00123
00124 }
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134 void feti1::coarse_dofs (partop *ptop,FILE *out)
00135 {
00136 long i,j,k,m,n;
00137 long tnbn;
00138 long *bnmultip,**nbdofnd,**llnbn,**ldn,***lbcn;
00139
00140 if (myrank==0){
00141
00142
00143 tnbn = ptop->tnbn;
00144
00145 bnmultip = ptop->bnmultip;
00146
00147 nbdofnd = ptop->nbdofnd;
00148
00149 llnbn = ptop->llnbn;
00150
00151 ldn = ptop->ldn;
00152
00153 lbcn = ptop->lbcn;
00154
00155
00156
00157
00158 ndofncn = new long* [tnbn];
00159 for (i=0;i<tnbn;i++){
00160 ndofncn[i] = new long [bnmultip[i]];
00161 }
00162
00163 for (i=0;i<tnbn;i++){
00164 for (j=0;j<bnmultip[i];j++){
00165 k=llnbn[i][j];
00166 m=ldn[i][j];
00167 ndofncn[i][j]=nbdofnd[m][k];
00168 }
00169 }
00170
00171
00172
00173
00174 cncn = new long** [tnbn];
00175 for (i=0;i<tnbn;i++){
00176 cncn[i] = new long* [bnmultip[i]];
00177 for (j=0;j<bnmultip[i];j++){
00178 cncn[i][j] = new long [ndofncn[i][j]];
00179 }
00180 }
00181
00182
00183 for (i=0;i<tnbn;i++){
00184 for (j=0;j<bnmultip[i];j++){
00185 m=llnbn[i][j];
00186 n=ldn[i][j];
00187 for (k=0;k<ndofncn[i][j];k++){
00188 cncn[i][j][k]=lbcn[n][m][k];
00189 }
00190 }
00191 }
00192
00193
00194 ndofcp=1;
00195 for (i=0;i<tnbn;i++){
00196 for (j=0;j<bnmultip[i]-1;j++){
00197 for (k=0;k<ndofncn[i][j];k++){
00198 if (cncn[i][j][k]>0){
00199 cncn[i][j][k]=ndofcp;
00200 ndofcp++;
00201 }
00202 }
00203 }
00204 }
00205 ndofcp--;
00206
00207 fprintf (out,"\n\n number of unknowns (DOFs) in the coarse problem is %ld",ndofcp);
00208 fprintf (stdout,"\n\n number of unknowns (DOFs) in the coarse problem is %ld",ndofcp);
00209 }
00210 }
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 void feti1::number_contributing_nodes_dofs (gtopology *top,partop *ptop,long *domproc,FILE *out)
00232 {
00233 long i,j,k,l,m,ii,jj,ndofn,buffsize;
00234 long *buff;
00235 long *lnbn,*bnmultip,*nbnd,**ldn,**cnbn;
00236 MPI_Status stat;
00237
00238 buffsize=3;
00239 buff = new long [buffsize];
00240
00241
00242 lnbn = ptop->lnbn;
00243
00244 if (myrank==0){
00245
00246 bnmultip = ptop->bnmultip;
00247
00248 nbnd = ptop->nbnd;
00249
00250 ldn = ptop->ldn;
00251
00252 cnbn = ptop->cnbn;
00253
00254
00255 ncnd = new long [nproc];
00256
00257 ncdofd = new long [nproc];
00258
00259
00260 maxncnd = 0;
00261
00262 maxncdofd = 0;
00263
00264
00265 for (i=0;i<nproc;i++){
00266 ncnd[i]=0;
00267 ncdofd[i]=0;
00268 for (j=0;j<nbnd[i];j++){
00269 m=cnbn[i][j];
00270 for (k=0;k<bnmultip[m];k++){
00271 if (ldn[m][k]==i){
00272 if (k==0){
00273 for (ii=0;ii<bnmultip[m]-1;ii++){
00274 ncnd[i]++;
00275 for (jj=0;jj<ndofncn[m][ii];jj++){
00276 if (cncn[m][ii][jj]>0)
00277 ncdofd[i]++;
00278 }
00279 }
00280 }
00281 else{
00282 ncnd[i]++;
00283 for (jj=0;jj<ndofncn[m][k-1];jj++){
00284 if (cncn[m][k-1][jj]>0) ncdofd[i]++;
00285 }
00286
00287 }
00288 break;
00289 }
00290 }
00291 }
00292
00293 if (maxncnd<ncnd[i])
00294 maxncnd=ncnd[i];
00295 if (maxncdofd<ncdofd[i])
00296 maxncdofd=ncdofd[i];
00297 }
00298
00299 for (i=1;i<nproc;i++){
00300 j=domproc[i];
00301 buff[0]=maxncnd;
00302 buff[1]=ncnd[j];
00303 buff[2]=maxncdofd;
00304 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00305 }
00306
00307 j=domproc[0];
00308 buff[0]=maxncnd;
00309 buff[1]=ncnd[j];
00310 buff[2]=maxncdofd;
00311 }
00312 else{
00313 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00314 }
00315 MPI_Barrier (MPI_COMM_WORLD);
00316
00317 maxncnd=buff[0];
00318 ncn=buff[1];
00319 maxncdofd=buff[2];
00320
00321 delete [] buff;
00322
00323 fprintf (out,"\n\n KONTROLA 1. maxncnd je %ld ncn je %ld maxncdofd je %ld",maxncnd,ncn,maxncdofd);
00324
00325 if (myrank==0){
00326 fprintf (out,"\n\n kontrola poctu prispivajicich uzlu ncnd\n");
00327 for (i=0;i<nproc;i++){
00328 fprintf (out," %ld",ncnd[i]);
00329 }
00330 fprintf (out,"\n");
00331 for (i=0;i<nproc;i++){
00332 fprintf (out," %ld",ncdofd[i]);
00333 }
00334 fprintf (out,"\n");
00335 }
00336
00337 }
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350 void feti1::contributing_nodes_dofs (gtopology *top,partop *ptop,long *domproc,FILE *out)
00351 {
00352 long i,j,k,l,m,ii,jj,ndofn,buffsize;
00353 long tnbn;
00354 long *buff;
00355 long *lnbn,*bnmultip,**ldn,*nbnd,**llnbn,**cnbn;
00356 MPI_Status stat;
00357
00358
00359 lnbn = ptop->lnbn;
00360
00361 if (myrank==0){
00362
00363 tnbn = ptop->tnbn;
00364
00365 bnmultip = ptop->bnmultip;
00366
00367 ldn = ptop->ldn;
00368
00369 nbnd = ptop->nbnd;
00370
00371 llnbn = ptop->llnbn;
00372
00373 cnbn = ptop->cnbn;
00374 }
00375
00376
00377
00378 buffsize = maxncnd;
00379 buff = new long [buffsize];
00380
00381 if (myrank==0){
00382
00383 for (i=1;i<nproc;i++){
00384
00385 for (j=0;j<buffsize;j++){
00386 buff[j]=0;
00387 }
00388
00389 m=0;
00390 for (j=0;j<tnbn;j++){
00391 for (k=0;k<bnmultip[j];k++){
00392 if (ldn[j][k]==domproc[i]){
00393 if (k==0){
00394 for (l=0;l<bnmultip[j]-1;l++){
00395 buff[m]=llnbn[j][k];
00396 buff[m]=0-buff[m]-1;
00397 m++;
00398 }
00399 }
00400 else{
00401 buff[m]=llnbn[j][k];
00402 m++;
00403 }
00404 break;
00405 }
00406 }
00407 }
00408
00409 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00410 }
00411
00412 i=0;
00413 for (j=0;j<buffsize;j++){
00414 buff[j]=0;
00415 }
00416
00417 m=0;
00418 for (j=0;j<tnbn;j++){
00419 for (k=0;k<bnmultip[j];k++){
00420 if (ldn[j][k]==domproc[i]){
00421 if (k==0){
00422 for (l=0;l<bnmultip[j]-1;l++){
00423 buff[m]=llnbn[j][k];
00424 buff[m]=0-buff[m]-1;
00425 m++;
00426 }
00427 }
00428 else{
00429 buff[m]=llnbn[j][k];
00430 m++;
00431 }
00432 break;
00433 }
00434 }
00435 }
00436
00437
00438 }
00439 else{
00440 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00441 }
00442 MPI_Barrier (MPI_COMM_WORLD);
00443
00444 fprintf (out,"\n\n KONTROLA 2. kontrola hranicnich uzlu \n");
00445 for (i=0;i<ncn;i++){
00446 if (buff[i]<0){
00447 j=0-buff[i]-1;
00448 }
00449 else
00450 j=buff[i];
00451
00452 fprintf (out," %ld %ld",buff[i],lnbn[j]);
00453 }
00454 fprintf (out,"\n");
00455
00456
00457 ncdof=0;
00458 for (i=0;i<ncn;i++){
00459 if (buff[i]<0){
00460 j=0-buff[i]-1;
00461 k=lnbn[j];
00462 ndofn=top->give_ndofn (k);
00463 for (j=0;j<ndofn;j++){
00464 l=top->give_dof (k,j);
00465 if (l>0) ncdof++;
00466 }
00467 }
00468 else{
00469 k=lnbn[buff[i]];
00470 ndofn=top->give_ndofn (k);
00471 for (j=0;j<ndofn;j++){
00472 l=top->give_dof (k,j);
00473 if (l>0) ncdof++;
00474 }
00475 }
00476 }
00477
00478 fprintf (out,"\n\n KONTROLA 3. pocet stupnu volnosti prispivajicich do hrubeho problemu %ld",ncdof);
00479
00480
00481
00482
00483
00484
00485 edofs = new long [ncdof];
00486 m=0;
00487 for (i=0;i<ncn;i++){
00488 if (buff[i]<0){
00489 j=0-buff[i]-1;
00490 k=lnbn[j];
00491 ndofn=top->give_ndofn (k);
00492 for (j=0;j<ndofn;j++){
00493 l=top->give_dof (k,j);
00494 if (l>0){
00495 edofs[m]=0-l;
00496 m++;
00497 }
00498 }
00499 }
00500 else{
00501 k=lnbn[buff[i]];
00502 ndofn=top->give_ndofn (k);
00503 for (j=0;j<ndofn;j++){
00504 l=top->give_dof (k,j);
00505 if (l>0){
00506 edofs[m]=l;
00507 m++;
00508 }
00509 }
00510 }
00511 }
00512
00513 delete [] buff;
00514
00515
00516 fprintf (out,"\n\n\n KONTROLA 4. extract code numbers edofs\n");
00517 for (i=0;i<ncdof;i++){
00518 fprintf (out,"\n%ld",edofs[i]);
00519 }
00520 fprintf (out,"\n");
00521
00522
00523
00524
00525 if (myrank==0){
00526
00527
00528
00529 ccn = new long* [nproc];
00530 for (i=0;i<nproc;i++){
00531 ccn[i] = new long [ncdofd[i]];
00532 }
00533
00534 for (i=0;i<nproc;i++){
00535 l=0;
00536 for (j=0;j<nbnd[i];j++){
00537 m=cnbn[i][j];
00538 for (k=0;k<bnmultip[m];k++){
00539 if (ldn[m][k]==i){
00540 if (k==0){
00541 for (ii=0;ii<bnmultip[m]-1;ii++){
00542 for (jj=0;jj<ndofncn[m][ii];jj++){
00543 if (cncn[m][ii][jj]>0){
00544
00545 ccn[i][l]=cncn[m][ii][jj];
00546 l++;
00547 }
00548 }
00549 }
00550 }
00551 else{
00552 for (jj=0;jj<ndofncn[m][k-1];jj++){
00553 if (cncn[m][k-1][jj]>0){
00554 ccn[i][l]=cncn[m][k-1][jj];
00555 l++;
00556 }
00557 }
00558
00559 }
00560 break;
00561 }
00562 }
00563 }
00564
00565 }
00566
00567 fprintf (out,"\n\n\n kontrola coarse code numbers ccn\n");
00568 for (i=0;i<nproc;i++){
00569 fprintf (out,"\n subdomain %ld ncdofd %ld",i,ncdofd[i]);
00570 for (j=0;j<ncdofd[i];j++){
00571 fprintf (out,"\n%ld",ccn[i][j]);
00572 }
00573 }
00574 fprintf (out,"\n");
00575
00576 }
00577
00578 }
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590 void feti1::coarse_problem_ordering (gtopology *top,partop *ptop,long *domproc,FILE *out)
00591 {
00592
00593
00594 coarse_dofs (ptop,out);
00595
00596
00597 number_contributing_nodes_dofs (top,ptop,domproc,out);
00598
00599
00600
00601 contributing_nodes_dofs (top,ptop,domproc,out);
00602 }
00603
00604
00605
00606
00607
00608
00609
00610
00611 void feti1::subdomain_ordering (gtopology *top,long *ltg,FILE *out)
00612 {
00613 if(fetiprecond == dirichlet){
00614
00615 long i,j,k,ndofn;
00616
00617
00618 ndof=1;
00619 for (i=0;i<top->nn;i++){
00620 if (ltg[i]==-1){
00621
00622 ndofn=top->give_ndofn (i);
00623 for (j=0;j<ndofn;j++){
00624 k=top->give_dof(i,j);
00625 if (k<0) continue;
00626 if (k==0) continue;
00627 if (k>0){
00628 top->save_dof (i,j,ndof); ndof++;
00629 }
00630 }
00631 }
00632 }
00633
00634 for (i=0;i<top->nn;i++){
00635 if (ltg[i]>-1){
00636
00637 ndofn=top->give_ndofn (i);
00638 for (j=0;j<ndofn;j++){
00639 k=top->give_dof (i,j);
00640 if (k<0) continue;
00641 if (k==0) continue;
00642 if (k>0){
00643 top->save_dof (i,j,ndof); ndof++;
00644 }
00645 }
00646 }
00647 }
00648 ndof--;
00649
00650 }
00651 else{
00652
00653 ndof = top->codenum_generation (out);
00654 }
00655 }
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673 void feti1::local_buff (double *lv,double *buff)
00674 {
00675 long i,j;
00676
00677 for (i=0;i<ncdof;i++){
00678 j=edofs[i];
00679 if (j<0){
00680 j=0-j-1;
00681 buff[i]=0.0-lv[j];
00682 }
00683 else{
00684 j=j-1;
00685 buff[i]=lv[j];
00686 }
00687 }
00688 }
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699 void feti1::buff_coarse (double *cv,double *buff,long nd)
00700 {
00701 long i,j;
00702
00703 for (i=0;i<ncdofd[nd];i++){
00704 j=ccn[nd][i]-1;
00705 cv[j]+=buff[i];
00706 }
00707 }
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718 void feti1::coarse_buff (double *cv,double *buff,long nd)
00719 {
00720 long i,j;
00721
00722 for (i=0;i<ncdofd[nd];i++){
00723 j=ccn[nd][i]-1;
00724 buff[i]=cv[j];
00725 }
00726 }
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736 void feti1::buff_local (double *lv,double *buff)
00737 {
00738 long i,j;
00739
00740 for (i=0;i<ncdof;i++){
00741 j=edofs[i];
00742 if (j<0){
00743 j=0-j-1;
00744 lv[j]-=buff[i];
00745 }
00746 else{
00747 j=j-1;
00748 lv[j]+=buff[i];
00749 }
00750 }
00751 }
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
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
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
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
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495 void feti1::hmatrixsize (double *rbm,long *domproc,FILE *out)
01496 {
01497 long i,j,k;
01498 MPI_Status stat;
01499
01500 if (myrank==0){
01501
01502 nrbmdom = new long [nproc];
01503
01504 rbmadr = new long [nproc+1];
01505 for (i=0;i<nproc+1;i++){
01506 rbmadr[i]=0;
01507 }
01508
01509 j=domproc[0];
01510 maxnrbm=nrbm; nrbmdom[j]=nrbm;
01511 for (i=1;i<nproc;i++){
01512 MPI_Recv(&k,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01513 if (maxnrbm<k) maxnrbm=k;
01514 j=domproc[stat.MPI_TAG];
01515 nrbmdom[j]=k;
01516 }
01517
01518
01519 rbmadr[0]=0;
01520 for (i=0;i<nproc;i++){
01521 rbmadr[i+1]=rbmadr[i]+nrbmdom[i];
01522 }
01523
01524 for (i=1;i<nproc;i++){
01525 MPI_Send(&maxnrbm,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01526 }
01527
01528
01529 hsize=rbmadr[nproc];
01530 }
01531 else{
01532 MPI_Send(&nrbm,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01533 MPI_Recv(&maxnrbm,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01534 }
01535
01536 }
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550 void feti1::hmatrix (double *h,double *rbm,long *domproc)
01551 {
01552 long i,j,k,l,buffsize;
01553 double *buff,*av;
01554 MPI_Status stat;
01555
01556 buffsize = maxnrbm*maxncdofd;
01557 buff = new double [buffsize];
01558 for (i=0;i<buffsize;i++){
01559 buff[i]=0.0;
01560 }
01561
01562
01563 for (i=0;i<nrbm;i++){
01564 local_buff (rbm+i*ndof,buff+i*maxncdofd);
01565 }
01566
01567
01568 if (myrank==0){
01569
01570 av = new double [ndofcp];
01571
01572
01573 j=domproc[0];
01574 for (k=0;k<nrbmdom[j];k++){
01575 nullv (av,ndofcp);
01576 buff_coarse (av,buff+k*maxncdofd,j);
01577 for (l=0;l<ndofcp;l++){
01578 h[l*hsize+rbmadr[j]+k]=0.0-av[l];
01579 }
01580 }
01581
01582
01583 for (i=1;i<nproc;i++){
01584 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01585
01586
01587 j=domproc[stat.MPI_TAG];
01588 for (k=0;k<nrbmdom[j];k++){
01589 nullv (av,ndofcp);
01590 buff_coarse (av,buff+k*maxncdofd,j);
01591 for (l=0;l<ndofcp;l++){
01592 h[l*hsize+rbmadr[j]+k]=0.0-av[l];
01593 }
01594 }
01595
01596 }
01597
01598 delete [] av;
01599 }
01600 else{
01601 MPI_Send(buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01602 }
01603 MPI_Barrier (MPI_COMM_WORLD);
01604
01605 delete [] buff;
01606 }
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621 void feti1::qvector (double *q,double *rbm,double *f,long *domproc)
01622 {
01623 long i,j,k,l;
01624 double *buff;
01625 MPI_Status stat;
01626
01627 buff = new double [maxnrbm];
01628
01629 for (i=0;i<nrbm;i++){
01630 buff[i]=0.0-ss(rbm+i*ndof,f,ndof);
01631 }
01632
01633 if (myrank==0){
01634
01635 k=domproc[0];
01636 l=rbmadr[k];
01637 for (j=0;j<nrbmdom[k];j++){
01638 q[l]=buff[j]; l++;
01639 }
01640
01641 for (i=1;i<nproc;i++){
01642 MPI_Recv(buff,maxnrbm,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01643 k=domproc[stat.MPI_TAG];
01644 l=rbmadr[k];
01645 for (j=0;j<nrbmdom[k];j++){
01646 q[l]=buff[j]; l++;
01647 }
01648 }
01649
01650 }
01651 else{
01652 MPI_Send(buff,maxnrbm,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01653 }
01654
01655 delete [] buff;
01656 }
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673 void feti1::feti_projection (double *v,double *h,double *h1)
01674 {
01675 long i;
01676 double *p,*q,*w;
01677
01678 p = new double [hsize];
01679 q = new double [hsize];
01680 w = new double [ndofcp];
01681
01682
01683
01684 mtxv (h,v,p,ndofcp,hsize);
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696 mxv (h1,p,q,hsize,hsize);
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708 mxv (h,q,w,ndofcp,hsize);
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719 for (i=0;i<ndofcp;i++){
01720 v[i]-=w[i];
01721 }
01722
01723 delete [] w; delete [] q; delete [] p;
01724 }
01725
01726
01727
01728
01729 void feti1::subdomain_matrix (gmatrix *gm,long *domproc,FILE *out)
01730 {
01731 long i,j,k;
01732 long *aux;
01733 MPI_Status stat;
01734 long buffsize;
01735 double *buff,*auxscal;
01736
01737 fprintf(out,"\n\n\nkontrola subdomain_matrix\n\n\n\n kontrola ndof %ld\n\n\n",ndof);
01738 aux = new long [ndof];
01739 for(i = 0; i < ndof; i++){
01740 aux[i]=0;
01741 }
01742
01743 fprintf(out,"kontrola ncdof %ld\n\n\n",ncdof);
01744
01745 for (i = 0; i < ncdof; i++){
01746 k = edofs[i];
01747
01748 if (k<0)
01749 k=0-k-1;
01750 else
01751 k--;
01752 if (k>-1)
01753 aux[k]++;
01754 }
01755
01756
01757 fprintf (out,"aux domene\n");
01758 for (i = 0; i < ndof; i++){
01759 fprintf (out," %ld\n",aux[i]);
01760 }
01761
01762 ndofprec = 0;
01763 for (i = 0;i < ndof; i++){
01764 if (aux[i] > 0)
01765 ndofprec++;
01766 }
01767 fprintf (out,"ndofprec %ld\n",ndofprec);
01768
01769
01770 cnprec = new long [ndofprec];
01771
01772
01773 k=0;
01774 for (i = 0; i < ndof;i++){
01775 if (aux[i]>0){
01776 cnprec[k]=i;
01777
01778 k++;
01779 }
01780 }
01781 delete [] aux;
01782
01783
01784
01785
01786
01787
01788
01789
01790 fprintf (out,"Kontrola cnprec\n");
01791 for (i = 0; i < ndofprec; i++){
01792 fprintf (out," cnprec[%ld] = %ld\n",i,cnprec[i]);
01793 }
01794
01795 switch (fetiprecond){
01796 case lumped:{
01797 smsky = new skyline;
01798 gm->sky->select_submatrix(smsky,ndofprec,cnprec,out);
01799 break;
01800 }
01801 case dirichlet:{
01802
01803
01804
01805 smsky = new skyline;
01806 smsky->n=gm->sky->n;
01807
01808
01809 smsky->adr = new long [smsky->n+1];
01810 for(i = 0; i < smsky->n+1;i++){
01811 smsky->adr[i]=gm->sky->adr[i];
01812 }
01813
01814
01815
01816
01817
01818
01819 smsky->adr[smsky->n]= gm->sky->negm;
01820 smsky->negm = gm->sky->negm;
01821
01822 smsky->a = new double [smsky->negm];
01823
01824 for(i = 0; i < smsky->negm;i++){
01825 smsky->a[i]=gm->sky->a[i];
01826 }
01827
01828
01829
01830
01831
01832
01833 smdm = new densemat;
01834 double *x,*y,*av;
01835 x = new double[ndof];
01836 y = new double[ndof];
01837 smdm->a = new double [ndofprec*ndofprec];
01838 av = new double [ndofprec];
01839 smdm->negm=ndofprec*ndofprec;
01840 smdm->n=ndofprec;
01841
01842 smsky->ldlkon_sky (smdm->a,av,x,y,ndofprec,1);
01843
01844 smsky->~skyline ();
01845
01846 for(i = 0; i < smdm->n*smdm->n; i++){
01847 fprintf(out,"%le\n",smdm->a[i]);
01848 }
01849
01850 delete []x;
01851 delete []y;
01852 delete [] av;
01853
01854 break;
01855 }
01856 }
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938 }
01939
01940
01941
01942
01943
01944 void feti1::scaling(double *invect,double *outvect,long n,FILE *out)
01945 {
01946 double *aux;
01947 long i;
01948 aux = new double[n];
01949
01950 for(i = 0; i < n; i++){
01951 aux[i] = invect[i]/wscalmat[i];
01952 }
01953
01954
01955 fprintf (out,"\n\n\n kontrola skalovani \n");
01956 for(i = 0; i < n; i++){
01957 fprintf (out,"%ld pred %le po %le wscalmat %lf\n",i+1,invect[i],aux[i],wscalmat[i]);
01958 }
01959
01960
01961 for(i = 0; i < n; i++){
01962 outvect[i] = aux[i];
01963 }
01964
01965 delete []aux;
01966 }
01967
01968 void feti1::locscaling (double *invect,double *outvect,FILE *out)
01969 {
01970 long i;
01971 double *aux;
01972
01973 aux = new double[ndofprec];
01974
01975 for(i = 0; i < ndofprec; i++){
01976 aux[i] = invect[i]/lwscalmat[i];
01977 }
01978
01979
01980 fprintf (out,"\n\n\n kontrola skalovani \n");
01981 for(i = 0; i < ndofprec; i++){
01982 fprintf (out,"%ld pred %le po %le lwscalmat %lf\n",i+1,invect[i],aux[i],lwscalmat[i]);
01983 }
01984
01985
01986 for(i = 0; i < ndofprec; i++){
01987 outvect[i] = aux[i];
01988 }
01989
01990 delete []aux;
01991 }
01992
01993
01994
01995
01996 void feti1::lumpedprec (gmatrix *gm,double *dd,double *pp,FILE *out)
01997 {
01998
01999 long i;
02000 double *d,*p;
02001
02002 d = new double [ndofprec];
02003 p = new double [ndofprec];
02004
02005 for (i = 0; i < ndofprec; i++){
02006 d[i]=0.0;
02007 p[i]=0.0;
02008 }
02009
02010 for (i=0;i<ndofprec;i++){
02011 d[i]=dd[cnprec[i]];
02012 }
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036 smsky->mxv_sky (d,p);
02037
02038
02039
02040
02041 for (i = 0;i < ndofprec; i++){
02042 pp[cnprec[i]]=p[i];
02043 }
02044
02045 delete [] p;
02046 delete [] d;
02047
02048 }
02049
02050
02051 void feti1::dirichletprec (double *dd,double *pp,FILE *out)
02052 {
02053 long i;
02054 double *d,*p;
02055
02056 d = new double [ndofprec];
02057 p = new double [ndofprec];
02058
02059 for (i=0;i<ndofprec;i++){
02060 d[i]=0.0;
02061 p[i]=0.0;
02062 }
02063
02064 for (i=0;i<ndofprec;i++){
02065 d[i]=dd[cnprec[i]];
02066 }
02067
02068
02069
02070
02071
02072
02073
02074 smdm->mxv_dm (d,p);
02075
02076
02077
02078 for (i=0;i<ndofprec;i++){
02079 pp[cnprec[i]]=p[i];
02080 }
02081
02082 delete [] p;
02083 delete [] d;
02084 }
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098
02099
02100
02101
02102 void feti1::mpcg (gtopology *top,gmatrix *gm,long *domproc,
02103 double *w,double *rhs,double *q,double *h,double *h1,long *rbmi,FILE *out)
02104 {
02105 long i,j,k,buffsize;
02106 long l;
02107 double nom,denom,alpha,beta;
02108 double *d,*dd,*g,*p,*pp,*hq,*buff;
02109 MPI_Status stat;
02110
02111
02112 dd = new double [ndof];
02113 pp = new double [ndof];
02114
02115
02116 buffsize=maxncdofd+1;
02117 buff = new double [buffsize];
02118 buff[maxncdofd]=0.0;
02119
02120 if (myrank==0){
02121
02122
02123 d = new double [ndofcp];
02124
02125 g = new double [ndofcp];
02126
02127 p = new double [ndofcp];
02128
02129
02130
02131
02132
02133 hq = new double [hsize];
02134
02135 mxv (h1,q,hq,hsize,hsize);
02136
02137 mxv (h,hq,w,ndofcp,hsize);
02138 delete [] hq;
02139
02140 for (j=1;j<nproc;j++){
02141 k=domproc[j];
02142 coarse_buff (w,buff,k);
02143 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02144 }
02145
02146 k=domproc[0];
02147 coarse_buff (w,buff,k);
02148 nullv (dd,ndof);
02149 buff_local (dd,buff);
02150
02151 subv (dd,rhs,ndof);
02152
02153 fprintf (out,"\n");
02154 for (i=0;i<ndof;i++){
02155 fprintf (out,"\n zoufalec %4ld %20.15e",i,dd[i]);
02156 }
02157
02158
02159 gm->ldl_feti (pp,dd,nrbm,rbmi,zero);
02160
02161 fprintf (out,"\n");
02162 for (i=0;i<nrbm;i++){
02163 fprintf (out,"\n zoufalec se %4ld",rbmi[i]);
02164 }
02165 fprintf (out,"\n");
02166 for (i=0;i<ndof;i++){
02167 fprintf (out,"\n zoufalec %4ld %20.15e",i,pp[i]);
02168 }
02169
02170
02171 local_buff (pp,buff);
02172
02173
02174 nullv (g,ndofcp);
02175 k=domproc[0];
02176 buff_coarse (g,buff,k);
02177
02178 for (j=1;j<nproc;j++){
02179 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02180 k=domproc[stat.MPI_TAG];
02181 buff_coarse (g,buff,k);
02182 }
02183
02184
02185 fprintf (out,"\n\n\n");
02186 for (i=0;i<ndofcp;i++){
02187 fprintf (out,"\n zoufalec %4ld %20.15e",i,g[i]);
02188 }
02189
02190
02191
02192 feti_projection (g,h,h1);
02193
02194
02195 for (i=0;i<ndofcp;i++){
02196 d[i]=0.0-g[i];
02197 }
02198
02199
02200 nom = ss (g,g,ndofcp);
02201
02202
02203 }
02204 else{
02205
02206
02207 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02208 nullv (dd,ndof);
02209
02210 buff_local (dd,buff);
02211
02212 subv (dd,rhs,ndof);
02213
02214 fprintf (out,"\n");
02215 for (i=0;i<ndof;i++){
02216 fprintf (out,"\n zoufalec %4ld %20.15e",i,dd[i]);
02217 }
02218
02219
02220 gm->ldl_feti (pp,dd,nrbm,rbmi,zero);
02221
02222 fprintf (out,"\n");
02223 for (i=0;i<nrbm;i++){
02224 fprintf (out,"\n zoufalec se %4ld",rbmi[i]);
02225 }
02226 fprintf (out,"\n");
02227 for (i=0;i<ndof;i++){
02228 fprintf (out,"\n zoufalec %4ld %20.15e",i,pp[i]);
02229 }
02230
02231
02232 local_buff (pp,buff);
02233 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
02234
02235 }
02236 MPI_Barrier (MPI_COMM_WORLD);
02237
02238
02239
02240 if (myrank==0){
02241
02242
02243
02244
02245 for (i=0;i<nicg;i++){
02246
02247
02248
02249
02250 for (j=1;j<nproc;j++){
02251
02252 k=domproc[j];
02253 coarse_buff (d,buff,k);
02254 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02255 }
02256
02257
02258
02259
02260 k=domproc[0];
02261 coarse_buff (d,buff,k);
02262 nullv (dd,ndof);
02263 buff_local (dd,buff);
02264
02265 gm->ldl_feti (pp,dd,nrbm,rbmi,zero);
02266
02267
02268 local_buff (pp,buff);
02269
02270 nullv (p,ndofcp);
02271 k=domproc[0];
02272 buff_coarse (p,buff,k);
02273
02274 for (j=1;j<nproc;j++){
02275 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02276 k=domproc[stat.MPI_TAG];
02277 buff_coarse (p,buff,k);
02278 }
02279
02280
02281
02282 denom = ss (d,p,ndofcp);
02283
02284
02285
02286 fprintf (out,"\n\n kontrola citatele a jmenovatele pred alpha %e / %e",nom,denom);
02287
02288
02289 if (fabs(denom)<zero){
02290 fprintf (stderr,"\n\n zero denominator in modified conjugate gradient method (%s, line %d).\n",__FILE__,__LINE__);
02291
02292
02293 buff[maxncdofd]=1.0;
02294 for (j=1;j<nproc;j++){
02295 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02296 }
02297
02298 break;
02299 }
02300
02301
02302
02303
02304
02305
02306 alpha = nom/denom;
02307
02308
02309
02310
02311 for (j=0;j<ndofcp;j++){
02312 w[j]+=alpha*d[j];
02313 g[j]+=alpha*p[j];
02314 }
02315
02316 feti_projection (g,h,h1);
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328 nullv (p,ndofcp);
02329 switch (fetiprecond){
02330
02331 case nofetiprecond:{
02332
02333 copyv(g,p,ndofcp);
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368 break;
02369 }
02370 case lumped:{
02371
02372 scaling(g,g,ndofcp,out);
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383 for (j=1;j<nproc;j++){
02384
02385 k=domproc[j];
02386 coarse_buff (g,buff,k);
02387 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02388
02389
02390
02391
02392
02393
02394
02395
02396 }
02397
02398
02399 k=domproc[0];
02400 coarse_buff (g,buff,k);
02401 nullv (dd,ndof);
02402 buff_local (dd,buff);
02403
02404
02405
02406
02407
02408
02409 nullv (pp,ndof);
02410 lumpedprec (gm,dd,pp,out);
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424 local_buff (pp,buff);
02425
02426 nullv (p,ndofcp);
02427 k=domproc[0];
02428 buff_coarse (p,buff,k);
02429
02430 for (j=1;j<nproc;j++){
02431 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02432 k=domproc[stat.MPI_TAG];
02433 buff_coarse (p,buff,k);
02434
02435
02436
02437
02438 }
02439
02440
02441
02442
02443
02444
02445 scaling(p,p,ndofcp,out);
02446
02447 break;
02448 }
02449 case dirichlet:{
02450
02451 scaling(g,g,ndofcp,out);
02452
02453
02454
02455
02456
02457
02458
02459
02460 for (j=1;j<nproc;j++){
02461
02462 k=domproc[j];
02463 coarse_buff (g,buff,k);
02464 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02465
02466
02467
02468
02469
02470
02471
02472
02473 }
02474
02475
02476 k=domproc[0];
02477 coarse_buff (g,buff,k);
02478 nullv (dd,ndof);
02479 buff_local (dd,buff);
02480
02481
02482
02483
02484
02485
02486 nullv (pp,ndof);
02487 dirichletprec (dd,pp,out);
02488
02489
02490
02491
02492
02493
02494
02495 for (l = 0; l < ndof; l++){
02496 fprintf(out,"pp[%ld] %le\n",l,pp[l]);
02497 }
02498
02499
02500
02501 local_buff (pp,buff);
02502
02503 nullv (p,ndofcp);
02504 k=domproc[0];
02505 buff_coarse (p,buff,k);
02506
02507 for (j=1;j<nproc;j++){
02508 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02509 k=domproc[stat.MPI_TAG];
02510 buff_coarse (p,buff,k);
02511
02512
02513
02514
02515 }
02516
02517 for (l = 0; l < ndofcp; l++){
02518 fprintf(out,"p[%ld] %le\n",l,p[l]);
02519 }
02520
02521
02522 scaling(p,p,ndofcp,out);
02523 break;
02524 }
02525 }
02526
02527
02528
02529
02530 feti_projection (p,h,h1);
02531
02532
02533
02534
02535
02536
02537
02538
02539 denom = nom;
02540 if (fabs(denom)<zero){
02541 fprintf (stderr,"\n\n zero denominator of beta in function (%s, line %d).\n",__FILE__,__LINE__);
02542
02543 if (i!=nicg-1){
02544 buff[maxncdofd]=1.0;
02545 for (j=1;j<nproc;j++){
02546 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02547 }
02548 }
02549 break;
02550 }
02551
02552
02553 nom=ss(p,g,ndofcp);
02554
02555 fprintf (stdout,"\n iteration %4ld norm g %e",i,sqrt(nom));
02556 fprintf (out,"\n iteration %4ld norm g %e",i,sqrt(nom));
02557
02558 fprintf (out,"\n kontrola citatele a jmenovatele pred betou %e / %e",nom,denom);
02559
02560
02561 if (sqrt(nom)<errcg){
02562 if (i!=nicg-1){
02563 buff[maxncdofd]=1.0;
02564 for (j=1;j<nproc;j++){
02565 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02566 }
02567 }
02568 break;
02569 }
02570
02571
02572
02573 beta = nom/denom;
02574
02575
02576 for (j=0;j<ndofcp;j++){
02577
02578 d[j]=beta*d[j]-p[j];
02579 }
02580
02581 }
02582
02583 anicg=i; aerrcg=nom;
02584
02585 fprintf (out,"\n\n\n\n kontrola Lagrangeovych multiplikatoru \n");
02586 for (i=0;i<ndofcp;i++){
02587 fprintf (out,"\n lagr. mult %4ld %e",i,w[i]);
02588 }
02589
02590 delete [] d;
02591 delete [] p;
02592 delete [] g;
02593 }
02594 else{
02595
02596 for (i=0;i<nicg;i++){
02597
02598 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02599
02600 if (buff[maxncdofd]>0.5){ break; }
02601
02602 nullv (dd,ndof);
02603
02604 buff_local (dd,buff);
02605
02606
02607 gm->ldl_feti (pp,dd,nrbm,rbmi,zero);
02608
02609
02610 local_buff (pp,buff);
02611
02612 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
02613
02614
02615
02616
02617
02618
02619
02620 switch (fetiprecond){
02621
02622 case nofetiprecond:{
02623
02624
02625
02626
02627
02628
02629
02630
02631
02632
02633
02634
02635
02636
02637
02638
02639
02640 break;
02641 }
02642
02643 case lumped:{
02644 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02645
02646 if (buff[maxncdofd]>0.5){ break; }
02647
02648 nullv (dd,ndof);
02649 buff_local (dd,buff);
02650
02651
02652
02653
02654
02655 nullv (pp,ndof);
02656 lumpedprec (gm,dd,pp,out);
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667 local_buff (pp,buff);
02668 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
02669 break;
02670 }
02671
02672 case dirichlet:{
02673 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02674
02675 if (buff[maxncdofd]>0.5){ break; }
02676
02677 nullv (dd,ndof);
02678 buff_local (dd,buff);
02679
02680
02681
02682
02683
02684 nullv (pp,ndof);
02685 dirichletprec (dd,pp,out);
02686
02687 for (l = 0; l < ndof; l++){
02688 fprintf(out,"pp[%ld] %le\n",l,pp[l]);
02689 }
02690
02691
02692
02693
02694
02695
02696
02697 local_buff (pp,buff);
02698 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
02699 break;
02700 }
02701 }
02702
02703
02704
02705
02706
02707
02708 }
02709
02710 }
02711 MPI_Barrier (MPI_COMM_WORLD);
02712
02713 delete [] dd;
02714 delete [] pp;
02715 delete [] buff;
02716 }
02717
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728
02729
02730
02731
02732
02733
02734 void feti1::lagrmultdispl (gtopology *top,gmatrix *gm,long *domproc,
02735 double *w,double *d,double *f,double *rbm,long *rbmi,
02736 double *h,double *h1,FILE *out)
02737 {
02738 long i,j,k,l,buffsize;
02739 double *g,*pp,*av,*gamma,*buff,*u;
02740 MPI_Status stat;
02741
02742 pp = new double [ndof];
02743 av = new double [maxnrbm];
02744
02745 buffsize=maxncdofd;
02746 buff = new double [buffsize];
02747
02748 if (myrank==0){
02749 g = new double [ndofcp];
02750
02751 for (j=1;j<nproc;j++){
02752 coarse_buff (w,buff,domproc[j]);
02753 MPI_Send (buff,buffsize,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02754 }
02755
02756 coarse_buff (w,buff,domproc[0]);
02757 nullv (pp,ndof);
02758 buff_local (pp,buff);
02759 subv (pp,f,ndof);
02760
02761 gm->ldl_feti (d,pp,nrbm,rbmi,zero);
02762
02763 local_buff (d,buff);
02764
02765 nullv (g,ndofcp);
02766 buff_coarse (g,buff,domproc[0]);
02767
02768 for (j=1;j<nproc;j++){
02769 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02770 buff_coarse (g,buff,domproc[stat.MPI_TAG]);
02771 }
02772 }
02773 else{
02774 MPI_Recv (buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02775 nullv (pp,ndof);
02776 buff_local (pp,buff);
02777
02778 subv (pp,f,ndof);
02779
02780 gm->ldl_feti (d,pp,nrbm,rbmi,zero);
02781
02782 local_buff (d,buff);
02783
02784 MPI_Send (buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
02785 }
02786 MPI_Barrier (MPI_COMM_WORLD);
02787
02788
02789 if (myrank==0){
02790
02791 fprintf (out,"\n");
02792 for (i=0;i<ndofcp;i++){
02793 fprintf (out,"\n g %20.15le",g[i]);
02794 }
02795
02796 gamma = new double [hsize];
02797 u = new double [hsize];
02798
02799 mtxv (h,g,u,ndofcp,hsize);
02800 mxv (h1,u,gamma,hsize,hsize);
02801 cmulv (-1.0,gamma,hsize);
02802
02803 for (j=1;j<nproc;j++){
02804 k=domproc[j];
02805 l=rbmadr[k];
02806 for (i=0;i<nrbmdom[k];i++){
02807 av[i] = gamma[l]; l++;
02808 }
02809
02810 fprintf (out,"\n");
02811 for (i=0;i<nrbmdom[k];i++){
02812 fprintf (out,"\n av %20.15le",av[i]);
02813 }
02814
02815
02816 MPI_Send (av,maxnrbm,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD);
02817 }
02818
02819
02820
02821
02822 k=domproc[0];
02823 l=rbmadr[k];
02824 for (i=0;i<nrbmdom[k];i++){
02825 av[i] = gamma[l]; l++;
02826 }
02827
02828 fprintf (out,"\n");
02829 for (i=0;i<nrbmdom[k];i++){
02830 fprintf (out,"\n av %20.15le",av[i]);
02831 }
02832
02833 delete [] u;
02834 delete [] gamma;
02835 delete [] g;
02836 }
02837 else{
02838 MPI_Recv (av,maxnrbm,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02839 }
02840 MPI_Barrier (MPI_COMM_WORLD);
02841
02842
02843
02844 cmulv (-1.0,d,ndof);
02845
02846 for (i=0;i<ndof;i++){
02847 for (j=0;j<nrbm;j++){
02848 d[i]+=rbm[j*ndof+i]*av[j];
02849 }
02850 }
02851
02852
02853 delete [] buff;
02854 delete [] av; delete [] pp;
02855 }
02856
02857
02858
02859
02860
02861
02862 void feti1::solve_system (gtopology *top,gmatrix *gm,
02863 long *domproc,double *lhs,double *rhs,FILE *out)
02864 {
02865 long i,j;
02866 long *rbmi;
02867 double *rbm,*h,*q,*hh,*ih,*ee,*lm;
02868 time_t t1,t2,t3,t4,t5,t6;
02869
02870
02871
02872 rbm = new double [enrbm*ndof];
02873 memset (rbm,0,enrbm*ndof*sizeof(double));
02874
02875 rbmi = new long [enrbm];
02876 for (i=0;i<enrbm;i++){
02877 rbmi[i]=-1;
02878 }
02879
02880
02881 switch (fetiprecond){
02882 case nofetiprecond:{
02883 fprintf(out,"\n\n\n No preconditioning is used for FETI method\n\n");
02884 break;
02885 }
02886 case lumped:{
02887 fprintf(out,"\n\n\n Lumped preconditioning is used for FETI method\n\n");
02888 subdomain_matrix (gm,domproc,out);
02889 break;
02890 }
02891 case dirichlet:{
02892 fprintf(out,"\n\n\n Dirichlet preconditioning is used for FETI method\n\n");
02893 subdomain_matrix (gm,domproc,out);
02894 break;
02895 }
02896 }
02897
02898
02899 t1 = time (NULL);
02900
02901
02902
02903
02904 gm->kernel (rbm,nrbm,rbmi,enrbm,lithr,3);
02905
02906
02907
02908 t2 = time (NULL);
02909
02910 fprintf (stdout,"\n proc %d nrbm %ld",myrank,nrbm);
02911 fprintf (stdout,"\n proc %d enrbm %ld",myrank,enrbm);
02912 fprintf (stdout,"\n proc %d lithr %e",myrank,lithr);
02913 fprintf (stdout,"\n proc %d ndof %ld",myrank,ndof);
02914
02915 fprintf (out,"\n proc %d nrbm %ld",myrank,nrbm);
02916 fprintf (out,"\n proc %d enrbm %ld",myrank,enrbm);
02917 fprintf (out,"\n proc %d lithr %e",myrank,lithr);
02918 fprintf (out,"\n proc %d ndof %ld",myrank,ndof);
02919
02920
02921
02922 hmatrixsize (rbm,domproc,out);
02923
02924
02925
02926 fprintf (out,"\n\n\n kontrola RBM");
02927 for (i=0;i<ndof;i++){
02928 fprintf (out,"\n");
02929 for (j=0;j<nrbm;j++){
02930 fprintf (out," %f",rbm[j*ndof+i]);
02931 }
02932 }
02933
02934
02935
02936
02937
02938
02939
02940
02941
02942
02943 if (myrank==0){
02944 hsize=rbmadr[nproc];
02945
02946 if (hsize==0){
02947 fprintf (stdout,"\n\n hsize = 0 \n");
02948 fprintf (stderr,"\n\n hsize = 0 \n");
02949 abort ();
02950 }
02951
02952 h = new double [ndofcp*hsize];
02953 q = new double [hsize];
02954 lm = new double [ndofcp];
02955
02956 fprintf (out,"\n hsize %ld",hsize);
02957 fprintf (out,"\n ndofcp %ld",ndofcp);
02958
02959 for (i=0;i<=nproc;i++){
02960 fprintf (out,"\n rbmadr %4ld %4ld",i,rbmadr[i]);
02961 }
02962
02963 }
02964
02965
02966
02967 hmatrix (h,rbm,domproc);
02968
02969
02970 if (myrank==0){
02971 fprintf (out,"\n\n\n kontrola matice H \n");
02972 for (i=0;i<ndofcp;i++){
02973 fprintf (out,"\n %4ld",i);
02974 for (j=0;j<hsize;j++){
02975 fprintf (out," %f",h[i*hsize+j]);
02976 }
02977 }
02978 }
02979
02980
02981 qvector (q,rbm,rhs,domproc);
02982
02983
02984 if (myrank==0){
02985 fprintf (out,"\n\n\n kontrola vektoru q \n");
02986 for (j=0;j<hsize;j++){
02987 fprintf (out,"\n %4ld %f",j,q[j]);
02988 }
02989 }
02990
02991 t3 = time (NULL);
02992
02993
02994 if (myrank==0){
02995 ih = new double [hsize*hsize];
02996 hh = new double [hsize*hsize];
02997 ee = new double [hsize*hsize];
02998 for (i=0;i<hsize;i++){
02999 for (j=0;j<hsize;j++){
03000 ee[i*hsize+j]=0.0;
03001 }
03002 ee[i*hsize+i]=1.0;
03003 }
03004
03005 mtxm (h,h,hh,ndofcp,hsize,hsize);
03006
03007 fprintf (out,"\n\n\n kontrola matice HH \n");
03008 for (i=0;i<hsize;i++){
03009 fprintf (out,"\n %4ld ",i);
03010 for (j=0;j<hsize;j++){
03011 fprintf (out," %lf",hh[i*hsize+j]);
03012 }
03013 }
03014
03015 gemp (hh,ih,ee,hsize,hsize,zero,1);
03016
03017
03018
03019
03020
03021 fprintf (out,"\n\n\n kontrola matice HHinv \n");
03022 for (i=0;i<hsize;i++){
03023 fprintf (out,"\n %4ld ",i);
03024 for (j=0;j<hsize;j++){
03025 fprintf (out," %lf",ih[i*hsize+j]);
03026 }
03027 }
03028
03029 delete [] ee;
03030 }
03031
03032
03033
03034 t4 = time (NULL);
03035
03036
03037 mpcg (top,gm,domproc,lm,rhs,q,h,ih,rbmi,out);
03038
03039 t5 = time (NULL);
03040
03041 lagrmultdispl (top,gm,domproc,lm,lhs,rhs,rbm,rbmi,h,ih,out);
03042
03043 t6 = time (NULL);
03044
03045
03046 fprintf (out,"\n\n evaluation of kernel %ld",t2-t1);
03047 fprintf (out,"\n\n matrix H computation %ld",t3-t2);
03048 fprintf (out,"\n\n decomposition of H %ld",t4-t3);
03049 fprintf (out,"\n\n modified conjug. gradient m. %ld",t5-t4);
03050 fprintf (out,"\n\n displacements computation %ld",t6-t5);
03051
03052 fprintf (stdout,"\n\n evaluation of kernel %ld",t2-t1);
03053 fprintf (stdout,"\n\n matrix H computation %ld",t3-t2);
03054 fprintf (stdout,"\n\n decomposition of H %ld",t4-t3);
03055 fprintf (stdout,"\n\n modified conjug. gradient m. %ld",t5-t4);
03056 fprintf (stdout,"\n\n displacements computation %ld",t6-t5);
03057
03058
03059 MPI_Status stat;
03060 long *buff;
03061 buff = new long[5];
03062
03063 buff[0]=t2-t1;
03064 buff[1]=t3-t2;
03065 buff[2]=t4-t3;
03066 buff[3]=t5-t4;
03067 buff[4]=t6-t5;
03068
03069 if (myrank==0){
03070
03071 fprintf (out,"\n\n\n\n\n");
03072 j=domproc[0];
03073 fprintf (out,"\n\n\n Domain %ld",j);
03074 fprintf (out,"\n evaluation of kernel %ld",buff[0]);
03075 fprintf (out,"\n matrix H computation %ld",buff[1]);
03076 fprintf (out,"\n decomposition of H %ld",buff[2]);
03077 fprintf (out,"\n modified conjug. gradient m. %ld",buff[3]);
03078 fprintf (out,"\n displacements computation %ld",buff[4]);
03079
03080
03081 j=domproc[0];
03082 fprintf (stdout,"\n\n\n Domain %ld",j);
03083 fprintf (stdout,"\n evaluation of kernel %ld",buff[0]);
03084 fprintf (stdout,"\n matrix H computation %ld",buff[1]);
03085 fprintf (stdout,"\n decomposition of H %ld",buff[2]);
03086 fprintf (stdout,"\n modified conjug. gradient m. %ld",buff[3]);
03087 fprintf (stdout,"\n displacements computation %ld",buff[4]);
03088
03089 for (i=1;i<nproc;i++){
03090 MPI_Recv (buff,5,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
03091 j=domproc[stat.MPI_TAG];
03092
03093 fprintf (out,"\n\n\n Domain %ld",j);
03094 fprintf (out,"\n evaluation of kernel %ld",buff[0]);
03095 fprintf (out,"\n matrix H computation %ld",buff[1]);
03096 fprintf (out,"\n decomposition of H %ld",buff[2]);
03097 fprintf (out,"\n modified conjug. gradient m. %ld",buff[3]);
03098 fprintf (out,"\n displacements computation %ld",buff[4]);
03099
03100 fprintf (stdout,"\n\n\n Domain %ld",j);
03101 fprintf (stdout,"\n evaluation of kernel %ld",buff[0]);
03102 fprintf (stdout,"\n matrix H computation %ld",buff[1]);
03103 fprintf (stdout,"\n decomposition of H %ld",buff[2]);
03104 fprintf (stdout,"\n modified conjug. gradient m. %ld",buff[3]);
03105 fprintf (stdout,"\n displacements computation %ld",buff[4]);
03106
03107 }
03108 }
03109 else{
03110 MPI_Send(buff,5,MPI_LONG,0,myrank,MPI_COMM_WORLD);
03111 }
03112 MPI_Barrier (MPI_COMM_WORLD);
03113
03114 delete [] buff;
03115
03116
03117
03118
03119 }