00001 #include "dpfeti.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <time.h>
00005 #include <limits.h>
00006
00007 dpfeti::dpfeti(int np,int mr,long nd)
00008 {
00009
00010 nproc=np;
00011
00012 myrank=mr;
00013
00014 ndom=nd;
00015
00016 ssle = new slesolv();
00017 }
00018
00019 dpfeti::~dpfeti()
00020 {
00021 long i;
00022
00023 delete [] cgrhs;
00024 delete [] cglhs;
00025 delete [] displ;
00026 delete [] br;
00027 delete [] bm;
00028
00029 delete [] krc;
00030
00031 for (i=0;i<nproc;i++){
00032 delete [] loccn[i];
00033 delete [] globcn[i];
00034 }
00035 delete [] loccn;
00036 delete [] globcn;
00037
00038 delete [] ncdofdom;
00039 delete [] nlbdofdom;
00040
00041 delete [] localcn;
00042 delete [] globalcn;
00043
00044 }
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060 void dpfeti::cornodedetection (gtopology *top,long *ltg,long *domproc,FILE *out)
00061 {
00062 long i,j,k,nbn,maxnbn,totnbn,totncn,totnrn,stop;
00063 long *buff,**gnn,*cnid,*nbndom;
00064 MPI_Status stat;
00065
00066 nbn=0;
00067 for (i=0;i<top->nn;i++){
00068 if (ltg[i]>-1) nbn++;
00069 }
00070
00071 if (myrank==0){
00072 nbndom = new long [nproc];
00073 j=domproc[0];
00074 nbndom[j]=nbn;
00075
00076
00077 maxnbn=nbn;
00078 for (i=1;i<nproc;i++){
00079 MPI_Recv (&nbn,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00080 j=domproc[stat.MPI_TAG];
00081 nbndom[j]=nbn;
00082
00083 if (maxnbn<nbn) maxnbn=nbn;
00084 }
00085 for (i=1;i<nproc;i++){
00086 MPI_Send (&maxnbn,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00087 }
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 }
00098 else{
00099 MPI_Send (&nbn,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00100 MPI_Recv (&maxnbn,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00101 }
00102 MPI_Barrier (MPI_COMM_WORLD);
00103
00104
00105 buff = new long [maxnbn];
00106 j=0;
00107 for (i=0;i<top->nn;i++){
00108 if (ltg[i]>-1){
00109 buff[j]=ltg[i];
00110 j++;
00111 }
00112 }
00113
00114
00115
00116 if (myrank==0){
00117
00118 gnn = new long* [nproc];
00119 for (i=0;i<nproc;i++){
00120 gnn[i] = new long [nbndom[i]];
00121 }
00122
00123
00124
00125 j=domproc[0];
00126 for (k=0;k<nbndom[j];k++){
00127 gnn[j][k]=buff[k];
00128 }
00129
00130 for (i=1;i<nproc;i++){
00131 MPI_Recv (buff,maxnbn,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00132 j=domproc[stat.MPI_TAG];
00133
00134 for (k=0;k<nbndom[j];k++){
00135 gnn[j][k]=buff[k];
00136 }
00137
00138 }
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 totnbn=0;
00153 for (i=0;i<nproc;i++){
00154 for (j=0;j<nbndom[i];j++){
00155 if (totnbn<gnn[i][j]) totnbn=gnn[i][j];
00156 }
00157 }
00158
00159 totnbn++;
00160
00161
00162
00163
00164
00165 cnid = new long [totnbn];
00166 for (i=0;i<totnbn;i++){
00167 cnid[i]=0;
00168 }
00169
00170
00171 for (i=0;i<nproc;i++){
00172 for (j=0;j<nbndom[i];j++){
00173 cnid[gnn[i][j]]++;
00174 }
00175 }
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 totncn=1; totnrn=1;
00188 for (i=0;i<totnbn;i++){
00189 if (cnid[i]>2){
00190 cnid[i]=0-totncn;
00191 totncn++;
00192 }
00193 else{
00194 cnid[i]=totnrn;
00195 totnrn++;
00196 }
00197 }
00198 totncn--; totnrn--;
00199
00200
00201 fprintf (out,"\n\n pocet corner points %ld",totncn);
00202 fprintf (out,"\n\n pocet remaining points %ld",totnrn);
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 for (i=1;i<nproc;i++){
00213 j=domproc[i];
00214 for (k=0;k<nbndom[j];k++){
00215 buff[k]=cnid[gnn[j][k]];
00216 }
00217 MPI_Send (buff,maxnbn,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00218 }
00219
00220 j=domproc[0];
00221 for (k=0;k<nbndom[j];k++){
00222 buff[k]=cnid[gnn[j][k]];
00223 }
00224
00225 }
00226 else{
00227
00228 MPI_Send (buff,maxnbn,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00229 MPI_Recv (buff,maxnbn,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00230
00231 }
00232 MPI_Barrier (MPI_COMM_WORLD);
00233
00234
00235
00236
00237 j=0;
00238 for(i=0;i<top->nn;i++){
00239 if (ltg[i]>-1){
00240 ltg[i]=buff[j]-1;
00241 j++;
00242 }
00243
00244 if (ltg[i]>=-1) top->gnodes[i].ai = 0;
00245 if (ltg[i]<-1) top->gnodes[i].ai = 0-ltg[i]-1;
00246 }
00247
00248 j=0;
00249 for(i=0;i<top->nn;i++){
00250 if (ltg[i]<-1){
00251 j++;
00252 }
00253 }
00254
00255 fprintf (stdout,"\n pocet corner points %ld",j);
00256
00257
00258 stop=0;
00259 if (myrank==0){
00260 if (j<3) stop=1;
00261 for (i=1;i<nproc;i++){
00262 MPI_Recv (&j,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00263 if (j<3) stop=1;
00264 }
00265 for (i=1;i<nproc;i++){
00266 MPI_Send (&stop,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00267 }
00268 }
00269 else{
00270 MPI_Send (&j,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00271 MPI_Recv (&stop,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00272 }
00273 MPI_Barrier (MPI_COMM_WORLD);
00274
00275
00276 if (stop==1){
00277 fprintf (stderr,"\n\n insufficient number of corner points (domain %ld, number of corner points %ld).\n",ndom,j);
00278 fprintf (stderr,"\n\n insufficient number of corner points (domain %ld, number of corner points %ld).\n",ndom,j);
00279
00280 abort ();
00281 }
00282
00283
00284 if (myrank==0){
00285 for (i=0;i<nproc;i++){
00286 delete [] gnn[i];
00287 }
00288 delete [] gnn;
00289
00290 delete [] cnid;
00291 delete [] nbndom;
00292 }
00293 delete [] buff;
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 void dpfeti::dpfetiordering (gtopology *top,long *ltg)
00323 {
00324 long i,j,k,ndofn;
00325
00326
00327 ndof=1;
00328 for (i=0;i<top->nn;i++){
00329 if (ltg[i]>=-1){
00330
00331 ndofn=top->give_ndofn (i);
00332 for (j=0;j<ndofn;j++){
00333 k=top->give_dof(i,j);
00334 if (k<0) continue;
00335 if (k==0) continue;
00336 if (k>0){
00337 top->save_dof (i,j,ndof); ndof++;
00338 }
00339 }
00340 }
00341 }
00342 nidof=ndof-1;
00343 top->nidof=nidof;
00344
00345 for (i=0;i<top->nn;i++){
00346 if (ltg[i]<-1){
00347
00348 ndofn=top->give_ndofn (i);
00349 for (j=0;j<ndofn;j++){
00350 k=top->give_dof (i,j);
00351 if (k<0) continue;
00352 if (k==0) continue;
00353 if (k>0){
00354 top->save_dof (i,j,ndof); ndof++;
00355 }
00356 }
00357 }
00358 }
00359 ndof--;
00360
00361
00362 for (i=0;i<top->ne;i++){
00363 if (top->give_cne (i)==1){
00364 for (j=0;j<top->give_ndofe(i);j++){
00365 if (top->give_dof (i,j)>ndof) ndof=top->give_dof (i,j);
00366 }
00367 }
00368 }
00369
00370 ncdof=ndof-nidof;
00371 top->nbdof=ncdof;
00372 }
00373
00374
00375
00376
00377 void dpfeti::initiate (selnodes *selnodschur,selnodes *selnodfeti)
00378 {
00379
00380 nlbdof = selnodfeti->ndof;
00381
00382 maxnbdof = selnodfeti->maxndof;
00383
00384 localcn = selnodfeti->ldof;
00385
00386
00387 if (myrank==0){
00388
00389 tnmdof = selnodfeti->tndof;
00390
00391 nlbdofdom = selnodfeti->ndofdom;
00392
00393 globcn = selnodfeti->cndom;
00394 }
00395
00396
00397
00398
00399 maxncdof = selnodschur->maxndof;
00400
00401 ncdof = selnodschur->ndof;
00402
00403 if (myrank==0){
00404
00405 tncdof = selnodschur->tndof;
00406
00407 ncdofdom=selnodschur->ndofdom;
00408
00409 ptop = new gtopology;
00410 ptop->initiate (selnodschur->cndom,selnodschur->ndofdom,nproc);
00411 ptop->cnstate=1;
00412 }
00413 }
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423 void dpfeti::globcnnum_dpfeti (gtopology *top,long *ltg,long *domproc,FILE *out)
00424 {
00425 long i,j,k,l,m,n,ndofn,stop,nbn,ncn,maxnbn,maxncn,tnbn,tncn;
00426 long *nibn,*lgnbn,*llnbn,*lgncn,*llncn,***lcnbn,***lcncn,**gnbndom,**gncndom;
00427 long **cncn,***cnbn,*nbndom,*ncndom,**dominc,**supelemcn;
00428 long buffsize,*buff;
00429 MPI_Status stat;
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442 ndofn = top->give_ndofn (0);
00443
00444
00445
00446 nbn=0; ncn=0;
00447 for (i=0;i<top->nn;i++){
00448 if (ltg[i]>-1) nbn++;
00449 if (ltg[i]<-1) ncn++;
00450 }
00451
00452
00453
00454 lgnbn = new long [nbn];
00455 lgncn = new long [ncn];
00456
00457
00458 llnbn = new long [nbn];
00459 llncn = new long [ncn];
00460
00461 j=0; k=0;
00462 for (i=0;i<top->nn;i++){
00463 if (ltg[i]>-1){
00464 lgnbn[j]=ltg[i];
00465 llnbn[j]=i;
00466 j++;
00467 }
00468 if (ltg[i]<-1){
00469 lgncn[k]=0-ltg[i]-2;
00470 llncn[k]=i;
00471 k++;
00472 }
00473 }
00474
00475 buff = new long [2];
00476 buff[0]=nbn;
00477 buff[1]=ncn;
00478
00479 if (myrank==0){
00480
00481 nbndom = new long [nproc];
00482
00483 ncndom = new long [nproc];
00484
00485 nbndom[domproc[0]]=nbn;
00486 ncndom[domproc[0]]=ncn;
00487
00488
00489
00490 maxnbn=nbn; maxncn=ncn;
00491 for (i=1;i<nproc;i++){
00492 MPI_Recv (buff,2,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00493 j=domproc[stat.MPI_TAG];
00494
00495 nbndom[j]=buff[0];
00496 ncndom[j]=buff[1];
00497
00498 if (maxnbn<buff[0]) maxnbn=buff[0];
00499 if (maxncn<buff[1]) maxncn=buff[1];
00500 }
00501
00502 buff[0]=maxnbn;
00503 buff[1]=maxncn;
00504 for (i=1;i<nproc;i++){
00505 MPI_Send (buff,2,MPI_LONG,i,myrank,MPI_COMM_WORLD);
00506 }
00507 }
00508 else{
00509 MPI_Send (buff,2,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00510 MPI_Recv (buff,2,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00511 maxnbn=buff[0];
00512 maxncn=buff[1];
00513 }
00514 MPI_Barrier (MPI_COMM_WORLD);
00515 delete [] buff;
00516
00517
00518 buffsize = (maxnbn+maxncn)*(ndofn+1);
00519 buff = new long [buffsize];
00520 k=0;
00521 for (i=0;i<nbn;i++){
00522 buff[k]=lgnbn[i]; k++;
00523 for (j=0;j<ndofn;j++){
00524 buff[k]=top->give_dof (llnbn[i],j);
00525 k++;
00526 }
00527 }
00528 for (i=0;i<ncn;i++){
00529 buff[k]=lgncn[i]; k++;
00530 for (j=0;j<ndofn;j++){
00531 buff[k]=top->give_dof (llncn[i],j);
00532 k++;
00533 }
00534 }
00535
00536 delete [] llnbn;
00537 delete [] llncn;
00538 delete [] lgnbn;
00539 delete [] lgncn;
00540
00541 if (myrank==0){
00542
00543
00544 lcnbn = new long** [nproc];
00545
00546 lcncn = new long** [nproc];
00547
00548 gnbndom = new long* [nproc];
00549
00550 gncndom = new long* [nproc];
00551 for (i=0;i<nproc;i++){
00552 gnbndom[i] = new long [nbndom[i]];
00553 gncndom[i] = new long [ncndom[i]];
00554
00555 lcnbn[i] = new long* [nbndom[i]];
00556 lcncn[i] = new long* [ncndom[i]];
00557 for (j=0;j<nbndom[i];j++){
00558 lcnbn[i][j] = new long [ndofn];
00559 }
00560 for (j=0;j<ncndom[i];j++){
00561 lcncn[i][j] = new long [ndofn];
00562 }
00563 }
00564
00565
00566 l=0;
00567 for (j=0;j<nbndom[domproc[0]];j++){
00568 gnbndom[domproc[0]][j]=buff[l]; l++;
00569 for (k=0;k<ndofn;k++){
00570 lcnbn[domproc[0]][j][k]=buff[l];
00571 l++;
00572 }
00573 }
00574 for (j=0;j<ncndom[domproc[0]];j++){
00575 gncndom[domproc[0]][j]=buff[l]; l++;
00576 for (k=0;k<ndofn;k++){
00577 lcncn[domproc[0]][j][k]=buff[l];
00578 l++;
00579 }
00580 }
00581
00582
00583 for (i=1;i<nproc;i++){
00584 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00585 m=domproc[stat.MPI_TAG];
00586 l=0;
00587 for (j=0;j<nbndom[m];j++){
00588 gnbndom[m][j]=buff[l]; l++;
00589 for (k=0;k<ndofn;k++){
00590 lcnbn[m][j][k]=buff[l];
00591 l++;
00592 }
00593 }
00594 for (j=0;j<ncndom[m];j++){
00595 gncndom[m][j]=buff[l]; l++;
00596 for (k=0;k<ndofn;k++){
00597 lcncn[m][j][k]=buff[l];
00598 l++;
00599 }
00600 }
00601
00602 }
00603 }
00604 else{
00605 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00606 }
00607 MPI_Barrier (MPI_COMM_WORLD);
00608 delete [] buff;
00609
00610
00611
00612 if (myrank==0){
00613
00614
00615
00616 tnbn = 0;
00617 for (i=0;i<nproc;i++){
00618 for (j=0;j<nbndom[i];j++){
00619 if (gnbndom[i][j]>tnbn) tnbn=gnbndom[i][j];
00620 }
00621 }
00622 tnbn++;
00623
00624 tncn = 0;
00625 for (i=0;i<nproc;i++){
00626 for (j=0;j<ncndom[i];j++){
00627 if (gncndom[i][j]>tncn) tncn=gncndom[i][j];
00628 }
00629 }
00630 tncn++;
00631
00632
00633 cncn = new long* [tncn];
00634 for (i=0;i<tncn;i++){
00635 cncn[i] = new long [ndofn];
00636 for (j=0;j<ndofn;j++){
00637 cncn[i][j]=0;
00638 }
00639 }
00640
00641
00642 for (i=0;i<nproc;i++){
00643 for (j=0;j<ncndom[i];j++){
00644 for (k=0;k<ndofn;k++){
00645 cncn[gncndom[i][j]][k]=lcncn[i][j][k];
00646 }
00647 }
00648 }
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661 tncdof=1;
00662 for (i=0;i<tncn;i++){
00663 for (j=0;j<ndofn;j++){
00664 if (cncn[i][j]>0){
00665 cncn[i][j]=tncdof;
00666 tncdof++;
00667 }
00668 }
00669 }
00670 tncdof--;
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683 ncdofdom = new long [nproc];
00684 for (i=0;i<nproc;i++){
00685 ncdofdom[i]=0;
00686 }
00687
00688
00689 for (i=0;i<nproc;i++){
00690 for (j=0;j<ncndom[i];j++){
00691 for (k=0;k<ndofn;k++){
00692 if (cncn[gncndom[i][j]][k]>0) ncdofdom[i]++;
00693 }
00694 }
00695 }
00696
00697
00698
00699
00700
00701
00702
00703
00704 maxncdof=0;
00705 for (i=0;i<nproc;i++){
00706 if (maxncdof<ncdofdom[i]) maxncdof=ncdofdom[i];
00707 }
00708
00709
00710
00711 supelemcn = new long* [nproc];
00712 for (i=0;i<nproc;i++){
00713 supelemcn[i] = new long [ncdofdom[i]];
00714 for (j=0;j<ncdofdom[i];j++){
00715 supelemcn[i][j]=0;
00716 }
00717 }
00718
00719
00720 for (i=0;i<nproc;i++){
00721 l=0;
00722 for (j=0;j<ncndom[i];j++){
00723 for (k=0;k<ndofn;k++){
00724 m=cncn[gncndom[i][j]][k];
00725 if (m>0){
00726 supelemcn[i][l]=m;
00727 l++;
00728 }
00729 }
00730 }
00731 }
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743 ptop = new gtopology;
00744 ptop->initiate (supelemcn,ncdofdom,nproc);
00745
00746 for (i=0;i<nproc;i++){
00747 delete [] gncndom[i];
00748 }
00749 delete [] gncndom;
00750
00751 for (i=0;i<nproc;i++){
00752 delete [] supelemcn[i];
00753 }
00754 delete [] supelemcn;
00755
00756 for (i=0;i<tncn;i++){
00757 delete [] cncn[i];
00758 }
00759 delete [] cncn;
00760
00761
00762
00763
00764
00765
00766 nibn = new long [tnbn];
00767 for (i=0;i<tnbn;i++){
00768 nibn[i]=0;
00769 }
00770
00771
00772 for (i=0;i<nproc;i++){
00773 for (j=0;j<nbndom[i];j++){
00774 nibn[gnbndom[i][j]]++;
00775 }
00776 }
00777
00778
00779 dominc = new long* [tnbn];
00780 for (i=0;i<tnbn;i++){
00781 dominc[i] = new long [nibn[i]];
00782 for (j=0;j<nibn[i];j++){
00783 dominc[i][j]=-1;
00784 }
00785 }
00786
00787 for (i=0;i<tnbn;i++){
00788 nibn[i]=0;
00789 }
00790
00791
00792 for (i=0;i<nproc;i++){
00793 for (j=0;j<nbndom[i];j++){
00794 dominc[gnbndom[i][j]][nibn[gnbndom[i][j]]]=i;
00795 nibn[gnbndom[i][j]]++;
00796 }
00797 }
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811 cnbn = new long** [tnbn];
00812 for (i=0;i<tnbn;i++){
00813 cnbn[i] = new long* [nibn[i]-1];
00814 for (j=0;j<nibn[i]-1;j++){
00815 cnbn[i][j] = new long [ndofn];
00816 for (k=0;k<ndofn;k++){
00817 cnbn[i][j][k]=0;
00818 }
00819 }
00820 }
00821
00822
00823 for (i=0;i<nproc;i++){
00824 for (j=0;j<nbndom[i];j++){
00825 m=gnbndom[i][j];
00826 for (k=0;k<nibn[m]-1;k++){
00827 for (l=0;l<ndofn;l++){
00828 cnbn[m][k][l]=lcnbn[i][j][l];
00829 }
00830 }
00831 }
00832 }
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848 tnmdof=1;
00849 for (i=0;i<tnbn;i++){
00850 for (j=0;j<nibn[i]-1;j++){
00851 for (k=0;k<ndofn;k++){
00852 if (cnbn[i][j][k]>0){
00853 cnbn[i][j][k]=tnmdof;
00854 tnmdof++;
00855 }
00856 }
00857 }
00858 }
00859 tnmdof--;
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874 nlbdofdom = new long [nproc];
00875 for (i=0;i<nproc;i++){
00876 nlbdofdom[i]=0;
00877 }
00878 for (i=0;i<nproc;i++){
00879 for (j=0;j<nbndom[i];j++){
00880 m=gnbndom[i][j];
00881 if (dominc[m][0]==i){
00882 for (k=0;k<nibn[m]-1;k++){
00883 for (l=0;l<ndofn;l++){
00884 if (cnbn[m][k][l]>0) nlbdofdom[i]++;
00885 }
00886 }
00887 }
00888 else{
00889 for (l=0;l<ndofn;l++){
00890 if (cnbn[m][0][l]>0) nlbdofdom[i]++;
00891 }
00892 }
00893 }
00894 }
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905 loccn = new long* [nproc];
00906 globcn = new long* [nproc];
00907 maxnbdof=0;
00908 for (i=0;i<nproc;i++){
00909 loccn[i] = new long [nlbdofdom[i]];
00910 globcn[i] = new long [nlbdofdom[i]];
00911 if (maxnbdof<nlbdofdom[i]) maxnbdof=nlbdofdom[i];
00912 }
00913
00914 for (i=0;i<nproc;i++){
00915 n=0;
00916 for (j=0;j<nbndom[i];j++){
00917 m=gnbndom[i][j];
00918 if (dominc[m][0]==i){
00919 for (k=0;k<nibn[m]-1;k++){
00920 for (l=0;l<ndofn;l++){
00921 if (cnbn[m][k][l]>0){
00922 globcn[i][n]=0-cnbn[m][k][l];
00923 loccn[i][n]=lcnbn[i][j][l];
00924 n++;
00925 }
00926 }
00927 }
00928 }
00929 else{
00930 stop=0;
00931 for (k=0;k<nibn[m]-1;k++){
00932 for (l=0;l<ndofn;l++){
00933 if (cnbn[m][k][l]==0){
00934 cnbn[m][k][k]=-3;
00935 stop=1;
00936 }
00937 if (cnbn[m][k][l]>0){
00938 globcn[i][n]=cnbn[m][k][l];
00939 loccn[i][n]=lcnbn[i][j][l];
00940 n++;
00941 cnbn[m][k][k]=-3;
00942 stop=1;
00943 }
00944 }
00945 if (stop==1) break;
00946 }
00947 }
00948 }
00949 }
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966 for (i=0;i<nproc;i++){
00967 for (j=0;j<nbndom[i];j++){
00968 delete [] lcnbn[i][j];
00969 }
00970 delete [] lcnbn[i];
00971 }
00972 delete [] lcnbn;
00973
00974 for (i=0;i<nproc;i++){
00975 for (j=0;j<ncndom[i];j++){
00976 delete [] lcncn[i][j];
00977 }
00978 delete [] lcncn[i];
00979 }
00980 delete [] lcncn;
00981
00982 delete [] nbndom;
00983
00984 delete [] ncndom;
00985
00986 for (i=0;i<nproc;i++){
00987 delete [] gnbndom[i];
00988 }
00989 delete [] gnbndom;
00990
00991 for (i=0;i<tnbn;i++){
00992 delete [] dominc[i];
00993 }
00994 delete [] dominc;
00995
00996 for (i=0;i<tnbn;i++){
00997 for (j=0;j<nibn[i]-1;j++){
00998 delete [] cnbn[i][j];
00999 }
01000 delete [] cnbn[i];
01001 }
01002 delete [] cnbn;
01003
01004 delete [] nibn;
01005
01006 }
01007
01008 buff = new long [3];
01009 if (myrank==0){
01010 buff[0]=maxncdof;
01011 buff[1]=maxnbdof;
01012
01013 for (i=1;i<nproc;i++){
01014 buff[2]=nlbdofdom[domproc[i]];
01015 MPI_Send (buff,3,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01016 }
01017 nlbdof=nlbdofdom[domproc[0]];
01018 }
01019 else{
01020 MPI_Recv (buff,3,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01021 maxncdof=buff[0];
01022 maxnbdof=buff[1];
01023 nlbdof=buff[2];
01024 }
01025 MPI_Barrier (MPI_COMM_WORLD);
01026 delete [] buff;
01027
01028 localcn = new long [nlbdof];
01029 globalcn = new long [nlbdof];
01030
01031 buffsize=maxnbdof;
01032 buff = new long [buffsize];
01033 if (myrank==0){
01034 for (i=1;i<nproc;i++){
01035 k=domproc[i];
01036 for (j=0;j<nlbdofdom[k];j++){
01037 buff[j]=loccn[k][j];
01038 }
01039 MPI_Send (buff,buffsize,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01040 }
01041 k=domproc[0];
01042 for (j=0;j<nlbdofdom[k];j++){
01043 buff[j]=loccn[k][j];
01044 }
01045 }
01046 else{
01047 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01048 }
01049 MPI_Barrier (MPI_COMM_WORLD);
01050
01051 for (j=0;j<nlbdof;j++){
01052 localcn[j]=buff[j];
01053 }
01054
01055 delete [] buff;
01056
01057
01058
01059
01060
01061
01062
01063 }
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075 void dpfeti::extract_from_local_vector (double *ev,double *lv)
01076 {
01077 long i,lcn;
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087 for (i=0;i<nlbdof;i++){
01088 lcn=localcn[i];
01089 if (lcn<0){
01090 lcn=0-lcn-1;
01091 ev[i]=0.0-lv[lcn];
01092 }
01093 else{
01094 lcn=lcn-1;
01095 ev[i]=lv[lcn];
01096 }
01097 }
01098
01099 }
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110 void dpfeti::put_into_local_vector (double *ev,double *lv)
01111 {
01112 long i,lcn;
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122 for (i=0;i<nlbdof;i++){
01123 lcn=localcn[i];
01124 if (lcn<0){
01125 lcn=0-lcn-1;
01126 lv[lcn]-=ev[i];
01127 }
01128 else{
01129 lcn=lcn-1;
01130 lv[lcn]+=ev[i];
01131 }
01132 }
01133
01134 }
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147 void dpfeti::extract_from_global_vector (double *ev,double *gv,long nsub)
01148 {
01149 long i,gcn;
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161 for (i=0;i<nlbdofdom[nsub];i++){
01162 gcn=globcn[nsub][i]-1;
01163 ev[i]=gv[gcn];
01164 }
01165
01166 }
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178 void dpfeti::put_into_global_vector (double *ev,double *gv,long nsub)
01179 {
01180 long i,gcn;
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192 for (i=0;i<nlbdofdom[nsub];i++){
01193 gcn=globcn[nsub][i]-1;
01194 gv[gcn]+=ev[i];
01195 }
01196
01197 }
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213 void dpfeti::arrmatrix (double *condmat,long *domproc,FILE *out)
01214 {
01215 long i,j,k,l,ii,gndofe,*cn;
01216 double *buff;
01217 MPI_Status stat;
01218
01219 buff = new double [maxncdof*maxncdof];
01220 ii=0;
01221 for (i=0;i<ncdof;i++){
01222 for (j=0;j<ncdof;j++){
01223 buff[ii]=condmat[i*ncdof+j]; ii++;
01224 }
01225 }
01226
01227 if (myrank==0){
01228 arr = new gmatrix;
01229
01230 arr->setval (ssle);
01231 arr->initiate (ptop,tncdof,rsmstor,1,out);
01232 arr->prepmat (0.0,1);
01233
01234 gndofe = ptop->give_ndofe (domproc[0]);
01235 cn = new long [gndofe];
01236 ptop->give_code_numbers (domproc[0],cn);
01237 arr->localized (condmat,cn,gndofe,domproc[0]);
01238 delete [] cn;
01239
01240 for (i=1;i<nproc;i++){
01241 MPI_Recv(buff,maxncdof*maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01242 l=domproc[stat.MPI_TAG];
01243
01244 ii=0;
01245 for (j=0;j<ncdofdom[l];j++){
01246 for (k=0;k<ncdofdom[l];k++){
01247 condmat[j*ncdofdom[l]+k]=buff[ii]; ii++;
01248 }
01249 }
01250
01251
01252
01253
01254
01255
01256
01257
01258 gndofe = ptop->give_ndofe (l);
01259 cn = new long [gndofe];
01260 ptop->give_code_numbers (l,cn);
01261 arr->localized (condmat,cn,gndofe,l);
01262 delete [] cn;
01263
01264 }
01265
01266 arr->decompose_matrix ();
01267
01268
01269 }
01270 else{
01271 MPI_Send(buff,maxncdof*maxncdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01272 }
01273 MPI_Barrier (MPI_COMM_WORLD);
01274 delete [] buff;
01275 }
01276
01277
01278 void dpfeti::vectors_br_bm (gtopology *top,gmatrix *gm,
01279 double *condmat,double *condvect,double *rhs,long *domproc,FILE *out)
01280 {
01281 long i,j,l,gndofe,*cn;
01282 double *lhs,*buff,*zaloha;
01283 MPI_Status stat;
01284
01285
01286
01287
01288 buff = new double [maxncdof];
01289 for (i=0;i<ncdof;i++){
01290 buff[i]=condvect[i];
01291 }
01292
01293 if (myrank==0){
01294 br = new double [tncdof];
01295 nullv (br,tncdof);
01296
01297 gndofe = ptop->give_ndofe (domproc[0]);
01298 cn = new long [gndofe];
01299 ptop->give_code_numbers (domproc[0],cn);
01300 locglob (br,condvect,cn,gndofe);
01301 delete [] cn;
01302
01303 for (i=1;i<nproc;i++){
01304 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01305 l=domproc[stat.MPI_TAG];
01306
01307 for (j=0;j<ncdofdom[l];j++){
01308 condvect[j]=buff[j];
01309 }
01310
01311 gndofe = ptop->give_ndofe (l);
01312 cn = new long [gndofe];
01313 ptop->give_code_numbers (l,cn);
01314 locglob (br,condvect,cn,gndofe);
01315 delete [] cn;
01316
01317 }
01318
01319
01320
01321
01322
01323
01324
01325
01326 }
01327 else{
01328 MPI_Send(buff,maxncdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01329 }
01330 MPI_Barrier (MPI_COMM_WORLD);
01331 delete [] buff;
01332
01333
01334
01335
01336
01337
01338
01339 buff = new double [maxnbdof];
01340
01341
01342
01343
01344
01345
01346
01347 lhs = new double [ndof];
01348 zaloha = new double [ndof];
01349 for (i=0;i<nidof;i++){
01350 zaloha[i]=rhs[i];
01351 }
01352
01353
01354
01355
01356
01357
01358
01359 gm->condense (top,condmat,buff,lhs,zaloha,ncdof,3,out);
01360 delete [] zaloha;
01361
01362
01363
01364 nullv (buff,maxnbdof);
01365 extract_from_local_vector (buff,lhs);
01366
01367 if (myrank==0){
01368 bm = new double [tnmdof];
01369 nullv (bm,tnmdof);
01370 put_into_global_vector (buff,bm,domproc[0]);
01371 for (j=1;j<nproc;j++){
01372 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01373 put_into_global_vector (buff,bm,domproc[stat.MPI_TAG]);
01374 }
01375
01376
01377
01378
01379
01380
01381
01382
01383 }
01384 else{
01385 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01386 }
01387 MPI_Barrier (MPI_COMM_WORLD);
01388
01389
01390 delete [] buff;
01391 delete [] lhs;
01392
01393
01394 }
01395
01396
01397 void dpfeti::rhs_dpfeti (gtopology *top,long *domproc,gmatrix *gm,FILE *out)
01398 {
01399 long i,j,gndofe,*cn;
01400 double *ur,*wr,*utc,*vtc,*buff;
01401 MPI_Status stat;
01402
01403
01404
01405 ur = new double [ndof];
01406 wr = new double [ndof];
01407 buff = new double [maxncdof];
01408
01409 if (myrank==0){
01410 utc = new double [tncdof];
01411 vtc = new double [tncdof];
01412
01413 for (i=0;i<tncdof;i++){
01414 utc[i]=br[i];
01415 }
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425 arr->back_substitution (vtc,utc);
01426
01427
01428
01429
01430
01431
01432
01433
01434 for (i=1;i<nproc;i++){
01435
01436 nullv (buff,maxncdof);
01437 gndofe = ptop->give_ndofe (domproc[i]);
01438 cn = new long [gndofe];
01439 ptop->give_code_numbers (domproc[i],cn);
01440 globloc (vtc,buff,cn,gndofe);
01441 delete [] cn;
01442
01443 MPI_Send(buff,maxncdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01444 }
01445
01446 nullv (buff,maxncdof);
01447 gndofe = ptop->give_ndofe (domproc[0]);
01448 cn = new long [gndofe];
01449 ptop->give_code_numbers (domproc[0],cn);
01450 globloc (vtc,buff,cn,gndofe);
01451 delete [] cn;
01452
01453 delete [] utc;
01454 delete [] vtc;
01455 }
01456 else{
01457 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01458 }
01459 MPI_Barrier (MPI_COMM_WORLD);
01460
01461
01462
01463 mxv (krc,buff,ur,nidof,ncdof);
01464 delete [] buff;
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474 gm->condense (top,buff,buff,wr,ur,ncdof,3,out);
01475
01476
01477
01478
01479
01480
01481
01482
01483 buff = new double [maxnbdof];
01484
01485
01486 nullv (buff,maxnbdof);
01487 extract_from_local_vector (buff,wr);
01488
01489 if (myrank==0){
01490 nullv (cgrhs,tnmdof);
01491 put_into_global_vector (buff,cgrhs,domproc[0]);
01492 for (j=1;j<nproc;j++){
01493 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01494 put_into_global_vector (buff,cgrhs,domproc[stat.MPI_TAG]);
01495 }
01496 }
01497 else{
01498 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01499 }
01500 MPI_Barrier (MPI_COMM_WORLD);
01501
01502 delete [] buff;
01503 delete [] ur;
01504 delete [] wr;
01505
01506 if (myrank==0){
01507 for (i=0;i<tnmdof;i++){
01508 cgrhs[i]=bm[i]-cgrhs[i];
01509 }
01510
01511
01512
01513
01514
01515
01516
01517
01518 }
01519
01520 }
01521
01522 void dpfeti::matxvect (gtopology *top,long *domproc,
01523 gmatrix *gm,double *input,double *output,FILE *out)
01524 {
01525 long i,j,gndofe,*cn;
01526 double *utc,*vtc,*ur,*vr,*wr,*buff;
01527 MPI_Status stat;
01528
01529 ur = new double [ndof];
01530 vr = new double [ndof];
01531 wr = new double [ndof];
01532
01533
01534
01535
01536
01537 buff = new double [maxnbdof];
01538
01539 if (myrank==0){
01540 for (i=1;i<nproc;i++){
01541 extract_from_global_vector (buff,input,domproc[i]);
01542 MPI_Send(buff,maxnbdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01543 }
01544 extract_from_global_vector (buff,input,domproc[0]);
01545 }
01546 else{
01547 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01548 }
01549 MPI_Barrier (MPI_COMM_WORLD);
01550
01551
01552 nullv (ur,nidof);
01553 put_into_local_vector (buff,ur);
01554 delete [] buff;
01555
01556
01557 gm->condense (top,buff,buff,vr,ur,ncdof,3,out);
01558
01559 buff = new double [maxncdof];
01560
01561
01562 mtxv (krc,vr,buff,nidof,ncdof);
01563
01564 if (myrank==0){
01565 utc = new double [tncdof];
01566 vtc = new double [tncdof];
01567 nullv (utc,tncdof);
01568
01569
01570 gndofe = ptop->give_ndofe (domproc[0]);
01571 cn = new long [gndofe];
01572 ptop->give_code_numbers (domproc[0],cn);
01573 locglob (utc,buff,cn,gndofe);
01574 delete [] cn;
01575
01576 for (i=1;i<nproc;i++){
01577 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01578
01579 j=domproc[stat.MPI_TAG];
01580 gndofe = ptop->give_ndofe (j);
01581 cn = new long [gndofe];
01582 ptop->give_code_numbers (j,cn);
01583 locglob (utc,buff,cn,gndofe);
01584 delete [] cn;
01585 }
01586
01587
01588 arr->back_substitution (vtc,utc);
01589
01590 for (i=1;i<nproc;i++){
01591
01592 nullv (buff,maxncdof);
01593 gndofe = ptop->give_ndofe (domproc[i]);
01594 cn = new long [gndofe];
01595 ptop->give_code_numbers (domproc[i],cn);
01596 globloc (vtc,buff,cn,gndofe);
01597 delete [] cn;
01598
01599 MPI_Send(buff,maxncdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01600 }
01601
01602 nullv (buff,maxncdof);
01603 gndofe = ptop->give_ndofe (domproc[0]);
01604 cn = new long [gndofe];
01605 ptop->give_code_numbers (domproc[0],cn);
01606 globloc (vtc,buff,cn,gndofe);
01607 delete [] cn;
01608
01609 delete [] utc;
01610 delete [] vtc;
01611 }
01612 else{
01613 MPI_Send(buff,maxncdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01614 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01615 }
01616 MPI_Barrier (MPI_COMM_WORLD);
01617
01618
01619 mxv (krc,buff,ur,nidof,ncdof);
01620 delete [] buff;
01621
01622
01623 gm->condense (top,buff,buff,wr,ur,ncdof,3,out);
01624
01625
01626
01627
01628
01629 for (i=0;i<nidof;i++){
01630 wr[i]=vr[i]+wr[i];
01631 }
01632
01633 buff = new double [maxnbdof];
01634
01635
01636 extract_from_local_vector (buff,wr);
01637
01638 if (myrank==0){
01639 nullv (output,tnmdof);
01640 put_into_global_vector (buff,output,domproc[0]);
01641 for (j=1;j<nproc;j++){
01642 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01643 put_into_global_vector (buff,output,domproc[stat.MPI_TAG]);
01644 }
01645 }
01646 else{
01647 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01648 }
01649 delete [] buff;
01650 delete [] ur;
01651 delete [] vr;
01652 delete [] wr;
01653 }
01654
01655
01656
01657
01658
01659
01660 void dpfeti::cg (gtopology *top,long *domproc,gmatrix *gm,long iv,FILE *out)
01661 {
01662 long i,j,stop;
01663 double nom,denom,norrhs,alpha,beta;
01664 double *d,*r,*p;
01665 MPI_Status stat;
01666
01667 if (myrank==0){
01668 d = new double [tnmdof];
01669 r = new double [tnmdof];
01670 p = new double [tnmdof];
01671
01672
01673 if (iv==0){
01674 for (i=0;i<tnmdof;i++)
01675 cglhs[i]=0.0;
01676 }
01677 }
01678 matxvect (top,domproc,gm,cglhs,p,out);
01679
01680 if (myrank==0){
01681 norrhs=0.0; nom=0.0;
01682 for (i=0;i<tnmdof;i++){
01683 norrhs+=cgrhs[i]*cgrhs[i];
01684 r[i]=cgrhs[i]-p[i];
01685 nom+=r[i]*r[i];
01686 d[i]=r[i];
01687 }
01688
01689 if (norrhs<zero){
01690 fprintf (stderr,"\n\n norm of right hand side in conjugate gradient method is smaller than %e",zero);
01691 fprintf (stderr,"\n see file %s, line %d.\n",__FILE__,__LINE__);
01692 aerrcgdpfeti=norrhs; anicgdpfeti=0;
01693 return;
01694 }
01695 }
01696
01697
01698
01699 stop=0;
01700 for (i=0;i<nicgdpfeti;i++){
01701
01702
01703
01704 matxvect (top,domproc,gm,d,p,out);
01705
01706 if (myrank==0){
01707 denom = ss (d,p,tnmdof);
01708 if (fabs(denom)<1.0e-25){
01709 fprintf (stdout,"\n there is zero denominator in alpha computation in conjugate gradient method (%s, line %d)\n",__FILE__,__LINE__);
01710 stop=1;
01711 }
01712
01713 alpha = nom/denom;
01714
01715
01716 for (j=0;j<tnmdof;j++){
01717 cglhs[j]+=alpha*d[j];
01718 r[j]-=alpha*p[j];
01719 }
01720
01721 denom=nom;
01722
01723 nom = ss (r,r,tnmdof);
01724
01725 fprintf (out,"\n iteration %ld norres/norrhs %e",i,sqrt(nom/norrhs));
01726 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,sqrt(nom/norrhs));
01727
01728
01729 if (sqrt(nom/norrhs)<errcgdpfeti) stop=1;
01730
01731
01732
01733 beta = nom/denom;
01734
01735
01736 for (j=0;j<tnmdof;j++){
01737 d[j]=beta*d[j]+r[j];
01738 }
01739
01740 for (j=1;j<nproc;j++){
01741 MPI_Send(&stop,1,MPI_LONG,j,myrank,MPI_COMM_WORLD);
01742 }
01743 }
01744 else{
01745 MPI_Recv(&stop,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01746 }
01747 MPI_Barrier (MPI_COMM_WORLD);
01748
01749 if (stop==1) break;
01750 }
01751
01752 anicgdpfeti=i; aerrcgdpfeti=nom;
01753
01754 if (myrank==0){
01755 delete [] p; delete [] r; delete [] d;
01756 }
01757
01758 }
01759
01760
01761 void dpfeti::corner_displ (gtopology *top,long *domproc,gmatrix *gm,FILE *out)
01762 {
01763 long i,j,gndofe,*cn;
01764 double *utc,*ur,*vr,*buff;
01765 MPI_Status stat;
01766
01767
01768
01769 ur = new double [ndof];
01770 vr = new double [ndof];
01771
01772 buff = new double [maxnbdof];
01773
01774 if (myrank==0){
01775 for (i=1;i<nproc;i++){
01776 extract_from_global_vector (buff,cglhs,domproc[i]);
01777 MPI_Send(buff,maxnbdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01778 }
01779 extract_from_global_vector (buff,cglhs,domproc[0]);
01780 }
01781 else{
01782 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01783 }
01784 MPI_Barrier (MPI_COMM_WORLD);
01785
01786
01787 nullv (ur,nidof);
01788 put_into_local_vector (buff,ur);
01789
01790
01791 gm->condense (top,buff,buff,vr,ur,ncdof,3,out);
01792
01793 delete [] buff;
01794 buff = new double [maxncdof];
01795
01796
01797 mtxv (krc,vr,buff,nidof,ncdof);
01798
01799 if (myrank==0){
01800 utc = new double [tncdof];
01801 nullv (utc,tncdof);
01802
01803
01804 gndofe = ptop->give_ndofe (domproc[0]);
01805 cn = new long [gndofe];
01806 ptop->give_code_numbers (domproc[0],cn);
01807 locglob (utc,buff,cn,gndofe);
01808 delete [] cn;
01809
01810 for (i=1;i<nproc;i++){
01811 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01812
01813 j=domproc[stat.MPI_TAG];
01814 gndofe = ptop->give_ndofe (j);
01815 cn = new long [gndofe];
01816 ptop->give_code_numbers (j,cn);
01817 locglob (utc,buff,cn,gndofe);
01818 delete [] cn;
01819 }
01820
01821
01822
01823
01824
01825 for (i=0;i<tncdof;i++){
01826 utc[i]=br[i]+utc[i];
01827 }
01828
01829
01830 arr->back_substitution (displ,utc);
01831
01832 delete [] utc;
01833 }
01834 else{
01835 MPI_Send(buff,maxncdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
01836 }
01837 MPI_Barrier (MPI_COMM_WORLD);
01838 delete [] buff;
01839 delete [] ur;
01840 delete [] vr;
01841 }
01842
01843 void dpfeti::compute_displ (gtopology *top,gmatrix *gm,long *domproc,
01844 double *subdispl,double *rhs,FILE *out)
01845 {
01846 long i,j,gndofe,*cn;
01847 double *ur,*vr,*buff;
01848 MPI_Status stat;
01849
01850 ur = new double [nidof];
01851 vr = new double [nidof];
01852
01853 buff = new double [maxncdof];
01854 if (myrank==0){
01855 for (i=1;i<nproc;i++){
01856
01857 nullv (buff,maxncdof);
01858 gndofe = ptop->give_ndofe (domproc[i]);
01859 cn = new long [gndofe];
01860 ptop->give_code_numbers (domproc[i],cn);
01861 globloc (displ,buff,cn,gndofe);
01862 delete [] cn;
01863
01864 MPI_Send(buff,maxncdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01865 }
01866
01867 nullv (buff,maxncdof);
01868 gndofe = ptop->give_ndofe (domproc[0]);
01869 cn = new long [gndofe];
01870 ptop->give_code_numbers (domproc[0],cn);
01871 globloc (displ,buff,cn,gndofe);
01872 delete [] cn;
01873
01874 }
01875 else{
01876 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01877 }
01878 MPI_Barrier (MPI_COMM_WORLD);
01879
01880 j=nidof;
01881 for (i=0;i<ncdof;i++){
01882 subdispl[j]=buff[i]; j++;
01883 }
01884
01885
01886 mxv (krc,buff,ur,nidof,ncdof);
01887 delete [] buff;
01888
01889
01890 buff = new double [maxnbdof];
01891
01892 if (myrank==0){
01893 for (i=1;i<nproc;i++){
01894 extract_from_global_vector (buff,cglhs,domproc[i]);
01895 MPI_Send(buff,maxnbdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01896 }
01897 extract_from_global_vector (buff,cglhs,domproc[0]);
01898 }
01899 else{
01900 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01901 }
01902 MPI_Barrier (MPI_COMM_WORLD);
01903
01904
01905 nullv (vr,nidof);
01906 put_into_local_vector (buff,vr);
01907 delete [] buff;
01908
01909 for (i=0;i<nidof;i++){
01910 rhs[i]=rhs[i]-ur[i]-vr[i];
01911 }
01912
01913 gm->condense (top,buff,buff,subdispl,rhs,ncdof,3,out);
01914
01915 delete [] ur;
01916 delete [] vr;
01917 }
01918
01919
01920
01921
01922
01923
01924
01925 void dpfeti::solve_system (gtopology *top,gmatrix *gm,long *domproc,
01926 double *lhs,double *rhs,FILE *out,long mespr)
01927 {
01928 long i,j;
01929 double *condmat,*condvect,*rrhs,*crhs;
01930 time_t t1,t2,t3,t4,t5,t6,t7,t8;
01931
01932
01933 t1 = time (NULL);
01934
01935 if (myrank==0){
01936 condmat = new double [maxncdof*maxncdof];
01937 memset (condmat,0,maxncdof*maxncdof*sizeof(double));
01938 condvect = new double [maxncdof];
01939 memset (condvect,0,maxncdof*sizeof(double));
01940 }
01941 else{
01942 condmat = new double [ncdof*ncdof];
01943 memset (condmat,0,ncdof*ncdof*sizeof(double));
01944 condvect = new double [ncdof];
01945 memset (condvect,0,ncdof*sizeof(double));
01946 }
01947
01948 if (myrank==0){
01949
01950 fprintf (out,"\n\n tncdof %ld",tncdof);
01951
01952 cglhs = new double [tnmdof];
01953 cgrhs = new double [tnmdof];
01954 displ = new double [tncdof];
01955 }
01956
01957
01958 krc = new double [nidof*ncdof];
01959 for (i=0;i<nidof*ncdof;i++){
01960 krc[i]=0.0;
01961 }
01962
01963
01964 gm->a12block (top,krc,ncdof,out);
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976 crhs = new double [ncdof];
01977
01978 rrhs = new double [ndof];
01979
01980
01981 for (i=0;i<nidof;i++){
01982 rrhs[i]=rhs[i];
01983 }
01984 j=0;
01985 for (i=nidof;i<ndof;i++){
01986 crhs[j]=rhs[i]; j++;
01987 }
01988
01989 t2 = time (NULL);
01990
01991
01992 gm->condense (top,condmat,condvect,lhs,rhs,ncdof,1,out);
01993
01994 t3 = time (NULL);
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010 arrmatrix (condmat,domproc,out);
02011
02012
02013 t4 = time (NULL);
02014
02015
02016 vectors_br_bm (top,gm,condmat,condvect,rrhs,domproc,out);
02017
02018
02019
02020 rhs_dpfeti (top,domproc,gm,out);
02021
02022 t5 = time (NULL);
02023
02024
02025
02026 cg (top,domproc,gm,0,out);
02027
02028
02029 t6 = time (NULL);
02030
02031 if (myrank==0){
02032 fprintf (out,"\n\n L A G R A N G E M U L T I P L I E R S \n");
02033 for (i=0;i<tnmdof;i++){
02034 fprintf (out,"\n %7ld %30.15f",i,cglhs[i]);
02035 }
02036 }
02037
02038
02039 corner_displ (top,domproc,gm,out);
02040
02041 t7 = time (NULL);
02042
02043 if (myrank==0){
02044 fprintf (out,"\n\n C O R N E R D I S P L A C E M E N T S \n");
02045 fprintf (out,"\n\n tncdof %ld",tncdof);
02046 for (i=0;i<tncdof;i++){
02047 fprintf (out,"\n %7ld %30.15f",i,displ[i]);
02048 }
02049 fprintf (out,"\n\n konec tisku\n");
02050 }
02051
02052
02053 compute_displ (top,gm,domproc,lhs,rrhs,out);
02054
02055 t8 = time (NULL);
02056
02057 fprintf (out,"\n\n A_rc block assembling %ld",t2-t1);
02058 fprintf (out,"\n\n time of condensation %ld",t3-t2);
02059 fprintf (out,"\n\n A_rr assembling %ld",t4-t3);
02060 fprintf (out,"\n\n b_R and b_M assembling %ld",t5-t4);
02061 fprintf (out,"\n\n solution of reduced system %ld",t6-t5);
02062 fprintf (out,"\n\n corner displacements calcul. %ld",t7-t6);
02063 fprintf (out,"\n\n all displacements calcul. %ld",t8-t7);
02064
02065
02066 MPI_Status stat;
02067 long *buff,buffsize;
02068
02069 buffsize=7;
02070 buff = new long [buffsize];
02071
02072 buff[0]=t2-t1;
02073 buff[1]=t3-t2;
02074 buff[2]=t4-t3;
02075 buff[3]=t5-t4;
02076 buff[4]=t6-t5;
02077 buff[5]=t7-t6;
02078 buff[6]=t8-t7;
02079
02080 if (myrank==0){
02081
02082 long *fact,*cp,*arrasem,*tottime;
02083 fact=new long [nproc];
02084 cp=new long [nproc];
02085 arrasem = new long [nproc];
02086 tottime=new long [nproc];
02087
02088 fprintf (out,"\n\n\n\n\n");
02089 j=domproc[myrank];
02090 fprintf (out,"\n\n\n Domain %ld",j);
02091 fprintf (out,"\n A_rc block assembling %ld",buff[0]);
02092 fprintf (out,"\n time of condensation %ld",buff[1]);
02093 fprintf (out,"\n A_rr assembling %ld",buff[2]);
02094 fprintf (out,"\n b_R and b_M assembling %ld",buff[3]);
02095 fprintf (out,"\n solution of reduced system %ld",buff[4]);
02096 fprintf (out,"\n corner displacements calcul. %ld",buff[5]);
02097 fprintf (out,"\n all displacements calcul. %ld",buff[6]);
02098
02099
02100 j=domproc[myrank];
02101 fact[j]=buff[1];
02102 arrasem[j]=buff[2];
02103 cp[j]=buff[4];
02104 tottime[j]=buff[0]+buff[1]+buff[2]+buff[3]+buff[4]+buff[5]+buff[6];
02105
02106 fprintf (stdout,"\n\n\n Domain %ld",j);
02107 fprintf (stdout,"\n A_rc block assembling %ld",buff[0]);
02108 fprintf (stdout,"\n time of condensation %ld",buff[1]);
02109 fprintf (stdout,"\n A_rr assembling %ld",buff[2]);
02110 fprintf (stdout,"\n b_R and b_M assembling %ld",buff[3]);
02111 fprintf (stdout,"\n solution of reduced system %ld",buff[4]);
02112 fprintf (stdout,"\n corner displacements calcul. %ld",buff[5]);
02113 fprintf (stdout,"\n all displacements calcul. %ld",buff[6]);
02114
02115
02116 for (i=1;i<nproc;i++){
02117 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
02118 j=domproc[stat.MPI_TAG];
02119
02120
02121 fact[j]=buff[1];
02122 arrasem[j]=buff[2];
02123 cp[j]=buff[4];
02124 tottime[j]=buff[0]+buff[1]+buff[2]+buff[3]+buff[4]+buff[5]+buff[6];
02125
02126
02127 fprintf (out,"\n\n\n Domain %ld",j);
02128 fprintf (out,"\n A_rc block assembling %ld",buff[0]);
02129 fprintf (out,"\n time of condensation %ld",buff[1]);
02130 fprintf (out,"\n A_rr assembling %ld",buff[2]);
02131 fprintf (out,"\n b_R and b_M assembling %ld",buff[3]);
02132 fprintf (out,"\n solution of reduced system %ld",buff[4]);
02133 fprintf (out,"\n corner displacements calcul. %ld",buff[5]);
02134 fprintf (out,"\n all displacements calcul. %ld",buff[6]);
02135
02136
02137 fprintf (stdout,"\n\n\n Domain %ld",j);
02138 fprintf (stdout,"\n A_rc block assembling %ld",buff[0]);
02139 fprintf (stdout,"\n time of condensation %ld",buff[1]);
02140 fprintf (stdout,"\n A_rr assembling %ld",buff[2]);
02141 fprintf (stdout,"\n b_R and b_M assembling %ld",buff[3]);
02142 fprintf (stdout,"\n solution of reduced system %ld",buff[4]);
02143 fprintf (stdout,"\n corner displacements calcul. %ld",buff[5]);
02144 fprintf (stdout,"\n all displacements calcul. %ld",buff[6]);
02145
02146
02147 }
02148
02149 long minfact,maxfact,mincp,maxcp,minarr,maxarr,mintot,maxtot;
02150 maxfact=0;
02151 maxcp=0;
02152 maxarr=0;
02153 maxtot=0;
02154 mincp=LONG_MAX;
02155 minarr=LONG_MAX;
02156 minfact=LONG_MAX;
02157 mintot=LONG_MAX;
02158 for (i=0;i<nproc;i++){
02159 if (maxfact<fact[i])
02160 maxfact=fact[i];
02161 if (minfact>fact[i])
02162 minfact=fact[i];
02163
02164 if (mincp>cp[i])
02165 mincp=cp[i];
02166 if (maxcp<cp[i])
02167 maxcp=cp[i];
02168
02169 if (minarr>arrasem[i])
02170 minarr=arrasem[i];
02171 if (maxarr<arrasem[i])
02172 maxarr=arrasem[i];
02173
02174 if (mintot>tottime[i])
02175 mintot=tottime[i];
02176 if (maxtot<tottime[i])
02177 maxtot=tottime[i];
02178 }
02179
02180 fprintf (out,"\n maximum factorization time %ld",maxfact);
02181 fprintf (out,"\n minimum factorization time %ld",minfact);
02182 fprintf (out,"\n maximum coarse problem time %ld",maxcp);
02183 fprintf (out,"\n minimum coarse problem time %ld",mincp);
02184 fprintf (out,"\n maximum A_rr time %ld",maxarr);
02185 fprintf (out,"\n minimum A_rr time %ld",minarr);
02186 fprintf (out,"\n maximum total time %ld",maxtot);
02187 fprintf (out,"\n minimum total time %ld",mintot);
02188
02189 delete [] fact;
02190 delete [] cp;
02191
02192 }
02193 else{
02194 MPI_Send(buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
02195 }
02196 MPI_Barrier (MPI_COMM_WORLD);
02197 delete [] buff;
02198
02199
02200 delete [] condmat;
02201 delete [] condvect;
02202 delete [] rrhs;
02203 delete [] crhs;
02204
02205 }
02206