00001 #include "dpfeti_new.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <time.h>
00005
00006 dpfeti::dpfeti(int np,int mr,int nd)
00007 {
00008 nproc=np; myrank=mr; ndom=nd;
00009
00010 }
00011
00012 dpfeti::~dpfeti()
00013 {
00014 delete [] cgrhs;
00015 delete [] cglhs;
00016 delete [] displ;
00017 delete [] br;
00018 delete [] bm;
00019
00020 delete [] krc;
00021
00022 for (i=0;i<nproc;i++){
00023 delete [] loccn[i];
00024 delete [] globcn[i];
00025 }
00026 delete [] loccn;
00027 delete [] globcn;
00028
00029 delete [] nlbdof;
00030 delete [] ncdofdom;
00031
00032 }
00033
00034
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 dpfeti::dpfetiordering (gtopology *top,long *ltg)
00059 {
00060 long i,j,k,ndofn;
00061
00062
00063 ndof=1;
00064 for (i=0;i<top->nn;i++){
00065 if (ltg[i]>=-1){
00066
00067 ndofn=top->give_ndofn (i);
00068 for (j=0;j<ndofn;j++){
00069 k=top->give_dof(i,j);
00070 if (k<0) continue;
00071 if (k==0) continue;
00072 if (k>0){
00073 top->save_dof (i,j,ndof); ndof++;
00074 }
00075 }
00076 }
00077 }
00078 nidof=ndof-1;
00079
00080 for (i=0;i<top->nn;i++){
00081 if (ltg[i]<-1){
00082
00083 ndofn=top->give_ndofn (i);
00084 for (j=0;j<ndofn;j++){
00085 k=top->give_dof (i,j);
00086 if (k<0) continue;
00087 if (k==0) continue;
00088 if (k>0){
00089 top->save_dof (i,j,ndof); ndof++;
00090 }
00091 }
00092 }
00093 }
00094 ndof--;
00095
00096
00097 for (i=0;i<top->ne;i++){
00098 if (top->give_cne (i)==1){
00099 for (j=0;j<top->give_ndofe(i);j++){
00100 if (top->give_dof (i,j)>ndof) ndof=top->give_dof (i,j);
00101 }
00102 }
00103 }
00104
00105 ncdof=ndof-nidof;
00106 }
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 void dpfeti::globcnnum_dpfeti (gtopology *top,long *ltg,long *domproc,FILE *out)
00117 {
00118 long i,j,k,l,m,n,ndofn,stop,nbn,ncn,maxnbn,maxncn,tnbn,tncn;
00119 long *nibn,*lgnbn,*llnbn,*lgncn,*llncn,***lcnbn,***lcncn,**gnbndom,**gncndom;
00120 long **cncn,**cnbn,*nbndom,*ncndom,**dominc,**supelemcn;
00121 long buffsize,*buff;
00122 MPI_Status stat;
00123
00124
00125 ndofn = top->give_ndofn (0);
00126
00127
00128
00129 nbn=0; ncn=0;
00130 for (i=0;i<top->nn;i++){
00131 if (ltg[i]>-1) nbn++;
00132 if (ltg[i]<-1) ncn++;
00133 }
00134
00135
00136
00137 lgnbn = new long [nbn];
00138 lgncn = new long [ncn];
00139
00140
00141 llnbn = new long [nbn];
00142 llncn = new long [ncn];
00143
00144 j=0; k=0;
00145 for (i=0;i<top->nn;i++){
00146 if (ltg[i]>-1){
00147 lgnbn[j]=ltg[i];
00148 llnbn[j]=i;
00149 j++;
00150 }
00151 if (ltg[i]<-1){
00152 lgncn[k]=0-ltg[i]-2;
00153 llncn[k]=i;
00154 k++;
00155 }
00156 }
00157
00158 buff = new long [2];
00159 buff[0]=nbn;
00160 buff[1]=ncn;
00161
00162 if (myrank==0){
00163
00164 nbndom = new long [nproc];
00165
00166 ncndom = new long [nproc];
00167
00168 nbndom[ndom]=nbn;
00169 ncndom[ndom]=ncn;
00170
00171
00172
00173 maxnbn=nbn; maxncn=ncn;
00174 for (i=1;i<nproc;i++){
00175 MPI_Recv (buff,2,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00176
00177 nbndom[stat.MPI_TAG]=buff[0];
00178 ncndom[stat.MPI_TAG]=buff[1];
00179
00180 if (maxnbn<buff[0]) maxnbn=buff[0];
00181 if (maxncn<buff[1]) maxncn=buff[1];
00182 }
00183
00184 buff[0]=maxnbn;
00185 buff[1]=maxncn;
00186 for (i=1;i<nproc;i++){
00187 MPI_Send (buff,2,MPI_LONG,i,ndom,MPI_COMM_WORLD);
00188 }
00189 }
00190 else{
00191 MPI_Send (buff,2,MPI_LONG,0,ndom,MPI_COMM_WORLD);
00192 MPI_Recv (buff,2,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00193 maxnbn=buff[0];
00194 maxncn=buff[1];
00195 }
00196 MPI_Barrier (MPI_COMM_WORLD);
00197 delete [] buff;
00198
00199
00200 buffsize = (maxnbn+maxncn)*(ndofn+1);
00201 buff = new long [buffsize];
00202 k=0;
00203 for (i=0;i<nbn;i++){
00204 buff[k]=lgnbn[i]; k++;
00205 for (j=0;j<ndofn;j++){
00206 buff[k]=top->give_dof (llnbn[i],j);
00207 k++;
00208 }
00209 }
00210 for (i=0;i<ncn;i++){
00211 buff[k]=lgncn[i]; k++;
00212 for (j=0;j<ndofn;j++){
00213 buff[k]=top->give_dof (llncn[i],j);
00214 k++;
00215 }
00216 }
00217
00218 delete [] llnbn;
00219 delete [] llncn;
00220 delete [] lgnbn;
00221 delete [] lgncn;
00222
00223 if (myrank==0){
00224
00225
00226 lcnbn = new long** [nproc];
00227
00228 lcncn = new long** [nproc];
00229
00230 gnbndom = new long* [nproc];
00231
00232 gncndom = new long* [nproc];
00233 for (i=0;i<nproc;i++){
00234 gnbndom[i] = new long [nbndom[i]];
00235 gncndom[i] = new long [ncndom[i]];
00236
00237 lcnbn[i] = new long* [nbndom[i]];
00238 lcncn[i] = new long* [ncndom[i]];
00239 for (j=0;j<nbndom[i];j++){
00240 lcnbn[i][j] = new long [ndofn];
00241 }
00242 for (j=0;j<ncndom[i];j++){
00243 lcncn[i][j] = new long [ndofn];
00244 }
00245 }
00246
00247
00248 l=0;
00249 for (j=0;j<nbndom[ndom];j++){
00250 gnbndom[ndom][j]=buff[l]; l++;
00251 for (k=0;k<ndofn;k++){
00252 lcnbn[ndom][j][k]=buff[l];
00253 l++;
00254 }
00255 }
00256 for (j=0;j<ncndom[ndom];j++){
00257 gncndom[ndom][j]=buff[l]; l++;
00258 for (k=0;k<ndofn;k++){
00259 lcncn[ndom][j][k]=buff[l];
00260 l++;
00261 }
00262 }
00263
00264
00265 for (i=1;i<nproc;i++){
00266 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00267
00268 l=0;
00269 for (j=0;j<nbndom[stat.MPI_TAG];j++){
00270 gnbndom[stat.MPI_TAG][j]=buff[l]; l++;
00271 for (k=0;k<ndofn;k++){
00272 lcnbn[stat.MPI_TAG][j][k]=buff[l];
00273 l++;
00274 }
00275 }
00276 for (j=0;j<ncndom[stat.MPI_TAG];j++){
00277 gncndom[stat.MPI_TAG][j]=buff[l]; l++;
00278 for (k=0;k<ndofn;k++){
00279 lcncn[stat.MPI_TAG][j][k]=buff[l];
00280 l++;
00281 }
00282 }
00283
00284 }
00285 }
00286 else{
00287 MPI_Send (buff,buffsize,MPI_LONG,0,ndom,MPI_COMM_WORLD);
00288 }
00289 MPI_Barrier (MPI_COMM_WORLD);
00290 delete [] buff;
00291
00292
00293
00294 if (myrank==0){
00295
00296
00297
00298 tnbn = 0;
00299 for (i=0;i<nproc;i++){
00300 for (j=0;j<nbndom[i];j++){
00301 if (gnbndom[i][j]>tnbn) tnbn=gnbndom[i][j];
00302 }
00303 }
00304 tnbn++;
00305
00306 tncn = 0;
00307 for (i=0;i<nproc;i++){
00308 for (j=0;j<ncndom[i];j++){
00309 if (gncndom[i][j]>tncn) tncn=gncndom[i][j];
00310 }
00311 }
00312 tncn++;
00313
00314
00315 cncn = new long* [tncn];
00316 for (i=0;i<tncn;i++){
00317 cncn[i] = new long [ndofn];
00318 for (j=0;j<ndofn;j++){
00319 cncn[i][j]=0;
00320 }
00321 }
00322
00323
00324 for (i=0;i<nproc;i++){
00325 for (j=0;j<ncndom[i];j++){
00326 for (k=0;k<ndofn;k++){
00327 cncn[gncndom[i][j]][k]=lcncn[i][j][k];
00328 }
00329 }
00330 }
00331
00332
00333 tncdof=1;
00334 for (i=0;i<tncn;i++){
00335 for (j=0;j<ndofn;j++){
00336 if (cncn[i][j]>0){
00337 cncn[i][j]=tncdof;
00338 tncdof++;
00339 }
00340 }
00341 }
00342 tncdof--;
00343
00344
00345 ncdofdom = new long [nproc];
00346 for (i=0;i<nproc;i++){
00347 ncdofdom[i]=0;
00348 }
00349
00350
00351 for (i=0;i<nproc;i++){
00352 for (j=0;j<ncndom[i];j++){
00353 for (k=0;k<ndofn;k++){
00354 if (cncn[gncndom[i][j]][k]>0) ncdofdom[i]++;
00355 }
00356 }
00357 }
00358
00359 maxncdof=0;
00360 for (i=0;i<nproc;i++){
00361 if (maxncdof<ncdofdom[i]) maxncdof=ncdofdom[i];
00362 }
00363
00364
00365
00366 supelemcn = new long* [nproc];
00367 for (i=0;i<nproc;i++){
00368 supelemcn[i] = new long [ncdofdom[i]];
00369 for (j=0;j<ncdofdom[i];j++){
00370 supelemcn[i][j]=0;
00371 }
00372 }
00373
00374
00375 for (i=0;i<nproc;i++){
00376 l=0;
00377 for (j=0;j<ncndom[i];j++){
00378 for (k=0;k<ndofn;k++){
00379 m=cncn[gncndom[i][j]][k];
00380 if (m>0){
00381 supelemcn[i][l]=m;
00382 l++;
00383 }
00384 }
00385 }
00386 }
00387
00388
00389 ptop = new gtopology;
00390 ptop->initiate (supelemcn,ncdofdom,nproc);
00391
00392 for (i=0;i<nproc;i++){
00393 delete [] gncndom[i];
00394 }
00395 delete [] gncndom;
00396
00397 for (i=0;i<nproc;i++){
00398 delete [] supelemcn[i];
00399 }
00400 delete [] supelemcn;
00401
00402 for (i=0;i<tncn;i++){
00403 delete [] cncn[i];
00404 }
00405 delete [] cncn;
00406
00407
00408
00409
00410
00411
00412 nibn = new long [tnbn];
00413 for (i=0;i<tnbn;i++){
00414 nibn[i]=0;
00415 }
00416
00417
00418 for (i=0;i<nproc;i++){
00419 for (j=0;j<nbndom[i];j++){
00420 nibn[gnbndom[i][j]]++;
00421 }
00422 }
00423
00424
00425 dominc = new long* [tnbn];
00426 for (i=0;i<tnbn;i++){
00427 dominc[i] = new long [nibn[i]];
00428 for (j=0;j<nibn[i];j++){
00429 dominc[i][j]=-1;
00430 }
00431 }
00432
00433 for (i=0;i<tnbn;i++){
00434 nibn[i]=0;
00435 }
00436
00437
00438 for (i=0;i<nproc;i++){
00439 for (j=0;j<nbndom[i];j++){
00440 dominc[gnbndom[i][j]][nibn[gnbndom[i][j]]]=i;
00441 nibn[gnbndom[i][j]]++;
00442 }
00443 }
00444
00445
00446
00447 cnbn = new long** [tnbn];
00448 for (i=0;i<tnbn;i++){
00449 cnbn[i] = new long* [nibn[i]-1];
00450 for (j=0;j<nibn[i]-1;j++){
00451 cnbn[i][j] = new long [ndofn];
00452 for (k=0;k<ndofn;k++){
00453 cnbn[i][j][k]=0;
00454 }
00455 }
00456 }
00457
00458
00459 for (i=0;i<nproc;i++){
00460 for (j=0;j<nbndom[i];j++){
00461 m=gnbndom[i][j];
00462 for (k=0;k<nibn[m]-1;k++){
00463 for (l=0;l<ndofn;l++){
00464 cnbn[m][k][l]=lcnbn[i][j][l];
00465 }
00466 }
00467 }
00468 }
00469
00470
00471 tnmdof=1;
00472 for (i=0;i<tnbn;i++){
00473 for (j=0;j<nibn[i]-1;j++){
00474 for (k=0;k<ndofn;k++){
00475 if (cnbn[i][j][k]>0){
00476 cnbn[i][j][k]=tnmdof;
00477 tnmdof++;
00478 }
00479 }
00480 }
00481 }
00482 tnmdof--;
00483
00484
00485 nlbdofdom = new long [nproc];
00486 for (i=0;i<nproc;i++){
00487 nlbdofdom[i]=0;
00488 }
00489 for (i=0;i<nproc;i++){
00490 for (j=0;j<nbndom[i];j++){
00491 m=gnbndom[i][j];
00492 if (dominc[m][0]==i){
00493 for (k=0;k<nibn[m]-1;k++){
00494 for (l=0;l<ndofn;l++){
00495 if (cnbn[m][k][l]>0) nlbdofdom[i]++;
00496 }
00497 }
00498 }
00499 else{
00500 for (l=0;l<ndofn;l++){
00501 if (cnbn[m][0][l]>0) nlbdofdom[i]++;
00502 }
00503 }
00504 }
00505 }
00506
00507
00508
00509 loccn = new long* [nproc];
00510 globcn = new long* [nproc];
00511 maxnbdof=0;
00512 for (i=0;i<nproc;i++){
00513 loccn[i] = new long [nlbdofdom[i]];
00514 globcn[i] = new long [nlbdofdom[i]];
00515 if (maxnbdof<nlbdofdom[i]) maxnbdof=nlbdofdom[i];
00516 }
00517
00518 for (i=0;i<nproc;i++){
00519 n=0;
00520 for (j=0;j<nbndom[i];j++){
00521 m=gnbndom[i][j];
00522 if (dominc[m][0]==i){
00523 for (k=0;k<nibn[m]-1;k++){
00524 for (l=0;l<ndofn;l++){
00525 if (cnbn[m][k][l]>0){
00526 globcn[i][n]=0-cnbn[m][k][l];
00527 loccn[i][n]=lcnbn[i][j][l];
00528 n++;
00529 }
00530 }
00531 }
00532 }
00533 else{
00534 stop=0;
00535 for (k=0;k<nibn[m]-1;k++){
00536 for (l=0;l<ndofn;l++){
00537 if (cnbn[m][k][l]==0){
00538 cnbn[m][k][k]=-3;
00539 stop=1;
00540 }
00541 if (cnbn[m][k][l]>0){
00542 globcn[i][n]=cnbn[m][k][l];
00543 loccn[i][n]=lcnbn[i][j][l];
00544 n++;
00545 cnbn[m][k][k]=-3;
00546 stop=1;
00547 }
00548 }
00549 if (stop==1) break;
00550 }
00551 }
00552 }
00553 }
00554
00555 for (i=0;i<nproc;i++){
00556 for (j=0;j<nbndom[i];j++){
00557 delete [] lcnbn[i][j];
00558 }
00559 delete [] lcnbn[i];
00560 }
00561 delete [] lcnbn;
00562
00563 for (i=0;i<nproc;i++){
00564 for (j=0;j<ncndom[i];j++){
00565 delete [] lcncn[i][j];
00566 }
00567 delete [] lcncn[i];
00568 }
00569 delete [] lcncn;
00570
00571 for (i=0;i<nproc;i++){
00572 delete [] nbndom[i];
00573 }
00574 delete [] nbndom;
00575
00576 for (i=0;i<nproc;i++){
00577 delete [] ncndom[i];
00578 }
00579 delete [] ncndom;
00580
00581 for (i=0;i<nproc;i++){
00582 delete [] gnbndom[i];
00583 }
00584 delete [] gnbndom;
00585
00586 for (i=0;i<tnbn;i++){
00587 delete [] dominc[i];
00588 }
00589 delete [] dominc;
00590
00591 for (i=0;i<tnbn;i++){
00592 for (j=0;j<nibn[i]-1;j++){
00593 delete [] cnbn[i][j];
00594 }
00595 delete [] cnbn[i];
00596 }
00597 delete [] cnbn;
00598
00599 delete [] nibn;
00600
00601 }
00602
00603 buff = new long [3];
00604 if (myrank==0){
00605 buff[0]=maxncdof;
00606 buff[1]=maxnbdof;
00607
00608 for (i=0;i<nproc;i++){
00609 if (domproc[i]==0) continue;
00610 buff[2]=nlbdofdom[i];
00611 MPI_Send (buff,3,MPI_LONG,domproc[i],ndom,MPI_COMM_WORLD);
00612 }
00613 nlbdof=nlbdofdom[domproc[nproc]];
00614 }
00615 else{
00616 MPI_Recv (buff,3,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00617 maxncdof=buff[0];
00618 maxnbdof=buff[1];
00619 nlbdof=buff[2];
00620 }
00621 MPI_Barrier (MPI_COMM_WORLD);
00622 delete [] buff;
00623
00624 localcn = new long [nlbdof];
00625 globalcn = new long [nlbdof];
00626
00627 buffsize=maxnbdof;
00628 buff = new long [buffsize];
00629 if (myrank==0){
00630 for (i=0;i<nproc;i++){
00631 if (domproc[i]==0) continue;
00632 for (j=0;j<nlbdofdom[i];j++){
00633 buff[j]=loccn[i][j];
00634 }
00635 MPI_Send (buff,buffsize,MPI_LONG,domproc[i],ndom,MPI_COMM_WORLD);
00636 }
00637 }
00638 else{
00639 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00640 for (j=0;j<nlbdof;j++){
00641 localcn[j]=buff[j];
00642 }
00643 }
00644 MPI_Barrier (MPI_COMM_WORLD);
00645 delete [] buff;
00646
00647 }
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659 void dpfeti::extract_from_local_vector (double *ev,double *lv)
00660 {
00661 long i,lcn;
00662
00663 for (i=0;i<nlbdof;i++){
00664 lcn=localcn[i];
00665 if (lcn==0) continue;
00666 ev[i]=lv[lcn-1];
00667 }
00668
00669 }
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680 void dpfeti::put_into_local_vector (double *ev,double *lv)
00681 {
00682 long i,lcn;
00683
00684 for (i=0;i<nlbdof;i++){
00685 lcn=localcn[i];
00686 if (lcn==0) continue;
00687 lv[lcn-1]+=ev[i];
00688 }
00689
00690 }
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703 void dpfeti::extract_from_global_vector (double *ev,double *gv,long nsub)
00704 {
00705 long i,gcn;
00706
00707 for (i=0;i<nlbdofdom[nsub];i++){
00708 gcn=globcn[nsub][i];
00709 if (gcn>0)
00710 ev[i]=gv[gcn-1];
00711 if (gcn<0)
00712 ev[i]=0-gv[0-gcn-1];
00713 }
00714
00715 }
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727 void dpfeti::put_into_global_vector (double *ev,double *gv,long nsub)
00728 {
00729 long i,gcn;
00730
00731 for (i=0;i<nlbdofdom[nsub];i++){
00732 gcn=globcn[nsub][i];
00733 if (gcn>0)
00734 gv[gcn-1]+=ev[i];
00735 if (gcn<0)
00736 gv[gcn-1]-=ev[i];
00737 }
00738
00739 }
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755 void dpfeti::arrmatrix (double *condmat)
00756 {
00757 long i,j,k,l,ii,gndofe,*cn;
00758 double *buff;
00759 MPI_Status stat;
00760
00761 buff = new double [maxncdof*maxncdof];
00762 ii=0;
00763 for (i=0;i<ncdof;i++){
00764 for (j=0;j<ncdof;j++){
00765 buff[ii]=condmat[i*ncdof+j]; ii++;
00766 }
00767 }
00768
00769 if (myrank==0){
00770 arr = new gmatrix;
00771
00772 arr->ts = rsmstor;
00773 arr->tlinsol = tlinsol;
00774
00775 arr->zero = zero;
00776 arr->limit = limit;
00777
00778 arr->initiate (ptop,tncdof,1);
00779 arr->prepmat (0.0,1);
00780
00781 gndofe = ptop->give_ndofe (ndom);
00782 cn = new long [gndofe];
00783 ptop->give_code_numbers (ndom,cn);
00784 arr->localized (condmat,cn,gndofe);
00785 delete [] cn;
00786
00787 for (i=1;i<nproc;i++){
00788 MPI_Recv(buff,maxncdof*maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00789 l=stat.MPI_TAG;
00790
00791 ii=0;
00792 for (j=0;j<ncdofdom[l];j++){
00793 for (k=0;k<ncdofdom[l];k++){
00794 condmat[j*ncdofdom[l]+k]=buff[ii]; ii++;
00795 }
00796 }
00797
00798 gndofe = ptop->give_ndofe (l);
00799 cn = new long [gndofe];
00800 ptop->give_code_numbers (l,cn);
00801 arr->localized (condmat,cn,gndofe);
00802 delete [] cn;
00803
00804 }
00805
00806 arr->decompose_matrix ();
00807
00808 }
00809 else{
00810 MPI_Send(buff,maxncdof*maxncdof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
00811 }
00812 MPI_Barrier (MPI_COMM_WORLD);
00813 delete [] buff;
00814 }
00815
00816
00817 void dpfeti::vectors_br_bm (gtopology *top,gmatrix *gm,
00818 double *condmat,double *condvect,double *rhs)
00819 {
00820 long i,j,l,gndofe,*cn;
00821 double *lhs,*buff,*zaloha;
00822 MPI_Status stat;
00823
00824
00825
00826
00827 buff = new double [maxncdof];
00828 for (i=0;i<ncdof;i++){
00829 buff[i]=condvect[i];
00830 }
00831
00832 if (myrank==0){
00833 br = new double [tncdof];
00834 nullvr (br,tncdof);
00835
00836 gndofe = ptop->give_ndofe (ndom);
00837 cn = new long [gndofe];
00838 ptop->give_code_numbers (ndom,cn);
00839 locglob (br,condvect,cn,gndofe);
00840 delete [] cn;
00841
00842 for (i=1;i<nproc;i++){
00843 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00844 l=stat.MPI_TAG;
00845
00846 for (j=0;j<ncdofdom[l];j++){
00847 condvect[j]=buff[j];
00848 }
00849
00850 gndofe = ptop->give_ndofe (l);
00851 cn = new long [gndofe];
00852 ptop->give_code_numbers (l,cn);
00853 locglob (br,condvect,cn,gndofe);
00854 delete [] cn;
00855
00856 }
00857
00858 }
00859 else{
00860 MPI_Send(buff,maxncdof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
00861 }
00862 MPI_Barrier (MPI_COMM_WORLD);
00863 delete [] buff;
00864
00865
00866
00867
00868
00869 buff = new double [maxnbdof];
00870 lhs = new double [nidof];
00871 zaloha = new double [nidof];
00872 for (i=0;i<nidof;i++){
00873 zaloha[i]=rhs[i];
00874 }
00875
00876
00877
00878 gm->condense (condmat,lhs,zaloha,ncdof,3);
00879 delete [] zaloha;
00880
00881
00882
00883 nullvr (buff,maxnbdof);
00884 extract_from_local_vector (buff,lhs);
00885
00886 if (myrank==0){
00887 bm = new double [tnmdof];
00888 nullvr (bm,tnmdof);
00889 put_into_global_vector (buff,bm,ndom);
00890 for (j=1;j<nproc;j++){
00891 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00892 put_into_global_vector (buff,bm,stat.MPI_TAG);
00893 }
00894 }
00895 else{
00896 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
00897 }
00898 MPI_Barrier (MPI_COMM_WORLD);
00899 delete [] buff;
00900 delete [] lhs;
00901
00902 }
00903
00904
00905 void dpfeti::rhs_dpfeti (gtopology *top,long *domproc,gmatrix *gm)
00906 {
00907 long i,j,gndofe,*cn;
00908 double *ur,*wr,*utc,*vtc,*buff;
00909 MPI_Status stat;
00910
00911 ur = new double [nidof];
00912 wr = new double [nidof];
00913 buff = new double [maxncdof];
00914
00915 if (myrank==0){
00916 utc = new double [tncdof];
00917 vtc = new double [tncdof];
00918
00919 for (i=0;i<tncdof;i++){
00920 utc[i]=br[i];
00921 }
00922
00923
00924 arr->back_substitution (vtc,utc);
00925
00926 for (i=0;i<nproc;i++){
00927
00928 if (domproc[i]==0) continue;
00929 nullvr (buff,maxncdof);
00930 gndofe = ptop->give_ndofe (i);
00931 cn = new long [gndofe];
00932 ptop->give_code_numbers (i,cn);
00933 globloc (vtc,buff,cn,gndofe);
00934 delete [] cn;
00935
00936 MPI_Send(buff,maxncdof,MPI_DOUBLE,domproc[i],ndom,MPI_COMM_WORLD);
00937 }
00938
00939 nullvr (buff,maxncdof);
00940 gndofe = ptop->give_ndofe (ndom);
00941 cn = new long [gndofe];
00942 ptop->give_code_numbers (ndom,cn);
00943 globloc (vtc,buff,cn,gndofe);
00944 delete [] cn;
00945
00946 delete [] utc;
00947 delete [] vtc;
00948 }
00949 else{
00950 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00951 }
00952 MPI_Barrier (MPI_COMM_WORLD);
00953
00954
00955
00956 mv (krc,buff,ur,nidof,ncdof);
00957 delete [] buff;
00958
00959
00960 gm->condense (buff,wr,ur,ncdof,3);
00961
00962
00963 buff = new double [maxnbdof];
00964
00965
00966 nullvr (buff,maxnbdof);
00967 extract_from_local_vector (buff,wr);
00968
00969 if (myrank==0){
00970 nullvr (cgrhs,tnmdof);
00971 put_into_global_vector (buff,cgrhs,stat.MPI_TAG);
00972 for (j=1;j<nproc;j++){
00973 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00974 put_into_global_vector (buff,cgrhs,stat.MPI_TAG);
00975 }
00976 }
00977 else{
00978 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
00979 }
00980 MPI_Barrier (MPI_COMM_WORLD);
00981
00982 delete [] buff;
00983 delete [] ur;
00984 delete [] wr;
00985
00986 if (myrank==0){
00987 for (i=0;i<tnmdof;i++){
00988 cgrhs[i]=bm[i]-cgrhs[i];
00989 }
00990 }
00991
00992 }
00993
00994 void dpfeti::matxvect (gtopology *top,long *domproc,
00995 gmatrix *gm,double *input,double *output)
00996 {
00997 long i,j,gndofe,*cn;
00998 double *utc,*vtc,*ur,*vr,*wr,*buff;
00999 MPI_Status stat;
01000
01001 ur = new double [nidof];
01002 vr = new double [nidof];
01003 wr = new double [nidof];
01004
01005 buff = new double [maxnbdof];
01006
01007 if (myrank==0){
01008 for (i=0;i<nproc;i++){
01009 if (domproc[i]==0) continue;
01010 extract_from_global_vector (buff,input,i);
01011 MPI_Send(buff,maxnbdof,MPI_DOUBLE,domproc[i],ndom,MPI_COMM_WORLD);
01012 }
01013 extract_from_global_vector (buff,input,domproc[nproc]);
01014 }
01015 else{
01016 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01017 }
01018 MPI_Barrier (MPI_COMM_WORLD);
01019
01020
01021 nullvr (ur,nidof);
01022 put_into_local_vector (buff,ur);
01023 delete [] buff;
01024
01025
01026 gm->condense (buff,vr,ur,ncdof,3);
01027
01028 buff = new double [maxncdof];
01029
01030
01031 mtv (krc,vr,buff,nidof,ncdof);
01032
01033 if (myrank==0){
01034 utc = new double [tncdof];
01035 vtc = new double [tncdof];
01036 nullvr (utc,tncdof);
01037
01038
01039 gndofe = ptop->give_ndofe (ndom);
01040 cn = new long [gndofe];
01041 ptop->give_code_numbers (ndom,cn);
01042 locglob (utc,buff,cn,gndofe);
01043 delete [] cn;
01044
01045 for (i=1;i<nproc;i++){
01046 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01047
01048 j=stat.MPI_TAG;
01049 gndofe = ptop->give_ndofe (j);
01050 cn = new long [gndofe];
01051 ptop->give_code_numbers (j,cn);
01052 locglob (utc,buff,cn,gndofe);
01053 delete [] cn;
01054 }
01055
01056
01057 arr->back_substitution (vtc,utc);
01058
01059 for (i=0;i<nproc;i++){
01060
01061 if (domproc[i]==0) continue;
01062 nullvr (buff,maxncdof);
01063 gndofe = ptop->give_ndofe (i);
01064 cn = new long [gndofe];
01065 ptop->give_code_numbers (i,cn);
01066 globloc (vtc,buff,cn,gndofe);
01067 delete [] cn;
01068
01069 MPI_Send(buff,maxncdof,MPI_DOUBLE,domproc[i],ndom,MPI_COMM_WORLD);
01070 }
01071
01072 nullvr (buff,maxncdof);
01073 gndofe = ptop->give_ndofe (ndom);
01074 cn = new long [gndofe];
01075 ptop->give_code_numbers (ndom,cn);
01076 globloc (vtc,buff,cn,gndofe);
01077 delete [] cn;
01078
01079 delete [] utc;
01080 delete [] vtc;
01081 }
01082 else{
01083 MPI_Send(buff,maxncdof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01084 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01085 }
01086 MPI_Barrier (MPI_COMM_WORLD);
01087
01088
01089 mv (krc,buff,ur,nidof,ncdof);
01090 delete [] buff;
01091
01092
01093 gm->condense (buff,wr,ur,ncdof,3);
01094
01095
01096
01097
01098
01099 for (i=0;i<nidof;i++){
01100 wr[i]=vr[i]+wr[i];
01101 }
01102
01103 buff = new double [maxnbdof];
01104
01105
01106 extract_from_local_vector (buff,wr);
01107
01108 if (myrank==0){
01109 nullvr (output,tnmdof);
01110 put_into_global_vector (buff,output,ndom);
01111 for (j=1;j<nproc;j++){
01112 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01113 put_into_global_vector (buff,output,stat.MPI_TAG);
01114 }
01115 }
01116 else{
01117 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01118 }
01119 delete [] buff;
01120 delete [] ur;
01121 delete [] vr;
01122 delete [] wr;
01123 }
01124
01125
01126
01127
01128
01129
01130 void dpfeti::cg (gtopology *top,long *domproc,gmatrix *gm,long iv,FILE *out)
01131 {
01132 long i,j,stop;
01133 double nom,denom,norrhs,alpha,beta;
01134 double *d,*r,*p;
01135 MPI_Status stat;
01136
01137 if (myrank==0){
01138 d = new double [tnmdof];
01139 r = new double [tnmdof];
01140 p = new double [tnmdof];
01141
01142
01143 if (iv==0){
01144 for (i=0;i<tnmdof;i++)
01145 cglhs[i]=0.0;
01146 }
01147 }
01148 matxvect (top,domproc,gm,cglhs,p);
01149
01150 if (myrank==0){
01151 norrhs=0.0; nom=0.0;
01152 for (i=0;i<tnmdof;i++){
01153 norrhs+=cgrhs[i]*cgrhs[i];
01154 r[i]=p[i]-cgrhs[i];
01155 nom+=r[i]*r[i];
01156 d[i]=-1.0*r[i];
01157 }
01158
01159 if (norrhs<zero){
01160 fprintf (stderr,"\n\n norm of right hand side in conjugate gradient method is smaller than %e",zero);
01161 fprintf (stderr,"\n see file %s, line %d.\n",__FILE__,__LINE__);
01162 aerrcg=norrhs; anicg=0;
01163 return;
01164 }
01165 }
01166
01167
01168
01169 stop=0;
01170 for (i=0;i<nicg;i++){
01171
01172
01173
01174 matxvect (top,domproc,gm,d,p);
01175
01176 if (myrank==0){
01177 denom = ss (d,p,tnmdof);
01178 if (fabs(denom)<zero){
01179 fprintf (stdout,"\n there is zero denominator in alpha computation in conjugate gradient method (%s, line %d)\n",__FILE__,__LINE__);
01180 stop=1;
01181 }
01182
01183 alpha = nom/denom;
01184
01185
01186 for (j=0;j<tnmdof;j++){
01187 cglhs[j]+=alpha*d[j];
01188 r[j]+=alpha*p[j];
01189 }
01190
01191 denom=nom;
01192
01193 nom = ss (r,r,tnmdof);
01194
01195 fprintf (out,"\n iteration %ld norres/norrhs %e",i,sqrt(nom/norrhs));
01196 fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,sqrt(nom/norrhs));
01197
01198
01199 if (sqrt(nom/norrhs)<errcg) stop=1;
01200
01201
01202
01203 beta = nom/denom;
01204
01205
01206 for (j=0;j<tnmdof;j++){
01207 d[j]=beta*d[j]-r[j];
01208 }
01209
01210 for (j=1;j<nproc;j++){
01211 MPI_Send(&stop,1,MPI_LONG,j,ndom,MPI_COMM_WORLD);
01212 }
01213 }
01214 else{
01215 MPI_Recv(&stop,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01216 }
01217 MPI_Barrier (MPI_COMM_WORLD);
01218
01219 if (stop==1) break;
01220 }
01221
01222 anicg=i; aerrcg=nom;
01223
01224 if (myrank==0){
01225 delete [] p; delete [] r; delete [] d;
01226 }
01227
01228 }
01229
01230
01231 void dpfeti::corner_displ (gtopology *top,long *domproc,gmatrix *gm)
01232 {
01233 long i,j,gndofe,*cn;
01234 double *utc,*ur,*vr,*buff;
01235 MPI_Status stat;
01236
01237 ur = new double [nidof];
01238 vr = new double [nidof];
01239
01240 buff = new double [maxnbdof];
01241
01242 if (myrank==0){
01243 for (i=1;i<nproc;i++){
01244 if (domproc[i]==0) continue;
01245 extract_from_global_vector (buff,cglhs,i);
01246 MPI_Send(buff,maxnbdof,MPI_DOUBLE,i,ndom,MPI_COMM_WORLD);
01247 }
01248 extract_from_global_vector (buff,cglhs,domproc[nproc]);
01249 }
01250 else{
01251 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01252 }
01253 MPI_Barrier (MPI_COMM_WORLD);
01254
01255
01256 nullvr (ur,nidof);
01257 put_into_local_vector (buff,ur);
01258
01259
01260 gm->condense (buff,vr,ur,ncdof,3);
01261
01262 delete [] buff;
01263 buff = new double [maxncdof];
01264
01265
01266 mtv (krc,vr,buff,nidof,ncdof);
01267
01268 if (myrank==0){
01269 utc = new double [tncdof];
01270 nullvr (utc,tncdof);
01271
01272
01273 gndofe = ptop->give_ndofe (ndom);
01274 cn = new long [gndofe];
01275 ptop->give_code_numbers (ndom,cn);
01276 locglob (utc,buff,cn,gndofe);
01277 delete [] cn;
01278
01279 for (i=1;i<nproc;i++){
01280 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01281
01282 j=stat.MPI_TAG;
01283 gndofe = ptop->give_ndofe (j);
01284 cn = new long [gndofe];
01285 ptop->give_code_numbers (j,cn);
01286 locglob (utc,buff,cn,gndofe);
01287 delete [] cn;
01288 }
01289
01290
01291
01292
01293
01294 for (i=0;i<tncdof;i++){
01295 utc[i]=br[i]+utc[i];
01296 }
01297
01298
01299 arr->back_substitution (displ,utc);
01300
01301 delete [] utc;
01302 }
01303 else{
01304 MPI_Send(buff,maxncdof,MPI_DOUBLE,0,ndom,MPI_COMM_WORLD);
01305 }
01306 MPI_Barrier (MPI_COMM_WORLD);
01307 delete [] buff;
01308 delete [] ur;
01309 delete [] vr;
01310 }
01311
01312 void dpfeti::compute_displ (gtopology *top,gmatrix *gm,long *domproc,
01313 double *subdispl,double *rhs)
01314 {
01315 long i,j,gndofe,*cn;
01316 double *ur,*vr,*buff;
01317 MPI_Status stat;
01318
01319 ur = new double [nidof];
01320 vr = new double [nidof];
01321
01322 buff = new double [maxncdof];
01323 if (myrank==0){
01324 for (i=0;i<nproc;i++){
01325
01326 if (domproc[i]==0) continue;
01327 nullvr (buff,maxncdof);
01328 gndofe = ptop->give_ndofe (i);
01329 cn = new long [gndofe];
01330 ptop->give_code_numbers (i,cn);
01331 globloc (displ,buff,cn,gndofe);
01332 delete [] cn;
01333
01334 MPI_Send(buff,maxncdof,MPI_DOUBLE,domproc[i],ndom,MPI_COMM_WORLD);
01335 }
01336
01337 nullvr (buff,maxncdof);
01338 gndofe = ptop->give_ndofe (ndom);
01339 cn = new long [gndofe];
01340 ptop->give_code_numbers (ndom,cn);
01341 globloc (displ,buff,cn,gndofe);
01342 delete [] cn;
01343
01344 }
01345 else{
01346 MPI_Recv(buff,maxncdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01347 }
01348 MPI_Barrier (MPI_COMM_WORLD);
01349
01350 j=nidof;
01351 for (i=0;i<ncdof;i++){
01352 subdispl[j]=buff[i]; j++;
01353 }
01354
01355
01356 mv (krc,buff,ur,nidof,ncdof);
01357 delete [] buff;
01358
01359
01360 buff = new double [maxnbdof];
01361
01362 if (myrank==0){
01363 for (i=1;i<nproc;i++){
01364 if (domproc[i]==0) continue;
01365 extract_from_global_vector (buff,cglhs,i);
01366 MPI_Send(buff,maxnbdof,MPI_DOUBLE,i,ndom,MPI_COMM_WORLD);
01367 }
01368 extract_from_global_vector (buff,cglhs,domproc[nproc]);
01369 }
01370 else{
01371 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01372 }
01373 MPI_Barrier (MPI_COMM_WORLD);
01374
01375
01376 nullvr (vr,nidof);
01377 put_into_local_vector (buff,vr);
01378 delete [] buff;
01379
01380 for (i=0;i<nidof;i++){
01381 rhs[i]=rhs[i]-ur[i]-vr[i];
01382 }
01383
01384 gm->condense (buff,subdispl,rhs,ncdof,3);
01385
01386 delete [] ur;
01387 delete [] vr;
01388 }
01389
01390
01391
01392
01393
01394
01395
01396 void dpfeti::solve_system (gtopology *top,gmatrix *gm,long *domproc,
01397 double *lhs,double *rhs,FILE *out,long mespr)
01398 {
01399 long i,j;
01400 double *condmat,*condvect,*rrhs,*crhs;
01401 time_t t1,t2,t3,t4,t5,t6,t7,t8;
01402
01403 t1 = time (NULL);
01404
01405 if (myrank==0){
01406 condmat = new double [maxncdof*maxncdof];
01407 memset (condmat,0,maxncdof*maxncdof*sizeof(double));
01408 condvect = new double [maxncdof];
01409 memset (condvect,0,maxncdof*sizeof(double));
01410 }
01411 else{
01412 condmat = new double [ncdof*ncdof];
01413 memset (condmat,0,ncdof*ncdof*sizeof(double));
01414 condvect = new double [ncdof];
01415 memset (condvect,0,ncdof*sizeof(double));
01416 }
01417
01418 if (myrank==0){
01419
01420 fprintf (out,"\n\n tncdof %ld",tncdof);
01421
01422 cglhs = new double [tnmdof];
01423 cgrhs = new double [tnmdof];
01424 displ = new double [tncdof];
01425 }
01426
01427
01428 krc = new double [nidof*ncdof];
01429
01430
01431 gm->a12block (krc,ncdof);
01432
01433 crhs = new double [ncdof];
01434 rrhs = new double [nidof];
01435
01436
01437 for (i=0;i<nidof;i++){
01438 rrhs[i]=rhs[i];
01439 }
01440 j=0;
01441 for (i=nidof;i<ndof;i++){
01442 crhs[j]=rhs[i]; j++;
01443 }
01444
01445 t2 = time (NULL);
01446
01447
01448 gm->condense (condmat,lhs,rhs,ncdof,1);
01449
01450 j=0;
01451 for (i=nidof;i<ndof;i++){
01452 condvect[j]=rhs[i];
01453 j++;
01454 }
01455
01456 t3 = time (NULL);
01457
01458 arrmatrix (condmat);
01459
01460 t4 = time (NULL);
01461
01462
01463 vectors_br_bm (top,gm,condmat,condvect,rrhs);
01464
01465
01466 rhs_dpfeti (top,domproc,gm);
01467
01468 t5 = time (NULL);
01469
01470
01471 cg (top,domproc,gm,0,out);
01472
01473 t6 = time (NULL);
01474
01475 if (myrank==0){
01476 fprintf (out,"\n\n L A G R A N G E M U L T I P L I E R S \n");
01477 for (i=0;i<tnmdof;i++){
01478 fprintf (out,"\n %7ld %30.15f",i,cglhs[i]);
01479 }
01480 }
01481
01482
01483 corner_displ (top,domproc,gm);
01484
01485 t7 = time (NULL);
01486
01487 if (myrank==0){
01488 fprintf (out,"\n\n C O R N E R D I S P L A C E M E N T S \n");
01489 fprintf (out,"\n\n tncdof %ld",tncdof);
01490 for (i=0;i<tncdof;i++){
01491 fprintf (out,"\n %7ld %30.15f",i,displ[i]);
01492 }
01493 fprintf (out,"\n\n konec tisku\n");
01494 }
01495
01496
01497 compute_displ (top,gm,domproc,lhs,rrhs);
01498
01499 t8 = time (NULL);
01500
01501 fprintf (out,"\n\n A_rc block assembling %ld",t2-t1);
01502 fprintf (out,"\n\n time of condensation %ld",t3-t2);
01503 fprintf (out,"\n\n A_rr assembling %ld",t4-t3);
01504 fprintf (out,"\n\n b_R and b_M assembling %ld",t5-t4);
01505 fprintf (out,"\n\n solution of reduced system %ld",t6-t5);
01506 fprintf (out,"\n\n corner displacements calcul. %ld",t7-t6);
01507 fprintf (out,"\n\n all displacements calcul. %ld",t8-t7);
01508
01509
01510 delete [] condmat;
01511 delete [] condvect;
01512 delete [] rrhs;
01513 delete [] crhs;
01514
01515 }
01516