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