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