00001 #include "psolver.h"
00002 #include <string.h>
00003 #include <math.h>
00004 #include <limits.h>
00005
00006 #include "selnodes.h"
00007
00008
00009
00010 psolver::psolver (int np,int mr, char *nameproc, int nl)
00011 {
00012
00013 nproc=np;
00014
00015 myrank=mr;
00016
00017 strcpy (procName,nameproc);
00018 nameLength=nl;
00019
00020
00021 ndom=0;
00022
00023
00024
00025 tdd = no_par_solver;
00026
00027 trssol = (redsystsolver) 0;
00028
00029 rsmstor = (storagetype) 0;
00030
00031 md = (meshdescription) 0;
00032
00033
00034 zero=1.0e-15;
00035
00036 condfixing = -1;
00037 nmember = -1;
00038
00039 ltg = NULL;
00040
00041 domproc = NULL;
00042
00043
00044 nndom = NULL;
00045
00046 nedom = NULL;
00047
00048 negmdom = NULL;
00049
00050
00051
00052 schcom = NULL;
00053
00054 f1 = NULL;
00055
00056 dpf = NULL;
00057
00058 lpp = NULL;
00059
00060
00061
00062
00063
00064 selnodschur = NULL;
00065 selnodfeti = NULL;
00066
00067
00068 ptop=NULL;
00069
00070
00071 fixnodes = NULL;
00072 methodcondcor = nomethod;
00073 typecondcur = notype;
00074 nmembercur = -1;
00075 nuserdefnod = 0;
00076 userdefnod = NULL;
00077 typecondsurf = notype;
00078 nmembersurf = -1;
00079 nring = -1;
00080 ring = NULL;
00081
00082
00083 ptopjk = NULL;
00084
00085 ssle = new slesolv();
00086 }
00087
00088
00089
00090 psolver::~psolver()
00091 {
00092 delete [] ltg;
00093 delete [] domproc;
00094 delete [] nndom;
00095 delete [] nedom;
00096 delete [] negmdom;
00097
00098 delete schcom;
00099 delete f1;
00100 delete dpf;
00101 delete lpp;
00102
00103
00104
00105 delete selnodschur;
00106 delete selnodfeti;
00107
00108 delete ptop;
00109 delete ssle;
00110 delete ptopjk;
00111 }
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 void psolver::initiate (gtopology *top,int argc,const char**argv)
00123 {
00124
00125 ptop = new partop (nproc,myrank,ndom,md,procName,nameLength);
00126 ptopjk = new partopjk (nproc,myrank,md,procName,nameLength);
00127
00128 switch (tdd){
00129
00130
00131
00132 case no_par_solver:{
00133 break;
00134 }
00135
00136 case schurcompldd:{
00137 schcom = new schurcompl (nproc,myrank,ndom);
00138
00139 schcom->trssol = trssol;
00140 schcom->rsmstor = rsmstor;
00141 schcom->ssle->initiate (ssle);
00142 schcom->dmstor = dmstor;
00143 schcom->nicgsch = nicg;
00144 schcom->errcgsch = errcg;
00145 schcom->zero = zero;
00146
00147 break;
00148 }
00149 case fetidd:{
00150 f1 = new feti1 (nproc,myrank,ndom);
00151
00152 f1->enrbm = enrbm;
00153 f1->lithr = lithr;
00154 f1->nicg = nicg;
00155 f1->errcg = errcg;
00156 f1->zero = zero;
00157 f1->fetiprecond = fetiprecond;
00158 break;
00159 }
00160 case dpfetidd:{
00161 dpf = new dpfeti (nproc,myrank,ndom);
00162
00163 dpf->rsmstor = rsmstor;
00164 dpf->ssle->initiate (ssle);
00165 dpf->nicgdpfeti = nicg;
00166 dpf->errcgdpfeti = errcg;
00167 dpf->zero = zero;
00168
00169 if(fixnodes != NULL){
00170 delete fixnodes;
00171 }
00172 fixnodes = new fixnodesel (nproc,myrank,ndom,md,domproc,procName,nameLength,mespr);
00173 fixnodes->condfixing = condfixing;
00174 fixnodes->methodcondcor = methodcondcor;
00175 fixnodes->typecondcur = typecondcur;
00176 fixnodes->nmembercur = nmembercur;
00177 fixnodes->nuserdefnod = nuserdefnod;
00178 fixnodes->userdefnod = userdefnod;
00179 fixnodes->typecondsurf = typecondsurf;
00180 fixnodes->nmembersurf = nmembersurf;
00181 fixnodes->nring = nring;
00182 fixnodes->ring = ring;
00183 break;
00184 }
00185 case parconjuggrad:{
00186
00187
00188
00189
00190 break;
00191 }
00192 case boundparconjuggrad:{
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210 break;
00211 }
00212 case layered_plate:{
00213 lpp = new lplate (nproc,myrank,ndom,top);
00214
00215 lpp->nicg = nicg;
00216 lpp->errcg = errcg;
00217 break;
00218 }
00219 default:{
00220 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
00221 }
00222 }
00223
00224
00225 }
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237 void psolver::read (XFILE *in,gtopology *gt,long nd,storagetype tsm,long mes)
00238 {
00239
00240 mespr=mes;
00241
00242
00243 ndom=nd;
00244
00245
00246 dmstor = tsm;
00247
00248
00249
00250
00251
00252
00253
00254
00255 xfscanf (in,"%k%m","mesh_description",&meshdescription_kwdset,(int*)&md);
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267 xfscanf (in,"%k%m","type_of_domain_decomposition_solver",&domdectype_kwdset,(int*)&tdd);
00268
00269
00270 switch (tdd){
00271
00272
00273
00274 case no_par_solver:{
00275
00276
00277
00278
00279 break;
00280 }
00281
00282
00283
00284
00285 case schurcompldd:{
00286 if (mespr==1) fprintf (stdout,"\n primal domain decomposition (Schur complement method) is used");
00287
00288
00289
00290
00291
00292 xfscanf (in,"%k%m","type_of_reduced_solver",&redsystsolver_kwdset,(int*)&trssol);
00293 switch (trssol){
00294 case master_sol:{
00295 if (mespr==1) fprintf (stdout,"\n reduced system of equations is solved sequentially on master processor");
00296 xfscanf (in,"%k%m","storage_of_reduced_system_matrix",&storagetype_kwdset,(int*)&rsmstor);
00297
00298 ssle->read (gt,in,mespr);
00299
00300 break;
00301 }
00302 case paral_sol:{
00303 if (mespr==1) fprintf (stdout,"\n reduced system of equations is solved in parallel");
00304 xfscanf (in,"%k%ld","number_of_iterations",&nicg);
00305 xfscanf (in,"%k%le","error_of_computation",&errcg);
00306 break;
00307 }
00308 default:{
00309 par_print_err(myrank+1, procName, "unknown reduced system solver is required",__FILE__, __LINE__, __func__);
00310 }
00311 }
00312 break;
00313 }
00314
00315
00316
00317 case fetidd:{
00318 if (mespr==1) fprintf (stdout,"\n dual domain decomposition (Finite Element Tearing and Interconnecting method) is used");
00319 xfscanf (in,"%k%ld","est_number_of_RBM",&enrbm);
00320 xfscanf (in,"%k%lf","treshold_lin_dep",&lithr);
00321
00322 xfscanf (in,"%k%ld","number_of_iterations",&nicg);
00323 xfscanf (in,"%k%le","error_of_computation",&errcg);
00324
00325
00326
00327
00328
00329 xfscanf (in,"%k%m","type_of_preconditioning",&fetiprecond_kwdset,&fetiprecond);
00330 switch (fetiprecond){
00331 case nofetiprecond:{
00332 fprintf (stdout,"\n No preconditioning is used for FETI method");
00333 break;}
00334 case lumped:{
00335 fprintf (stdout,"\n Lumped preconditioning is used for FETI method");
00336 }
00337 case dirichlet:{
00338 fprintf (stdout,"\n Dirichlet preconditioning is used for FETI method");
00339 break;
00340 }
00341 }
00342 break;
00343 }
00344
00345
00346
00347 case dpfetidd:{
00348 long i;
00349 if (mespr==1) fprintf (stdout,"\n dual-primal finite element tearing and interconnecting method is used");
00350
00351
00352 xfscanf (in,"%k%m","type_of_reduced_solver",&redsystsolver_kwdset,(int*)&trssol);
00353 switch (trssol){
00354 case master_sol:{
00355 if (mespr==1) fprintf (stdout,"\n reduced system of equations is solved sequentially on master processor");
00356 xfscanf (in,"%k%m","storage_of_reduced_system_matrix",&storagetype_kwdset,(int*)&rsmstor);
00357 ssle->read (gt,in,mespr);
00358 xfscanf (in,"%k%ld","number_of_iterations",&nicg);
00359 xfscanf (in,"%k%le","error_of_computation",&errcg);
00360 break;
00361 }
00362 default:{
00363 par_print_err(myrank+1, procName, "unknown reduced system solver is required",__FILE__, __LINE__, __func__);
00364 break;
00365 }
00366 }
00367 xfscanf (in,"%k%m","condensation_of_fixing_nodes",&condfixing_kwdset,&condfixing);
00368
00369 switch(condfixing){
00370 case nocondconer:{
00371 if (mespr==1) fprintf (stdout,"\n Fixing nodes will not be condensed");
00372 methodcondcor = nomethod;
00373 break;
00374 }
00375 case automatic:{
00376 if (mespr==1) fprintf (stdout,"\n Fixing nodes will be condensed automatically");
00377 break;
00378 }
00379 case userdef:{
00380 if (mespr==1) fprintf (stdout,"\n Fixing nodes will be condensed by user definition");
00381 xfscanf (in,"%k%m","place_of_condensation",&methodcond_kwdset,&methodcondcor);
00382 if (methodcondcor == cursurf || methodcondcor == curvecond){
00383 if(methodcondcor == curvecond){
00384 if (mespr==1) fprintf (stdout,"\n Additional fixing nodes will be added on boundary curve only");
00385
00386 }
00387 else{
00388 if (mespr==1) fprintf (stdout,"\n Additional fixing nodes will be added on boundary curve and boundary surfaces");
00389 }
00390 xfscanf (in,"%k%m","method_of_condensation_curve",&typecondfixing_kwdset,&typecondcur);
00391 switch(typecondcur){
00392 case centroid:{
00393 if (mespr==1) fprintf (stdout,"\n New fixing nodes will be added into centroid between two fixing nodes");
00394 break;
00395 }
00396 case nth_memb:{
00397 xfscanf (in,"%ld",&nmembercur);
00398 if (mespr==1) fprintf (stdout,"\n New fixing nodes will be added into %ld nodes between two fixing nodes",nmembercur);
00399 break;
00400 }
00401 case all_memb:{
00402 if (mespr==1) fprintf (stdout,"\n New fixing nodes will be added into all nodes between two fixing nodes");
00403 break;
00404 }
00405 case rand_memb:{
00406 xfscanf (in,"%ld",&nmembercur);
00407 if (mespr==1) fprintf (stdout,"\n New fixing nodes will be added into %ld nodes between two fixing nodes on random positions",nmembercur);
00408 break;
00409 }
00410 case n_part_curve:{
00411 xfscanf (in,"%ld",&nmembercur);
00412 if (mespr==1) fprintf (stdout,"\n Boundary curve between two fixing nodes will be cut into %ld parts",nmembercur);
00413 break;
00414 }
00415 case userposdef:{
00416 xfscanf (in,"%ld",&nuserdefnod);
00417 if (mespr==1) fprintf (stdout,"\n New fixing nodes will be added into %ld nodes which user specified in input\n",nuserdefnod);
00418 userdefnod = new long[nuserdefnod];
00419 for(i = 0; i < nuserdefnod; i++){
00420 xfscanf (in,"%ld",&userdefnod[i]);
00421 }
00422 break;
00423 }
00424 }
00425 }
00426 if (methodcondcor == cursurf || methodcondcor == surfacecond){
00427 xfscanf (in,"%k%m","method_of_condensation_surface",&typecondfixing_kwdset,&typecondsurf);
00428 switch(typecondsurf){
00429 case centroid:{
00430 if (mespr==1) fprintf (stdout,"\n New fixing nodes will be added into centroid of boundary surface");
00431 break;
00432 }
00433 case rand_memb:{
00434 xfscanf (in,"%ld",&nmembersurf);
00435 if (mespr==1) fprintf (stdout,"\n %ld nodes will be added on boundary surface",nmembersurf);
00436 break;
00437 }
00438 case all_memb:{
00439 if (mespr==1) fprintf (stdout,"\n New fixing nodes will be added into all nodes between two fixing nodes");
00440 break;
00441 }
00442 case n_mark:{
00443 xfscanf (in,"%ld",&nmembersurf);
00444 break;
00445 }
00446 case chose_ring:{
00447 xfscanf (in,"%ld",&nring);
00448 ring = new long[nring];
00449 for(i = 0; i < nring; i++){
00450 xfscanf (in,"%ld",&ring[i]);
00451 }
00452 break;
00453 }
00454 case userposdef:{
00455 xfscanf (in,"%ld",&nuserdefnod);
00456 if (mespr==1) fprintf (stdout,"\n New fixing nodes will be added into %ld nodes which user specified in input\n",nuserdefnod);
00457 userdefnod = new long[nuserdefnod];
00458 for(i = 0; i < nuserdefnod; i++){
00459 xfscanf (in,"%ld",&userdefnod[i]);
00460 }
00461 break;
00462 }
00463 }
00464 }
00465 break;
00466 }
00467 }
00468 break;
00469 }
00470
00471
00472
00473
00474 case parconjuggrad:{
00475 if (mespr==1) fprintf (stdout,"\n parallel conjugate gradient method is used");
00476
00477 xfscanf (in,"%k%ld","number_of_iterations",&nicg);
00478 xfscanf (in,"%k%le","error_of_computation",&errcg);
00479
00480 break;
00481 }
00482
00483
00484
00485 case boundparconjuggrad:{
00486 if (mespr==1) fprintf (stdout,"\n boundary version of parallel conjugate gradient method is used");
00487
00488 xfscanf (in,"%k%ld","number_of_iterations",&nicg);
00489 xfscanf (in,"%k%le","error_of_computation",&errcg);
00490
00491
00492
00493
00494
00495
00496 xfscanf (in,"%k%m","type_of_preconditioning",&parcgprec_kwdset,&prec);
00497 switch(prec){
00498 case noprec:{
00499 if (mespr==1) fprintf (stdout,"\n boundary version of parallel conjugate gradient method will not be preconditioned");
00500 break;
00501 }
00502 case petscilu:{
00503 if (mespr==1) fprintf (stdout,"\n boundary version of parallel conjugate gradient method will be preconditioned by PETSC ILU");
00504 break;
00505 }
00506 case pardiagprec:{
00507 if (mespr==1) fprintf (stdout,"\n boundary version of parallel conjugate gradient method will be preconditioned by Jacobi preconditioning");
00508 }
00509 }
00510 break;
00511 }
00512
00513
00514
00515 case layered_plate:{
00516
00517 if (mespr==1) fprintf (stdout,"\n layered linear plate problem is solved with help of orthonormalization of constraints");
00518
00519
00520 ssle->read (gt,in,mespr);
00521
00522
00523 xfscanf (in,"%k%ld","number_of_iterations",&nicg);
00524 xfscanf (in,"%k%le","error_of_computation",&errcg);
00525
00526 break;
00527 }
00528
00529
00530
00531 default:{
00532 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
00533 }
00534 }
00535
00536 }
00537
00538
00539
00540
00541
00542
00543
00544
00545 void psolver::print (FILE *out)
00546 {
00547
00548 fprintf (out,"%d\n",md);
00549
00550 fprintf (out,"%d\n", tdd);
00551
00552 switch (tdd){
00553 case no_par_solver:{
00554
00555
00556
00557
00558 break;
00559 }
00560
00561
00562
00563 case schurcompldd:{
00564
00565
00566 fprintf (out,"%d ",(int)trssol);
00567
00568 switch (trssol){
00569 case master_sol:{
00570 fprintf (out,"%d\n",(int)rsmstor);
00571
00572 ssle->print (out);
00573 fprintf (out,"\n");
00574
00575 break;
00576 }
00577 case paral_sol:{
00578 fprintf (out,"%ld %e\n",nicg,errcg);
00579 break;
00580 }
00581 default:{
00582 par_print_err(myrank+1, procName, "unknown reduced system solver is required",__FILE__, __LINE__, __func__);
00583 }
00584 }
00585
00586 break;
00587 }
00588
00589
00590
00591 case fetidd:{
00592 fprintf (out,"%ld %e ",enrbm,lithr);
00593 fprintf (out,"%ld %e\n",nicg,errcg);
00594 fprintf (out,"%ld\n",fetiprecond);
00595 switch (fetiprecond){
00596 case nofetiprecond:{
00597 fprintf (stdout,"\n No preconditioning is used for FETI method");
00598 break;}
00599 case lumped:{
00600 fprintf (stdout,"\n Lumped preconditioning is used for FETI method");
00601 }
00602 case dirichlet:{
00603 fprintf (stdout,"\n Dirichlet preconditioning is used for FETI method");
00604 break;
00605 }
00606 }
00607 break;
00608 }
00609
00610
00611
00612 case dpfetidd:{
00613
00614 fprintf (out,"%d\n",(int)trssol);
00615 switch (trssol){
00616 case master_sol:{
00617 fprintf (out,"%d\n",(int)rsmstor);
00618 ssle->print (out);
00619 fprintf (out,"%ld %e\n",nicg,errcg);
00620 break;
00621 }
00622 default:{
00623 par_print_err(myrank+1, procName, "unknown reduced system solver is required",__FILE__, __LINE__, __func__);
00624 }
00625 }
00626 long i;
00627 fprintf (out,"%ld\n",condfixing);
00628 if(condfixing == userdef){
00629 fprintf (out,"%ld\n",methodcondcor);
00630 if (methodcondcor == cursurf || methodcondcor == curvecond){
00631 fprintf (out,"%ld\n",typecondcur);
00632 if(typecondcur == nth_memb || typecondcur == rand_memb || typecondcur == n_part_curve){
00633 fprintf (out,"%ld\n",nmembercur);
00634 }
00635 if(typecondcur == userposdef){
00636 fprintf (out,"%ld\n",nuserdefnod);
00637 for(i = 0; i < nuserdefnod; i++){
00638 fprintf (out,"%ld ",userdefnod[i]);
00639 }
00640 fprintf (out,"\n");
00641 }
00642 break;
00643 }
00644 if (methodcondcor == cursurf || methodcondcor == surfacecond){
00645 fprintf (out,"%ld",typecondsurf);
00646 switch(typecondsurf){
00647 case rand_memb:{
00648 fprintf (out," %ld\n",nmembersurf);
00649 break;
00650 }
00651 case n_mark:{
00652 fprintf (out," %ld",nmembersurf);
00653 break;
00654 }
00655 case chose_ring:{
00656 fprintf (out," %ld ",nring);
00657 for(i = 0; i < nring; i++){
00658 fprintf (out," %ld",ring[i]);
00659 }
00660 fprintf (out,"\n");
00661 break;
00662 }
00663 case userposdef:{
00664 fprintf (out," %ld\n",nuserdefnod);
00665 for(i = 0; i < nuserdefnod; i++){
00666 fprintf (out," %ld ",userdefnod[i]);
00667 }
00668 fprintf (out,"\n");
00669 break;
00670 }
00671 }
00672 }
00673 }
00674
00675
00676 break;
00677 }
00678
00679
00680
00681 case parconjuggrad:{
00682
00683 fprintf (out,"%ld %le",nicg,errcg);
00684
00685 break;
00686 }
00687
00688
00689
00690 case boundparconjuggrad:{
00691
00692 fprintf (out,"%ld %le",nicg,errcg);
00693 fprintf (out," %ld",prec);
00694
00695 break;
00696 }
00697
00698 case layered_plate:{
00699
00700 ssle->print (out);
00701 fprintf (out,"\n");
00702
00703
00704 fprintf (out,"%ld %e\n",nicg,errcg);
00705
00706 break;
00707 }
00708
00709 default:{
00710 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
00711 }
00712 }
00713
00714 fprintf (out,"\n");
00715
00716 }
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728 void psolver::read_ltg (XFILE *in,gtopology *top,FILE *out)
00729 {
00730 long i;
00731
00732 if (tdd != no_par_solver){
00733
00734 ltg = new long [top->nn];
00735
00736 for (i=0;i<top->nn;i++){
00737 xfscanf (in,"%ld",ltg+i);
00738
00739 if (ltg[i]<0)
00740 ltg[i]=0-ltg[i];
00741
00742 ltg[i]--;
00743 }
00744 }
00745
00746 }
00747
00748
00749
00750
00751
00752
00753
00754
00755 void psolver::procdomcorr ()
00756 {
00757 long i,j;
00758 MPI_Status stat;
00759
00760 if (myrank==0){
00761 domproc = new long [nproc];
00762 domproc[0]=ndom;
00763 for (i=1;i<nproc;i++){
00764 MPI_Recv (&j,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
00765 domproc[stat.MPI_TAG]=j;
00766 }
00767 }
00768 else{
00769 MPI_Send (&ndom,1,MPI_LONG,0,myrank,MPI_COMM_WORLD);
00770 }
00771 MPI_Barrier (MPI_COMM_WORLD);
00772
00773 }
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796 long psolver::ordering (gtopology *top,FILE *out,char *proc_name)
00797 {
00798 long i,ndof;
00799 long *schurnodes,*fetinodes;
00800
00801 top->ordering = new long [top->nn];
00802
00803
00804 switch (top->nodren){
00805 case no_renumbering:{
00806 for (i=0;i<top->nn;i++){
00807 top->ordering[i]=i;
00808 }
00809 break;
00810 }
00811 case cuthill_mckee:
00812 case rev_cuthill_mckee:{
00813 top->cuthill_mckee_renumb (out);
00814 break;
00815 }
00816 case sloan:{
00817 top->sloan_renumb(out);
00818 break;
00819 }
00820 default:{
00821 print_err("unknown type of renumbering is required",__FILE__,__LINE__,__func__);
00822 }
00823 }
00824
00825
00826 switch (tdd){
00827
00828
00829
00830 case no_par_solver:{
00831 ndof = top->codenum_generation (out);
00832 break;
00833 }
00834
00835
00836
00837
00838 case schurcompldd:{
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883 ndof = ptopjk->vse (ltg,top,domproc,out,proc_name);
00884 schcom->initiate (ptopjk,out);
00885
00886 break;
00887 }
00888
00889
00890
00891
00892 case fetidd:{
00893
00894 ptop->initiation (top,ltg);
00895 ptop->numbers_of_all_nodes_on_subdomains (domproc,out);
00896 ptop->compute_multiplicity (ltg,domproc,out);
00897 ptop->find_boundary_nodes (ltg,domproc,out);
00898 ptop->rewrite_ltg (ltg);
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908 f1->subdomain_ordering (top,ltg,out);
00909 fetinodes = new long [top->nn];
00910
00911 fprintf (out,"\n\n kontrola fetinodes");
00912 for (i=0;i<top->nn;i++){
00913 if (ltg[i]>-1)
00914 fetinodes[i]=ltg[i];
00915 else
00916 fetinodes[i]=-1;
00917
00918 fprintf (out,"\n fetinodes %ld %ld ltg %ld",i,fetinodes[i],ltg[i]);
00919 }
00920
00921
00922 selnodfeti = new selnodes (nproc,myrank,ndom,top->nn,fetinodes,md,out,mespr);
00923 selnodfeti->number_of_selected_nodes (domproc,out);
00924 selnodfeti->nodes_on_master (domproc,out);
00925 selnodfeti->node_multiplicity (out);
00926 selnodfeti->number_all_dofs (top,domproc,out);
00927 selnodfeti->ndofn_on_master (top,domproc,out);
00928 selnodfeti->dof_indicators (top,domproc,out);
00929 selnodfeti->group_local_nodes (out);
00930 selnodfeti->dof_feti (out);
00931 selnodfeti->number_contrib (domproc,out);
00932 selnodfeti->contrib_dofs (top,domproc,out);
00933
00934
00935 selnodfeti->dof_multiplicity (out);
00936
00937 f1->initiate (selnodfeti,out);
00938
00939 ndof=f1->ndof;
00940
00941 delete [] fetinodes;
00942 break;
00943 }
00944
00945
00946
00947 case dpfetidd:{
00948
00949 ptop->prepare_data (domproc,ltg,top,out,proc_name);
00950 bool selection;
00951
00952 selection=fixnodes->check_ltg(ltg,ptop->nn,out);
00953 switch(selection){
00954 case true:{
00955
00956
00957 fixnodes->initiate(top,ptop,out);
00958
00959 fixnodes->fixing_detection(ltg);
00960
00961 md = bound_nodes;
00962 delete fixnodes;
00963 break;
00964 }
00965 case false:{
00966 delete fixnodes;
00967 break;
00968 }
00969 }
00970
00971
00972
00973 for (i=0;i<top->nn;i++){
00974 if (ltg[i]>-2)
00975 top->gnodes[i].ai=-1;
00976 else
00977 top->gnodes[i].ai = 0-ltg[i]-2;
00978 }
00979
00980 fprintf (out,"\n\n kontrola pole ltg");
00981 for (i=0;i<top->nn;i++){
00982 fprintf (out,"\n ltg %ld %ld",i+1,ltg[i]);
00983 }
00984 fprintf (out,"\n");
00985
00986 dpf->dpfetiordering (top,ltg);
00987 ndof=dpf->ndof;
00988
00989 fetinodes = new long [top->nn];
00990
00991
00992 for (i=0;i<top->nn;i++){
00993 if (ltg[i]>-1)
00994 fetinodes[i]=ltg[i];
00995 else
00996 fetinodes[i]=-1;
00997
00998
00999 }
01000
01001 selnodfeti = new selnodes (nproc,myrank,ndom,top->nn,fetinodes,md,out,mespr);
01002 selnodfeti->number_of_selected_nodes (domproc,out);
01003 selnodfeti->nodes_on_master (domproc,out);
01004 selnodfeti->node_multiplicity (out);
01005 selnodfeti->number_all_dofs (top,domproc,out);
01006 selnodfeti->ndofn_on_master (top,domproc,out);
01007 selnodfeti->dof_indicators (top,domproc,out);
01008 selnodfeti->group_local_nodes (out);
01009 selnodfeti->dof_feti (out);
01010 selnodfeti->number_contrib (domproc,out);
01011 selnodfeti->contrib_dofs (top,domproc,out);
01012
01013
01014 selnodfeti->dof_multiplicity (out);
01015
01016 fprintf (out,"\n\n\n KONEC FETI, ZACATEK SCHURA \n\n\n");
01017 fprintf (out,"\n\n kontrola schurnodes");
01018
01019 schurnodes = new long [top->nn];
01020 for (i=0;i<top->nn;i++){
01021 if (ltg[i]<-1)
01022 schurnodes[i]=0-ltg[i]-2;
01023 else
01024 schurnodes[i]=-1;
01025
01026 }
01027 selnodschur = new selnodes (nproc,myrank,ndom,top->nn,schurnodes,md,out,mespr);
01028 selnodschur->number_of_selected_nodes (domproc,out);
01029 selnodschur->nodes_on_master (domproc,out);
01030 selnodschur->node_multiplicity (out);
01031 selnodschur->number_all_dofs (top,domproc,out);
01032
01033 selnodschur->ndofn_on_master (top,domproc,out);
01034 selnodschur->dof_indicators (top,domproc,out);
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046 delete [] fetinodes;
01047 break;
01048 }
01049 case parconjuggrad:{
01050
01051
01052
01053
01054
01055
01056 break;
01057 }
01058
01059
01060 case boundparconjuggrad:{
01061 ptop->prepare_data (domproc,ltg,top,out,proc_name);
01062
01063
01064 schurnodes = new long [top->nn];
01065
01066
01067 for (i=0;i<top->nn;i++){
01068 if (ptop->nodmultip[i]>1)
01069 schurnodes[i]=1;
01070 else
01071 schurnodes[i]=-1;
01072 }
01073
01074 selnodschur = new selnodes (nproc,myrank,md,domproc,ptop->nnsd,schurnodes,
01075 ptop->tnbn,ptop->icmultip,ptop->nodmultip,ptop->icnbn,ptop->gnbndom,
01076 ptop->gnbncn,ptop->sid,
01077 out,proc_name,mespr);
01078
01079 selnodschur->node_coarse_numbers (out);
01080 selnodschur->number_all_dofs (top,domproc,out);
01081
01082 selnodschur->ndofn_on_master (top,domproc,out);
01083 selnodschur->dof_indicators (top,domproc,out);
01084
01085
01086 if (md==all_nodes){
01087 selnodschur->schur_ordering (top,ptop->dofind,out);
01088 }
01089 if (md==bound_nodes){
01090 selnodschur->schur_ordering (top,out);
01091 }
01092
01093 ndof = ptop->schur_ordering (top,domproc,out,proc_name);
01094
01095
01096
01097 delete [] schurnodes;
01098 break;
01099 }
01100 case layered_plate:{
01101 lpp->ndof = top->gencodnum ();
01102 lpp->globcnnum_lpp (top,domproc,out);
01103 ndof=lpp->ndof;
01104 break;
01105 }
01106 default:{
01107 par_print_err(myrank,"unknown type of domain decomposition solver is required",__FILE__,__LINE__,__func__);
01108 }
01109 }
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148 return ndof;
01149 }
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160 void psolver::constr_mat (double *th,FILE *out)
01161 {
01162 switch (tdd){
01163 case layered_plate:{
01164 lpp->constraintmat (th,domproc,out);
01165 break;
01166 }
01167 default:{
01168 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
01169 }
01170 }
01171
01172 }
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186 void psolver::par_linear_solver (gtopology *top,gmatrix *gm,
01187 double *lhs,double *rhs,FILE *out,long mespr)
01188 {
01189 switch (tdd){
01190 case schurcompldd:{
01191 schcom->solve_system (top,gm,domproc,lhs,rhs,out);
01192 break;
01193 }
01194 case fetidd:{
01195 f1->solve_system (top,gm,domproc,lhs,rhs,out);
01196 break;
01197 }
01198 case dpfetidd:{
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211 break;
01212 }
01213 case parconjuggrad:{
01214
01215 break;
01216 }
01217 case boundparconjuggrad:{
01218
01219 break;
01220 }
01221 case layered_plate:{
01222 lpp->solve_system (top,gm,domproc,lhs,rhs,out);
01223 break;
01224 }
01225 default:{
01226 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
01227 }
01228 }
01229
01230 }
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242 void psolver::gather_bound_vect (double *lv,double *gv)
01243 {
01244 switch (tdd){
01245 case schurcompldd:{
01246 schcom->gather_bound_vect (lv,gv,domproc);
01247 break;
01248 }
01249 default:{
01250 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
01251 }
01252 }
01253
01254 }
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266 void psolver::scatter_bound_vect (double *lv,double *gv)
01267 {
01268 switch (tdd){
01269 case schurcompldd:{
01270 schcom->scatter_bound_vect (lv,gv,domproc);
01271 break;
01272 }
01273 default:{
01274 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
01275 }
01276 }
01277
01278 }
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290 double psolver::unbalanced_values (double *lv,FILE *out)
01291 {
01292 double norm;
01293
01294 switch (tdd){
01295 case schurcompldd:{
01296 norm = schcom->unbalanced_values (lv,domproc,out);
01297 break;
01298 }
01299 default:{
01300 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
01301 }
01302 }
01303
01304 return norm;
01305 }
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315 double psolver::pss (double *lv1,double *lv2,FILE *out)
01316 {
01317 long i,n;
01318 double dp=0.0;
01319 double *gv1,*gv2;
01320 MPI_Status stat;
01321
01322 switch (tdd){
01323 case schurcompldd:{
01324 if (myrank==0){
01325
01326 n=schcom->ndofcp;
01327
01328 gv1 = new double [n];
01329 gv2 = new double [n];
01330 nullv (gv1,n);
01331 nullv (gv2,n);
01332
01333
01334 dp=schcom->pss_gather_bound_vect (lv1,gv1,lv2,gv2,domproc);
01335
01336 for (i=0;i<n;i++){
01337 dp+=gv1[i]*gv2[i];
01338 }
01339 delete [] gv1;
01340 delete [] gv2;
01341 }
01342 else{
01343
01344
01345 dp=schcom->pss_gather_bound_vect (lv1,gv1,lv2,gv2,domproc);
01346 }
01347
01348 if (myrank==0){
01349 for (i=1;i<nproc;i++){
01350 MPI_Send(&dp,1,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD);
01351 }
01352 }
01353 else{
01354 MPI_Recv (&dp,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01355 }
01356 MPI_Barrier (MPI_COMM_WORLD);
01357 break;
01358 }
01359 default:{
01360 par_print_err(myrank,"unknown type of domain decomposition is required",__FILE__,__LINE__,__func__);
01361 }
01362 }
01363
01364 return dp;
01365 }
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383 long psolver::selected_norm_calculation (long cid,double err,double thresh,gtopology *top,double *lv1,double *gv1,double *rhs,double *grhs,long n,FILE *out)
01384 {
01385 long stop;
01386
01387 switch (tdd){
01388 case schurcompldd:{
01389 stop = schcom->selected_norm_calculation (cid,err,thresh,ptopjk,top,lv1,gv1,rhs,grhs,domproc,n,out);
01390 break;
01391 }
01392 default:{
01393 par_print_err(myrank+1, procName, "unknown domain decomposition solver is required",__FILE__, __LINE__, __func__);
01394 }
01395 }
01396
01397 return stop;
01398 }
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424 void psolver::computation_stat (gtopology *top,gmatrix *gm)
01425 {
01426 long i,j,buffsize;
01427 long *buff;
01428 MPI_Status stat;
01429
01430 buffsize = 3;
01431 buff = new long [3];
01432
01433
01434 buff[0]=top->nn;
01435
01436 buff[1]=top->ne;
01437
01438 buff[2]=gm->give_negm ();
01439
01440 if (myrank==0){
01441 nndom = new long [nproc];
01442 nedom = new long [nproc];
01443 negmdom = new long [nproc];
01444
01445 j=domproc[0];
01446 nndom[j]=buff[0];
01447 nedom[j]=buff[1];
01448 negmdom[j]=buff[2];
01449
01450 for (i=1;i<nproc;i++){
01451 MPI_Recv (buff,buffsize,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01452 j=domproc[stat.MPI_TAG];
01453 nndom[j]=buff[0];
01454 nedom[j]=buff[1];
01455 negmdom[j]=buff[2];
01456 }
01457
01458 minnn=LONG_MAX;
01459 maxnn=0;
01460 for (i=0;i<nproc;i++){
01461 if (minnn>nndom[i])
01462 minnn=nndom[i];
01463 if (maxnn<nndom[i])
01464 maxnn=nndom[i];
01465 }
01466
01467 minne=LONG_MAX;
01468 maxne=0;
01469 for (i=0;i<nproc;i++){
01470 if (minne>nedom[i])
01471 minne=nedom[i];
01472 if (maxne<nedom[i])
01473 maxne=nedom[i];
01474 }
01475
01476 minnegm=LONG_MAX;
01477 maxnegm=0;
01478 for (i=0;i<nproc;i++){
01479 if (minnegm>negmdom[i])
01480 minnegm=negmdom[i];
01481 if (maxnegm<negmdom[i])
01482 maxnegm=negmdom[i];
01483 }
01484
01485 }
01486 else{
01487 MPI_Send (buff,buffsize,MPI_LONG,0,myrank,MPI_COMM_WORLD);
01488 }
01489 MPI_Barrier (MPI_COMM_WORLD);
01490
01491 delete [] buff;
01492 }
01493
01494
01495 void psolver::computation_stat_print (FILE *out)
01496 {
01497 long i;
01498
01499 if (myrank==0){
01500 fprintf (out,"\n\n\n informations about computation \n");
01501
01502 fprintf (out,"\n i NN NE NEGM");
01503 for (i=0;i<nproc;i++){
01504 fprintf (out,"\n %3ld %9ld %9ld %12ld",i+1,nndom[i],nedom[i],negmdom[i]);
01505 }
01506
01507 fprintf (out,"\n");
01508
01509 fprintf (out,"\n minimum number of nodes %ld",minnn);
01510 fprintf (out,"\n maximum number of nodes %ld",maxnn);
01511
01512 fprintf (out,"\n minimum number of elements %ld",minne);
01513 fprintf (out,"\n maximum number of elements %ld",maxne);
01514
01515 fprintf (out,"\n minimum number of matrix entries %ld",minnegm);
01516 fprintf (out,"\n maximum number of matrix entries %ld",maxnegm);
01517 }
01518 }
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535 long psolver::compare_on_master (double val, double lim)
01536 {
01537 long ret = -2;
01538 long i;
01539 MPI_Status stat;
01540
01541 if (myrank==0)
01542 {
01543 if (val > lim)
01544 ret = 1;
01545 else
01546 {
01547 if (val < lim)
01548 ret = -1;
01549 else
01550 ret = 0;
01551 }
01552 for (i=1;i<nproc;i++)
01553 MPI_Send(&ret,1,MPI_LONG,i,myrank,MPI_COMM_WORLD);
01554 }
01555 else
01556 MPI_Recv(&ret,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
01557
01558 MPI_Barrier (MPI_COMM_WORLD);
01559
01560 return ret;
01561 }