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