00001 #include "seqschur.h" 00002 #include "slesolv.h" 00003 #include "gmatrix.h" 00004 #include "skyline.h" 00005 #include <limits.h> 00006 #include <string.h> 00007 #include <math.h> 00008 #include <time.h> 00009 00010 /** 00011 constructor 00012 00013 JK, 13.5.2009 00014 */ 00015 seqschur::seqschur () 00016 { 00017 // number of subdomains 00018 ns=0; 00019 00020 // type of solver of reduced/coarse problem 00021 trssol = (redsystsolver) 0; 00022 // type of storage of matrix of reduced/coarse problem 00023 rsmstor = (storagetype) 0; 00024 // maximum number of iterations in the conjugate gradient method used for solution of reduced problem 00025 nicg=0; 00026 // required/prescribed norm of residuum 00027 errcg=0.0; 00028 // actual number of iterations (number of iterations performed during the solution) 00029 anicg=0; 00030 // attained norm of residuum 00031 aerrcg=0.0; 00032 00033 // computer zero 00034 zero=0.0; 00035 00036 00037 00038 00039 00040 00041 00042 // number of DOFs on subdomain 00043 ndof=0; 00044 // number of internal DOFs on subdomain 00045 nidof=0; 00046 // number of boundary/interface DOFs on subdomain 00047 nbdof=0; 00048 // maximum number of boundary/interface DOFs on subdomain 00049 maxnbdof=0; 00050 00051 00052 00053 00054 // array containing numbers of boundary/interface DOFs on subdomains 00055 nbdofmas=NULL; 00056 00057 00058 gtop=NULL; 00059 arr=NULL; 00060 00061 // solver of system of linear eqautions 00062 ssle = NULL; 00063 00064 } 00065 00066 seqschur::~seqschur() 00067 { 00068 delete ssle; 00069 } 00070 00071 00072 /** 00073 function reads basic data 00074 00075 @param i - type of solver of system of equations 00076 @param in - input file 00077 00078 JK, 13.5.2009 00079 */ 00080 void seqschur::read (gtopology *top,long mespr,XFILE *in) 00081 { 00082 // number of subdomains 00083 xfscanf (in,"%ld",&ns); 00084 00085 // type of reduced solver 00086 // trssol = master_sol = 1 - solution by a direct method 00087 // trssol = paral_sol = 2 - solution by an iterative method 00088 xfscanf (in,"%d",(int*)&trssol); 00089 switch (trssol){ 00090 case master_sol:{ 00091 if (mespr==1) fprintf (stdout,"\n reduced system of equations is solved by direct method"); 00092 xfscanf (in,"%d",(int*)&rsmstor); 00093 00094 if (ssle==NULL) 00095 ssle = new slesolv (); 00096 00097 ssle->read (top,in,mespr); 00098 00099 break; 00100 } 00101 case paral_sol:{ 00102 if (mespr==1) fprintf (stdout,"\n reduced system of equations is solved by the iterative method"); 00103 xfscanf (in,"%ld %le",&nicg,&errcg); 00104 break; 00105 } 00106 default:{ 00107 print_err("unknown reduced system solver is required",__FILE__,__LINE__,__func__); 00108 } 00109 } 00110 00111 } 00112 00113 /** 00114 function prints basic data 00115 00116 @param out - output file 00117 00118 JK, 13.5.2009 00119 */ 00120 void seqschur::print (FILE *out) 00121 { 00122 // number of subdomains 00123 fprintf (out,"%ld\n",ns); 00124 // type of reduced solver 00125 fprintf (out,"%d\n",trssol); 00126 00127 switch (trssol){ 00128 case master_sol:{ 00129 fprintf (out,"%d",rsmstor); 00130 ssle->print (out); 00131 break; 00132 } 00133 case paral_sol:{ 00134 fprintf (out,"%ld %le\n",nicg,errcg); 00135 break; 00136 } 00137 default:{ 00138 print_err("unknown reduced system solver is required",__FILE__,__LINE__,__func__); 00139 } 00140 } 00141 00142 } 00143 00144 00145 00146 00147 /** 00148 function generates ordering of unknowns on subdomains 00149 00150 @param top - pointer to the sequential general topology 00151 @param ltg - local to global array 00152 @param out - output stream 00153 00154 JK, 13.5.2009 00155 */ 00156 /* 00157 void seqschur::subdomain_ordering (gtopology *top,long *ltg,FILE *out) 00158 { 00159 long i,j,k,ii,ndofn,*aux,*idgnn,**gnnc, nbn; 00160 00161 // searching of maximum code number 00162 ndof=0; 00163 for (i=0;i<top->nn;i++){ 00164 ndofn=top->give_ndofn (i); 00165 for (j=0;j<ndofn;j++){ 00166 k=top->give_dof (i,j); 00167 if (k>ndof) ndof=k; 00168 } 00169 } 00170 ndof--; 00171 if (ndof<0) ndof=0; 00172 aux = new long [ndof]; 00173 for (i=0;i<ndof;i++){ 00174 aux[i]=-1; 00175 } 00176 00177 gnnc = new long*[top->nn]; 00178 memset(gnnc, 0, sizeof(*gnnc)*top->nn); 00179 idgnn = new long[top->nn]; 00180 memset(idgnn, 0, sizeof(*idgnn)*top->nn); 00181 j = 0; 00182 // array with indices of boundary nodes 00183 for (i=0;i<top->nn;i++){ 00184 if (ltg[i]>-1) 00185 { 00186 idgnn[j]=ltg[i]; 00187 j++; 00188 } 00189 } 00190 nbn = j; 00191 for (i=0;i<top->nn;i++){ 00192 ndofn=top->give_ndofn (i); 00193 gnnc[i] = new long[ndofn]; 00194 for (j=0;j<ndofn;j++){ 00195 gnnc[i][j]=-1; 00196 } 00197 } 00198 // contributions from nodes 00199 ndof=1; 00200 for (i=0;i<top->nn;i++){ 00201 ii=top->ordering[i]; 00202 if (ltg[ii]==-1){ 00203 // inner nodes 00204 ndofn=top->give_ndofn (ii); 00205 for (j=0;j<ndofn;j++){ 00206 k=top->give_dof(ii,j); 00207 if (k<0) continue; 00208 if (k==0) continue; 00209 if (k==1){ 00210 top->save_dof (ii,j,ndof); ndof++; 00211 } 00212 // common code numbers 00213 if (k>1){ 00214 // node has not assigned code number 00215 if ((aux[k-2]==-1) && (gnnc[ii][j]==-1)) 00216 { 00217 // search for the same common code number on the boundary nodes 00218 search_comcn(top,ltg,idgnn,nbn,gnnc,ii,k); 00219 } 00220 if ((aux[k-2]==-1) && (gnnc[ii][j]==-1)) 00221 { 00222 // node has not assigned common code number and it is not coupled with the domain boundary 00223 top->save_dof (ii,j,ndof); 00224 aux[k-2]=ndof; 00225 ndof++; 00226 } 00227 if ((aux[k-2]>-1) && (gnnc[ii][j]==-1)) 00228 { 00229 // node has already assigned common code number and it is not coupled with the domain boundary 00230 top->save_dof (ii,j,aux[k-2]); 00231 } 00232 } 00233 } 00234 } 00235 } 00236 nidof=ndof-1; 00237 top->nidof=nidof; 00238 00239 for (i=0;i<top->nn;i++){ 00240 ii=top->ordering[i]; 00241 if (ltg[ii]>-1){ 00242 // boundary nodes 00243 ndofn=top->give_ndofn (ii); 00244 for (j=0;j<ndofn;j++){ 00245 k=top->give_dof (ii,j); 00246 if (k<0) continue; 00247 if (k==0) continue; 00248 if (k==1){ 00249 top->save_dof (ii,j,ndof); ndof++; 00250 } 00251 // common code numbers 00252 if (k>1) 00253 { 00254 if (aux[k-2]==-1){ 00255 // node has not assigned common code number 00256 top->save_dof (ii,j,ndof); 00257 aux[k-2]=ndof; 00258 ndof++; 00259 } 00260 else{ 00261 // node has already assigned common code number 00262 top->save_dof (ii,j,aux[k-2]); 00263 } 00264 } 00265 } 00266 } 00267 } 00268 // numbering of nodes which are not on boundary and they have not assigned 00269 // common code number with the boundary nodes 00270 for (i=0;i<top->nn;i++){ 00271 ii=top->ordering[i]; 00272 ndofn=top->give_ndofn (ii); 00273 for (j=0;j<ndofn;j++){ 00274 k=top->give_dof(ii,j); 00275 if ((k>1) && (gnnc[ii][j]!=-1) && (ltg[ii]==-1)){ 00276 // node has already assigned common code number from boundary node and it does not lie on boundary 00277 top->save_dof (ii,j,aux[k-2]); 00278 } 00279 } 00280 } 00281 ndof--; 00282 00283 nbdof=ndof-nidof; 00284 top->nbdof=nbdof; 00285 00286 delete [] aux; 00287 delete [] idgnn; 00288 for (i=0;i<top->nn;i++) 00289 delete [] gnnc[i]; 00290 delete [] gnnc; 00291 00292 fprintf (out,"\n\n\n seqschur::subdomain_ordering \n\n"); 00293 for (i=0;i<top->nn;i++){ 00294 fprintf (out,"\n node %6ld ",i); 00295 ndofn=top->give_ndofn (i); 00296 for (j=0;j<ndofn;j++){ 00297 fprintf (out," %ld",top->give_dof(i,j)); 00298 } 00299 } 00300 fprintf (out,"\n\n"); 00301 } 00302 */ 00303 00304 00305 00306 /** 00307 function orders unknowns on subdomains with respect of Schur complement method 00308 inner variables are ordered at the beginning 00309 boundary variables are ordered at he end 00310 00311 @param top - pointer to domain topology 00312 @param ptop - pointer to topology for paralle computation 00313 00314 JK, 13.5.2009 00315 */ 00316 /* 00317 void seqschur::subdomain_ordering (gtopology *top,partop *ptop) 00318 { 00319 long i,j,k,nn,nbn,nin,ndofn,nid; 00320 00321 // number of all nodes on subdomain 00322 nn = ptop->nn; 00323 // number of internal nodes on subdomain 00324 nin = ptop->nin; 00325 // number of boundary nodes on subdomain 00326 nbn = ptop->nbn; 00327 00328 00329 // ******************************** 00330 // ordering of internal unknowns 00331 // ******************************** 00332 ndof=1; 00333 for (i=0;i<nin;i++){ 00334 // node id (in local ordering) 00335 nid=ptop->lnin[i]; 00336 // number of DOFs at required node 00337 ndofn = top->give_ndofn (nid); 00338 00339 for (j=0;j<ndofn;j++){ 00340 k=top->give_dof (nid,j); 00341 if (k>0){ 00342 top->save_dof (nid,j,ndof); 00343 ndof++; 00344 } 00345 } 00346 } 00347 // number of internal DOFs (unknowns) 00348 nidof = ndof-1; 00349 00350 00351 // ******************************** 00352 // ordering of boundary unknowns 00353 // ******************************** 00354 for (i=0;i<nbn;i++){ 00355 // node id (in local ordering) 00356 nid=ptop->lnbn[i]; 00357 // number of DOFs at required node 00358 ndofn = top->give_ndofn (nid); 00359 00360 for (j=0;j<ndofn;j++){ 00361 k=top->give_dof (nid,j); 00362 if (k>0){ 00363 top->save_dof (nid,j,ndof); 00364 ndof++; 00365 } 00366 } 00367 } 00368 ndof--; 00369 00370 00371 // contributions from elements 00372 for (i=0;i<top->ne;i++){ 00373 if (top->give_cne (i)==1){ 00374 for (j=0;j<top->give_ndofe(i);j++){ 00375 if (top->give_dof (i,j)>ndof) ndof=top->give_dof (i,j); 00376 } 00377 } 00378 } 00379 00380 // number of boundary DOFs (unknowns) 00381 nbdof = ndof-nidof; 00382 00383 // state of code numbers is changed 00384 // code numbers are generated 00385 top->cnstate=1; 00386 00387 } 00388 */ 00389 00390 00391 /** 00392 function initializes variables and arrays 00393 00394 namely: 00395 00396 ndofcp - number of coarse DOFs 00397 maxnbdof - maximum number of boundary/interface DOFs on subdomain 00398 nbdofmas - array of numbers of boundary/interface DOFs on subdomains 00399 00400 @param selnodschur - pointer to the object of the class seqselnodes 00401 00402 JK, 13.5.2009 00403 */ 00404 void seqschur::initiate (seqselnodes *selnodschur) 00405 { 00406 long i; 00407 00408 // number of coarse DOFs = total number of boundary/interface DOFs 00409 // number of DOFs of the coarse problem 00410 ndofcp = selnodschur->tndofsn; 00411 00412 // number of boundary/interface DOFs on subdomains 00413 nbdofmas = new long [ns]; 00414 00415 // array containing numbers of boundary/interface DOFs on subdomains 00416 // maximum number of reduced DOFs on subdomain 00417 maxnbdof=0; 00418 // loop over the number of subdomains 00419 for (i=0;i<ns;i++){ 00420 nbdofmas[i] = selnodschur->snndofmas[i]; 00421 if (maxnbdof<nbdofmas[i]){ 00422 maxnbdof=nbdofmas[i]; 00423 } 00424 } 00425 00426 // topology describing subdomains as superelements 00427 gtop = new gtopology; 00428 gtop->initiate (selnodschur->cndofmas,selnodschur->snndofmas,ns); 00429 gtop->cnstate=1; 00430 00431 } 00432 00433 /** 00434 function orders unknowns in coarse problem 00435 00436 function initializates object gtop of the class gtopology 00437 it represents subdomains as superelements 00438 00439 @param out - output file 00440 00441 JK, 8.7.2005 00442 */ 00443 /* 00444 void seqschur::coarse_problem_ordering (FILE *out) 00445 { 00446 00447 long i,j,k,m,tnbn; 00448 long *nbnd,*ndofncn; 00449 long **cnbn,**nbdofnd,**cncn,**bcn; 00450 long ***lbcn; 00451 00452 // array of numbers of boundary nodes on subdomains 00453 nbnd = gtop->stop->nbnd; 00454 // array containing coarse node numbers of boundary nodes 00455 cnbn = gtop->stop->lgnbn; 00456 // array containing numbers of DOFs on boundary nodes 00457 nbdofnd = gtop->stop->lbndofn; 00458 // total number of boundary nodes 00459 tnbn = gtop->stop->tnbn; 00460 // array containing local boundary code numbers 00461 lbcn = gtop->stop->lbcn; 00462 00463 00464 // number of DOFs at coarse nodes 00465 ndofncn = new long [tnbn]; 00466 00467 for (i=0;i<ns;i++){ 00468 for (j=0;j<nbnd[i];j++){ 00469 k=cnbn[i][j]; 00470 ndofncn[k]=nbdofnd[i][j]; 00471 } 00472 } 00473 00474 00475 // assignment of constraint indicators 00476 cncn = new long* [tnbn]; 00477 for (i=0;i<tnbn;i++){ 00478 cncn[i] = new long [ndofncn[i]]; 00479 } 00480 00481 00482 // number of boundary DOFs on subdomains 00483 // it must not be deleted, it is used in further subroutines 00484 nbdofmas = new long [ns]; 00485 00486 maxnbdof=0; 00487 for (i=0;i<ns;i++){ 00488 nbdofmas[i]=0; 00489 for (j=0;j<nbnd[i];j++){ 00490 m=cnbn[i][j]; 00491 for (k=0;k<nbdofnd[i][j];k++){ 00492 cncn[m][k]=lbcn[i][j][k]; 00493 if (lbcn[i][j][k]>0) 00494 nbdofmas[i]++; 00495 } 00496 } 00497 if (maxnbdof<nbdofmas[i]) 00498 maxnbdof=nbdofmas[i]; 00499 } 00500 00501 00502 // generation of code numbers (unknown ordering) 00503 ndofcp=1; 00504 for (i=0;i<tnbn;i++){ 00505 for (j=0;j<ndofncn[i];j++){ 00506 if (cncn[i][j]>0){ 00507 cncn[i][j]=ndofcp; 00508 ndofcp++; 00509 } 00510 } 00511 } 00512 ndofcp--; 00513 00514 00515 bcn = new long* [ns]; 00516 for (i=0;i<ns;i++){ 00517 bcn[i] = new long [nbdofmas[i]]; 00518 nbdofmas[i]=0; 00519 } 00520 00521 // assembling of lists of boundary code numbers on subdomains 00522 for (i=0;i<ns;i++){ 00523 for (j=0;j<nbnd[i];j++){ 00524 m=cnbn[i][j]; 00525 for (k=0;k<ndofncn[m];k++){ 00526 if (cncn[m][k]>0){ 00527 bcn[i][nbdofmas[i]]=cncn[m][k]; 00528 nbdofmas[i]++; 00529 } 00530 } 00531 } 00532 } 00533 00534 fprintf (out,"\n\n\n kontrola bcn \n"); 00535 for (i=0;i<ns;i++){ 00536 fprintf (out,"\n subdomain number %ld",i+1); 00537 for (j=0;j<nbdofmas[i];j++){ 00538 fprintf (out,"\n %ld %ld",j+1,bcn[i][j]); 00539 } 00540 } 00541 00542 // nasledujici cast jiz byla jednou volana 00543 //gtop = new gtopology; 00544 //gtop->initiate (bcn,nbdofmas,ns); 00545 //gtop->cnstate=1; 00546 00547 00548 00549 for (i=0;i<ns;i++){ 00550 delete [] bcn[i]; 00551 } 00552 delete [] bcn; 00553 00554 for (i=0;i<tnbn;i++){ 00555 delete [] cncn[i]; 00556 } 00557 delete [] cncn; 00558 00559 delete [] ndofncn; 00560 00561 } 00562 */ 00563 00564 /** 00565 function assembles list of unknowns belonging to the subdomains 00566 00567 ndofdom - numbers of DOFs on subdomains 00568 cndom - list of DOFs on subdomains 00569 00570 @param top - pointer to topology 00571 @param out - output file 00572 00573 JK, 21.5.2009 00574 */ 00575 void seqschur::assemble_subdom_unknowns (gtopology *top,FILE *out) 00576 { 00577 long i,j,l,k,g,ndofn,ii; 00578 long maxdof; 00579 long **aux; 00580 00581 aux = new long* [ns]; 00582 00583 // number of DOFs on subdomains 00584 ndofdom = new long [ns]; 00585 for (i=0;i<ns;i++){ 00586 ndofdom[i]=0; 00587 } 00588 00589 // loop over the number of subdomains 00590 for (i=0;i<ns;i++){ 00591 // maximum DOF id 00592 maxdof=0; 00593 00594 // loop over the number of internal nodes on subdomains 00595 for (j=0;j<top->stop->nind[i];j++){ 00596 // global glued number of the node 00597 g=top->stop->ggnin[i][j]; 00598 ndofn=top->give_ndofn (g); 00599 // loop over the number of DOFs 00600 for (k=0;k<ndofn;k++){ 00601 l=top->give_dof (g,k); 00602 if (l>maxdof) 00603 maxdof=l; 00604 } 00605 } 00606 // loop over the number of boundary/interface nodes on subdomains 00607 for (j=0;j<top->stop->nbnd[i];j++){ 00608 // global glued number of the node 00609 g=top->stop->ggnbn[i][j]; 00610 ndofn=top->give_ndofn (g); 00611 // loop over the number of DOFs 00612 for (k=0;k<ndofn;k++){ 00613 l=top->give_dof (g,k); 00614 if (l>maxdof) 00615 maxdof=l; 00616 } 00617 } 00618 00619 00620 aux[i] = new long [maxdof]; 00621 for (j=0;j<maxdof;j++){ 00622 aux[i][j]=-1; 00623 } 00624 00625 ii=0; 00626 // loop over the number of internal nodes on subdomains 00627 for (j=0;j<top->stop->nind[i];j++){ 00628 // global glued number of the node 00629 g=top->stop->ggnin[i][j]; 00630 ndofn=top->give_ndofn (g); 00631 // loop over the number of DOFs 00632 for (k=0;k<ndofn;k++){ 00633 l=top->give_dof (g,k); 00634 if (l>0){ 00635 if (aux[i][l-1]==-1){ 00636 aux[i][l-1]=ii; 00637 ii++; 00638 } 00639 } 00640 } 00641 } 00642 // loop over the number of boundary/interface nodes on subdomains 00643 for (j=0;j<top->stop->nbnd[i];j++){ 00644 // global glued number of the node 00645 g=top->stop->ggnbn[i][j]; 00646 ndofn=top->give_ndofn (g); 00647 // loop over the number of DOFs 00648 for (k=0;k<ndofn;k++){ 00649 l=top->give_dof (g,k); 00650 if (l>0){ 00651 if (aux[i][l-1]==-1){ 00652 aux[i][l-1]=ii; 00653 ii++; 00654 } 00655 } 00656 } 00657 } 00658 00659 ndofdom[i]=0; 00660 for (j=0;j<maxdof;j++){ 00661 if (aux[i][j]>-1){ 00662 ndofdom[i]++; 00663 } 00664 } 00665 } 00666 00667 // numbers of DOFs on subdomains are known now 00668 00669 // array of code numbers on subdomains 00670 cndom = new long* [ns]; 00671 for (i=0;i<ns;i++){ 00672 cndom[i] = new long [ndofdom[i]]; 00673 //ndofdom[i]=0; 00674 } 00675 00676 // loop over the number of subdomains 00677 for (i=0;i<ns;i++){ 00678 // loop over the number of internal nodes on subdomains 00679 for (j=0;j<top->stop->nind[i];j++){ 00680 // global glued number of the node 00681 g=top->stop->ggnin[i][j]; 00682 ndofn=top->give_ndofn (g); 00683 // loop over the number of DOFs 00684 for (k=0;k<ndofn;k++){ 00685 l=top->give_dof (g,k); 00686 if (l>0){ 00687 cndom[i][aux[i][l-1]]=l; 00688 //ndofdom[i]++; 00689 } 00690 } 00691 } 00692 // loop over the number of boundary/interface nodes on subdomains 00693 for (j=0;j<top->stop->nbnd[i];j++){ 00694 // global glued number of the node 00695 g=top->stop->ggnbn[i][j]; 00696 ndofn=top->give_ndofn (g); 00697 // loop over the number of DOFs 00698 for (k=0;k<ndofn;k++){ 00699 l=top->give_dof (g,k); 00700 if (l>0){ 00701 cndom[i][aux[i][l-1]]=l; 00702 //ndofdom[i]++; 00703 } 00704 } 00705 } 00706 } 00707 00708 00709 for (i=0;i<ns;i++){ 00710 delete [] aux[i]; 00711 } 00712 delete [] aux; 00713 00714 00715 // ******************* 00716 // auxiliary output 00717 // ******************* 00718 fprintf (out,"\n\n\n array ndofdom - number of DOFs on subdomains \n\n"); 00719 for (i=0;i<ns;i++){ 00720 fprintf (out,"\n subdomain %3ld %6ld",i,ndofdom[i]); 00721 } 00722 00723 fprintf (out,"\n\n\n array cndom - array of code numbers of DOFs on subdomains \n\n"); 00724 for (i=0;i<ns;i++){ 00725 fprintf (out,"\n\n subdomain %ld",i); 00726 for (j=0;j<ndofdom[i];j++){ 00727 fprintf (out,"\n DOF %6ld %ld",j,cndom[i][j]); 00728 } 00729 } 00730 00731 } 00732 00733 00734 /** 00735 function assembles subdomain matrices from the %matrix of the whole system 00736 00737 @param top - general topology 00738 @param gm - pointer to the %matrix of the system 00739 @param out - output file 00740 00741 JK, 13.5.2009 00742 */ 00743 void seqschur::subdomain_matrices (gmatrix *gm,FILE *out) 00744 { 00745 long i,j; 00746 /* 00747 smsky = new skyline [ns]; 00748 for (i=0;i<ns;i++){ 00749 gndofe = gtop->give_ndofe (i); 00750 cn = new long [gndofe]; 00751 gtop->give_code_numbers (i,cn); 00752 smsky[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,nbdofmas[i],cn); 00753 delete [] cn; 00754 } 00755 */ 00756 smsky = new skyline [ns]; 00757 for (i=0;i<ns;i++){ 00758 smsky[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofdom[i],cndom[i]); 00759 } 00760 00761 // ******************* 00762 // auxiliary output 00763 // ******************* 00764 00765 00766 fprintf (out,"\n\n kontrola adres ve skyline \n"); 00767 for (i=0;i<ns;i++){ 00768 fprintf (out,"\n\n podoblast %ld",i); 00769 for (j=0;j<smsky[i].n+1;j++){ 00770 fprintf (out,"\n adr %4ld %ld",j,smsky[i].adr[j]); 00771 } 00772 } 00773 00774 fprintf (out,"\n\n kontrola ndofdom\n"); 00775 for (i=0;i<ns;i++){ 00776 fprintf (out,"\n pocet DOFs na domene %ld %ld\n",i,ndofdom[i]); 00777 for (j=0;j<ndofdom[i];j++){ 00778 fprintf (out," %ld",cndom[i][j]); 00779 } 00780 } 00781 00782 fprintf (out,"\n\n\n"); 00783 for (i=0;i<smsky[1].n;i++){ 00784 fprintf (out,"\n\n%ld\n",i); 00785 for (j=0;j<i;j++){ 00786 if (smsky[1].adr[i]+i-j>=smsky[1].adr[i+1]) fprintf (out," %10.5f",0.0); 00787 else fprintf (out," %10.5f",smsky[1].a[smsky[1].adr[i]+i-j]); 00788 } 00789 for (j=i;j<smsky[1].n;j++){ 00790 if (smsky[1].adr[j]+j-i>=smsky[1].adr[j+1]) fprintf (out," %10.5f",0.0); 00791 else fprintf (out," %10.5f",smsky[1].a[smsky[1].adr[j]+j-i]); 00792 } 00793 } 00794 00795 } 00796 00797 00798 00799 00800 00801 /** 00802 function solves reduced system of equations from primal 00803 domain decomposition by conjugate gradient method 00804 %matrix of reduced system is not assembled 00805 iteration process is directed by the master processor 00806 00807 @param condmat - array containing reduced matrices (Schur complements) 00808 @param condvect - array containing reduced vectors 00809 @param out - output file (for auxiliary output) 00810 00811 JK, 13.5.2009 00812 */ 00813 void seqschur::solve_red_sys_iter (double **condmat,double **condvect,FILE *out) 00814 { 00815 /* 00816 int ind; 00817 long i,j,k,sizebuff,gndofe,*cn; 00818 double alpha,beta,nom,denom,norrhs; 00819 double *lhs,*rhs,*buff,*p,*d,*r; 00820 00821 rhs=NULL; 00822 lhs=NULL; 00823 sizebuff = maxnbdof+1; 00824 buff = new double [sizebuff]; 00825 memset (buff,0,sizebuff*sizeof(double)); 00826 00827 // right hand side assembling 00828 if (myrank==0){ 00829 lhs = new double [ndofcp]; 00830 memset (lhs,0,ndofcp*sizeof(double)); 00831 rhs = new double [ndofcp]; 00832 memset (rhs,0,ndofcp*sizeof(double)); 00833 00834 gndofe = gtop->give_ndofe (domproc[0]); 00835 cn = new long [gndofe]; 00836 gtop->give_code_numbers (domproc[0],cn); 00837 locglob (rhs,condvect,cn,gndofe); 00838 delete [] cn; 00839 00840 for (i=1;i<nproc;i++){ 00841 MPI_Recv(buff,sizebuff,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 00842 j=domproc[stat.MPI_TAG]; ind=0; 00843 00844 gndofe = gtop->give_ndofe (j); 00845 cn = new long [gndofe]; 00846 gtop->give_code_numbers (j,cn); 00847 locglob (rhs,buff,cn,gndofe); 00848 delete [] cn; 00849 00850 } 00851 } 00852 else{ 00853 copyv (condvect,buff,ndof-nidof); 00854 MPI_Send(buff,sizebuff,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD); 00855 } 00856 00857 // iteration loop block 00858 if (myrank==0){ 00859 r = new double [ndofcp]; 00860 d = new double [ndofcp]; 00861 p = new double [ndofcp]; 00862 00863 // iteration initialization 00864 norrhs=0.0; nom=0.0; 00865 for (i=0;i<ndofcp;i++){ 00866 lhs[i]=0.0; 00867 r[i]=rhs[i]; 00868 d[i]=rhs[i]; 00869 norrhs+=rhs[i]*rhs[i]; 00870 nom+=r[i]*r[i]; 00871 } 00872 00873 // iteration loop 00874 for (i=0;i<nicgsch;i++){ 00875 00876 // direction vector scattering 00877 for (j=1;j<nproc;j++){ 00878 memset (buff,0,sizebuff*sizeof(double)); 00879 gndofe = gtop->give_ndofe (domproc[j]); 00880 cn = new long [gndofe]; 00881 gtop->give_code_numbers (domproc[j],cn); 00882 globloc (d,buff,cn,gndofe); 00883 delete [] cn; 00884 00885 MPI_Send(buff,sizebuff,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD); 00886 } 00887 00888 memset (buff,0,sizebuff*sizeof(double)); 00889 gndofe = gtop->give_ndofe (domproc[0]); 00890 cn = new long [gndofe]; 00891 gtop->give_code_numbers (domproc[0],cn); 00892 globloc (d,buff,cn,gndofe); 00893 delete [] cn; 00894 00895 // matrix-vector multiplication 00896 mxv (condmat,buff,condvect,nbdof,nbdof); 00897 00898 // matrix-vector multiplication collection 00899 memset (p,0,ndofcp*sizeof(double)); 00900 00901 gndofe = gtop->give_ndofe (domproc[0]); 00902 cn = new long [gndofe]; 00903 gtop->give_code_numbers (domproc[0],cn); 00904 locglob (p,condvect,cn,gndofe); 00905 delete [] cn; 00906 00907 for (j=1;j<nproc;j++){ 00908 MPI_Recv(buff,sizebuff,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 00909 k=domproc[stat.MPI_TAG]; 00910 00911 gndofe = gtop->give_ndofe (k); 00912 cn = new long [gndofe]; 00913 gtop->give_code_numbers (k,cn); 00914 locglob (p,buff,cn,gndofe); 00915 delete [] cn; 00916 00917 } 00918 00919 // denominator 00920 denom=ss(d,p,ndofcp); 00921 00922 if (fabs(denom)<zero){ 00923 if (i!=nicgsch-1){ 00924 buff[maxnbdof]=1.0; 00925 for (j=1;j<nproc;j++){ 00926 MPI_Send(buff,sizebuff,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD); 00927 } 00928 } 00929 fprintf (stderr,"\n\n denominator in alpha expression in parallel conjugate gradient method is equal to zero.\n"); 00930 break; 00931 } 00932 00933 // coefficient alpha 00934 alpha=nom/denom; 00935 00936 // new global vectors 00937 for (j=0;j<ndofcp;j++){ 00938 r[j]-=alpha*p[j]; 00939 lhs[j]+=alpha*d[j]; 00940 } 00941 00942 00943 denom=nom; 00944 nom=ss(r,r,ndofcp); 00945 00946 if (i%100==0) fprintf (stdout,"\n iteration %ld norres/norrhs %e",i,nom/norrhs); 00947 //fprintf (out,"\n iteration %ld norres/norrhs %e",i,nom/norrhs); 00948 if (nom/norrhs<errcgsch){ 00949 fprintf (stdout,"\n last iteration %ld norres/norrhs %e",i,nom/norrhs); 00950 if (i!=nicgsch-1){ 00951 buff[maxnbdof]=1.0; 00952 for (j=1;j<nproc;j++){ 00953 MPI_Send(buff,sizebuff,MPI_DOUBLE,j,myrank,MPI_COMM_WORLD); 00954 } 00955 } 00956 break; 00957 } 00958 00959 // beta coefficient 00960 beta=nom/denom; 00961 00962 // new direction vector 00963 for (j=0;j<ndofcp;j++){ 00964 d[j]=r[j]+beta*d[j]; 00965 } 00966 00967 } 00968 anicgsch=i; 00969 aerrcgsch=nom/norrhs; 00970 00971 00972 //fprintf (out,"\n\n\n REDUCED SYSTEM SOLUTION (on master)\n"); 00973 //for (i=0;i<ndofcp;i++){ 00974 //fprintf (out,"\n lhs %ld %e",i,lhs[i]); 00975 //} 00976 00977 00978 } 00979 else{ 00980 // iteration loop 00981 for (i=0;i<nicgsch;i++){ 00982 // direction vector receipt 00983 MPI_Recv(buff,sizebuff,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 00984 00985 if (buff[maxnbdof]>0.5) break; 00986 00987 // matrix-vector multiplication 00988 mxv (condmat,buff,condvect,nbdof,nbdof); 00989 00990 // matrix-vector multiplication collection 00991 copyv (condvect,buff,nbdof); 00992 MPI_Send(buff,sizebuff,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD); 00993 00994 } 00995 anicgsch=i; 00996 } 00997 00998 // global solution distribution 00999 if (myrank==0){ 01000 for (i=1;i<nproc;i++){ 01001 01002 memset (buff,0,sizebuff*sizeof(double)); 01003 gndofe = gtop->give_ndofe (domproc[i]); 01004 cn = new long [gndofe]; 01005 gtop->give_code_numbers (domproc[i],cn); 01006 globloc (lhs,buff,cn,gndofe); 01007 delete [] cn; 01008 01009 MPI_Send(buff,sizebuff,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD); 01010 } 01011 01012 nullv (condvect,nbdofmas[domproc[0]]); 01013 01014 gndofe = gtop->give_ndofe (domproc[0]); 01015 cn = new long [gndofe]; 01016 gtop->give_code_numbers (domproc[0],cn); 01017 globloc (lhs,buff,cn,gndofe); 01018 delete [] cn; 01019 01020 } 01021 else{ 01022 MPI_Recv(buff,sizebuff,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 01023 } 01024 copyv (buff,condvect,nbdof); 01025 01026 delete [] buff; 01027 */ 01028 } 01029 01030 01031 /** 01032 function solves reduced system of linear algebraic equations on the master 01033 processor by a direct method 01034 01035 @param condmat - array containing reduced matrices (Schur complements) 01036 @param condvect - array containing reduced vectors 01037 @param out - output file (for auxiliary output) 01038 01039 JK, 13.5.2009 01040 */ 01041 void seqschur::solve_red_sys_fin (double **condmat,double **condvect,FILE *out) 01042 { 01043 long i,gndofe,*cn; 01044 double *lhs,*rhs; 01045 precond prec; 01046 01047 arr = new gmatrix; 01048 01049 //arr->ts = rsmstor; 01050 arr->setval (ssle); 01051 //arr->tlinsol=ldl; 01052 01053 01054 //arr->alloc (); 01055 arr->initiate (gtop,ndofcp,rsmstor,1,out); 01056 01057 lhs = new double [ndofcp]; 01058 memset (lhs,0,ndofcp*sizeof(double)); 01059 rhs = new double [ndofcp]; 01060 memset (rhs,0,ndofcp*sizeof(double)); 01061 01062 for (i=0;i<ns;i++){ 01063 gndofe = gtop->give_ndofe (i); 01064 cn = new long [gndofe]; 01065 gtop->give_code_numbers (i,cn); 01066 arr->localized (condmat[i],cn,gndofe,i); 01067 locglob (rhs,condvect[i],cn,gndofe); 01068 delete [] cn; 01069 } 01070 01071 arr->prepmat (0.0,1); 01072 01073 01074 // solution of the reduced problem 01075 arr->solve_system (gtop,prec,lhs,rhs,out); 01076 01077 01078 fprintf (out,"\n\n\n REDUCED SYSTEM SOLUTION (on master)\n"); 01079 for (i=0;i<ndofcp;i++){ 01080 fprintf (out,"\n lhs %ld %e",i,lhs[i]); 01081 } 01082 01083 for (i=0;i<ns;i++){ 01084 nullv (condvect[i],nbdofmas[i]); 01085 01086 gndofe = gtop->give_ndofe (i); 01087 cn = new long [gndofe]; 01088 gtop->give_code_numbers (i,cn); 01089 globloc (lhs,condvect[i],cn,gndofe); 01090 delete [] cn; 01091 } 01092 01093 01094 delete [] lhs; delete [] rhs; 01095 delete arr; 01096 01097 } 01098 01099 /** 01100 function solves reduced system of linear algebraic equations 01101 01102 @param condmat - array containing reduced matrices (Schur complements) 01103 @param condvect - array containing reduced vectors 01104 @param out - output file (for auxiliary output) 01105 01106 JK, 13.5.2009 01107 */ 01108 void seqschur::solve_red_sys (double **condmat,double **condvect,FILE *out) 01109 { 01110 switch (trssol){ 01111 case master_sol:{ 01112 solve_red_sys_fin (condmat,condvect,out); 01113 break; 01114 } 01115 case paral_sol:{ 01116 solve_red_sys_iter (condmat,condvect,out); 01117 break; 01118 } 01119 default:{ 01120 print_err ("unknown type of solver of reduced problem is required",__FILE__,__LINE__,__func__); 01121 } 01122 } 01123 } 01124 01125 01126 /** 01127 function solves system of linear algebraic equations 01128 by the Schur complment method 01129 01130 @param top - pointer to the sequential general topology 01131 @param gm - pointer to the subdomain %matrix stored in the form gmatrix 01132 @param lhs - array containing solution 01133 @param rhs - array containing right hand side 01134 @param out - output file (for auxiliary output) 01135 01136 JK, 13.5.2009 01137 */ 01138 void seqschur::solve_system (gtopology *top,gmatrix *gm, 01139 double *lhs,double *rhs,FILE *out) 01140 { 01141 long i,j; 01142 double **condmat,**condvect; 01143 time_t t1,t2,t3,t4; 01144 01145 // allocation of reduced matrix 01146 condmat = new double* [ns]; 01147 for (i=0;i<ns;i++){ 01148 condmat[i] = new double [nbdofmas[i]*nbdofmas[i]]; 01149 for (j=0;j<nbdofmas[i]*nbdofmas[i];j++){ 01150 condmat[i][j]=0.0; 01151 } 01152 } 01153 01154 // allocation of reduced vector 01155 condvect = new double* [ns]; 01156 for (i=0;i<ns;i++){ 01157 condvect[i] = new double [nbdofmas[i]]; 01158 for (j=0;j<nbdofmas[i];j++){ 01159 condvect[i][j]=0.0; 01160 } 01161 } 01162 01163 01164 // assembling of subdomain matrices 01165 subdomain_matrices (gm,out); 01166 01167 // allocation of vectors of right hand sides on subdomains 01168 double **ff; 01169 ff = new double* [ns]; 01170 for (i=0;i<ns;i++){ 01171 ff[i] = new double [ndofdom[i]]; 01172 for (j=0;j<ndofdom[i];j++){ 01173 ff[i][j]=0.0; 01174 } 01175 } 01176 01177 01178 // allocation of vectors of right hand sides on subdomains 01179 for (i=0;i<ns;i++){ 01180 globloc (rhs,ff[i],cndom[i],ndofdom[i]); 01181 } 01182 01183 double **dd; 01184 dd = new double* [ns]; 01185 for (i=0;i<ns;i++){ 01186 dd[i] = new double [ndofdom[i]]; 01187 for (j=0;j<ndofdom[i];j++){ 01188 dd[i][j]=0.0; 01189 } 01190 } 01191 01192 t1 = time (NULL); 01193 01194 // elimination of inner DOFs = static condensation of inner DOFs 01195 // the Schur complement matrix is stored in array condmat 01196 for (i=0;i<ns;i++){ 01197 smsky[i].ldlkon_sky (condmat[i],condvect[i],dd[i],ff[i],nbdofmas[i],1); 01198 } 01199 01200 t2 = time (NULL); 01201 01202 // solution of reduced problem 01203 solve_red_sys (condmat,condvect,out); 01204 01205 t3 = time (NULL); 01206 01207 01208 // back substitution on subdomains 01209 for (i=0;i<ns;i++){ 01210 smsky[i].ldlkon_sky (condmat[i],condvect[i],dd[i],ff[i],nbdofmas[i],2); 01211 } 01212 01213 for (i=0;i<ns;i++){ 01214 for (j=0;j<ndofdom[i];j++){ 01215 lhs[cndom[i][j]-1]=dd[i][j]; 01216 } 01217 } 01218 /* 01219 // allocation of vectors of right hand sides on subdomains 01220 for (i=0;i<ns;i++){ 01221 gndofe = gtop->give_ndofe (i); 01222 cn = new long [gndofe]; 01223 gtop->give_code_numbers (i,cn); 01224 locglob (lhs,dd[i],cn,gndofe); 01225 delete [] cn; 01226 } 01227 */ 01228 01229 t4 = time (NULL); 01230 01231 // **************************************************** 01232 // gathering of elapsed time on the master processor 01233 // **************************************************** 01234 long *buff; 01235 buff = new long[3]; 01236 buff[0]=t2-t1; 01237 buff[1]=t3-t2; 01238 buff[2]=t4-t3; 01239 01240 /* 01241 if (myrank==0){ 01242 01243 long *fact,*cp,*tottime; 01244 fact=new long [nproc]; 01245 cp=new long [nproc]; 01246 tottime = new long [nproc]; 01247 01248 fprintf (out,"\n\n\n\n\n"); 01249 j=domproc[myrank]; 01250 fact[j]=buff[0]; 01251 cp[j]=buff[1]; 01252 tottime[j]=buff[0]+buff[1]+buff[2]; 01253 01254 fprintf (out,"\n\n\n Domain %ld",j); 01255 fprintf (out,"\n time of condensation %ld",buff[0]); 01256 fprintf (out,"\n time of solution of red. sys. %ld",buff[1]); 01257 fprintf (out,"\n time of back substitution %ld",buff[2]); 01258 01259 j=domproc[myrank]; 01260 // fprintf (stdout,"\n\n\n Domain %ld",j); 01261 // fprintf (stdout,"\n time of condensation %ld",buff[0]); 01262 // fprintf (stdout,"\n time of solution of red. sys. %ld",buff[1]); 01263 // fprintf (stdout,"\n time of back substitution %ld",buff[2]); 01264 01265 for (i=1;i<nproc;i++){ 01266 MPI_Recv (buff,3,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 01267 j=domproc[stat.MPI_TAG]; 01268 fact[j]=buff[0]; 01269 cp[j]=buff[1]; 01270 tottime[j]=buff[0]+buff[1]+buff[2]; 01271 01272 fprintf (out,"\n\n\n Domain %ld",j); 01273 fprintf (out,"\n time of condensation %ld",buff[0]); 01274 fprintf (out,"\n time of solution of red. sys. %ld",buff[1]); 01275 fprintf (out,"\n time of back substitution %ld",buff[2]); 01276 01277 // fprintf (stdout,"\n\n\n Domain %ld",j); 01278 // fprintf (stdout,"\n time of condensation %ld",buff[0]); 01279 // fprintf (stdout,"\n time of solution of red. sys. %ld",buff[1]); 01280 // fprintf (stdout,"\n time of back substitution %ld",buff[2]); 01281 } 01282 01283 long maxfact,mincp,mintot,maxtot; 01284 maxfact=0; 01285 maxtot=0; 01286 mincp=LONG_MAX; 01287 mintot=LONG_MAX; 01288 for (i=0;i<nproc;i++){ 01289 if (maxfact<fact[i]) 01290 maxfact=fact[i]; 01291 01292 if (mincp>cp[i]) 01293 mincp=cp[i]; 01294 01295 if (mintot>tottime[i]) 01296 mintot=tottime[i]; 01297 if (maxtot<tottime[i]) 01298 maxtot=tottime[i]; 01299 } 01300 01301 fprintf (out,"\n\n maximum factorization time %ld",maxfact); 01302 fprintf (out,"\n coarse problem time %ld",mincp); 01303 fprintf (out,"\n minimum total time %ld",mintot); 01304 fprintf (out,"\n maximum total time %ld",maxtot); 01305 01306 delete [] fact; 01307 delete [] cp; 01308 01309 } 01310 else{ 01311 MPI_Send(buff,3,MPI_LONG,0,myrank,MPI_COMM_WORLD); 01312 } 01313 MPI_Barrier (MPI_COMM_WORLD); 01314 */ 01315 delete [] buff; 01316 01317 01318 delete [] condmat; 01319 delete [] condvect; 01320 01321 } 01322 01323 /** 01324 function gathers contributions from subdomains into one global %vector 01325 the global %vector must be deleted outside of the subroutine 01326 01327 @param lv - local %vector 01328 @param gv - global %vector, it is allocated only on the master processor 01329 the array must be deleted outside of the subroutine 01330 01331 JK, 6.11.2004 01332 */ 01333 /* 01334 void seqschur::gather_bound_vect (double *lv,double *gv,long *domproc) 01335 { 01336 long i,j,l,gndofe,*cn; 01337 double *buff; 01338 MPI_Status stat; 01339 01340 buff = new double [maxnbdof]; 01341 01342 // values from subdomain are selected 01343 j=0; 01344 for (i=nidof;i<ndof;i++){ 01345 buff[j]=lv[i]; j++; 01346 } 01347 01348 if (myrank==0){ 01349 // ****************************** 01350 // contribution from the master 01351 // ****************************** 01352 l=domproc[0]; 01353 // number of DOFs on appropriate subdomain 01354 gndofe = gtop->give_ndofe (l); 01355 cn = new long [gndofe]; 01356 // code numbers of subdomain 01357 gtop->give_code_numbers (l,cn); 01358 // localization of values from local vector to the global vector 01359 locglob (gv,buff,cn,gndofe); 01360 delete [] cn; 01361 01362 // *************************** 01363 // contributions from slaves 01364 // *************************** 01365 for (i=1;i<nproc;i++){ 01366 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 01367 l=domproc[stat.MPI_TAG]; 01368 01369 // number of DOFs on appropriate subdomain 01370 gndofe = gtop->give_ndofe (l); 01371 cn = new long [gndofe]; 01372 // code numbers of subdomain 01373 gtop->give_code_numbers (l,cn); 01374 // localization of values from local vector to the global vector 01375 locglob (gv,buff,cn,gndofe); 01376 delete [] cn; 01377 } 01378 } 01379 else{ 01380 MPI_Send(buff,maxnbdof,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD); 01381 } 01382 MPI_Barrier (MPI_COMM_WORLD); 01383 01384 delete [] buff; 01385 } 01386 */ 01387 01388 /** 01389 function gathers contributions from subdomains into one global %vector 01390 the global %vector must be deleted outside of the subroutine 01391 01392 @param lv - local %vector 01393 @param gv - global %vector, it is allocated only on the master processor 01394 the array must be deleted outside of the subroutine 01395 01396 JK, 6.11.2004 01397 */ 01398 /* 01399 double seqschur::pss_gather_bound_vect (double *lv,double *gv,long *domproc) 01400 { 01401 long i,j,l,gndofe,*cn,buffsize; 01402 double dp,idp,*buff; 01403 MPI_Status stat; 01404 01405 buffsize=maxnbdof+1; 01406 buff = new double [buffsize]; 01407 01408 idp=0.0; 01409 01410 // dot product of internal DOFs 01411 dp=0.0; 01412 for (i=0;i<nidof;i++){ 01413 dp+=lv[i]*lv[i]; 01414 } 01415 buff[maxnbdof]=dp; 01416 01417 // interface values from subdomain are selected 01418 j=0; 01419 for (i=nidof;i<ndof;i++){ 01420 buff[j]=lv[i]; j++; 01421 } 01422 01423 if (myrank==0){ 01424 // sum of contributions from internal DOFs 01425 idp=buff[maxnbdof]; 01426 01427 // ****************************** 01428 // contribution from the master 01429 // ****************************** 01430 l=domproc[0]; 01431 // number of DOFs on appropriate subdomain 01432 gndofe = gtop->give_ndofe (l); 01433 cn = new long [gndofe]; 01434 // code numbers of subdomain 01435 gtop->give_code_numbers (l,cn); 01436 // localization of values from local vector to the global vector 01437 locglob (gv,buff,cn,gndofe); 01438 delete [] cn; 01439 01440 // *************************** 01441 // contributions from slaves 01442 // *************************** 01443 for (i=1;i<nproc;i++){ 01444 MPI_Recv(buff,buffsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 01445 l=domproc[stat.MPI_TAG]; 01446 01447 // number of DOFs on appropriate subdomain 01448 gndofe = gtop->give_ndofe (l); 01449 cn = new long [gndofe]; 01450 // code numbers of subdomain 01451 gtop->give_code_numbers (l,cn); 01452 // localization of values from local vector to the global vector 01453 locglob (gv,buff,cn,gndofe); 01454 delete [] cn; 01455 01456 // sum of contributions from internal DOFs 01457 idp+=buff[maxnbdof]; 01458 } 01459 } 01460 else{ 01461 MPI_Send(buff,buffsize,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD); 01462 } 01463 MPI_Barrier (MPI_COMM_WORLD); 01464 01465 delete [] buff; 01466 01467 return idp; 01468 } 01469 */ 01470 01471 01472 01473 01474 /** 01475 function scatters contributions from one global vector into local vectors 01476 function rewrites part of the local vector by contributions from boundaries 01477 01478 @param lv - local vector 01479 @param gv - global vector, it is allocated only on the master processor 01480 01481 JK, 6.11.2004 01482 */ 01483 /* 01484 void seqschur::scatter_bound_vect (double *lv,double *gv,long *domproc) 01485 { 01486 long i,j,gndofe,*cn; 01487 double *buff; 01488 MPI_Status stat; 01489 01490 buff = new double [maxnbdof]; 01491 01492 if (myrank==0){ 01493 for (i=1;i<nproc;i++){ 01494 nullv (buff,nbdofmas[domproc[i]]); 01495 01496 gndofe = gtop->give_ndofe (domproc[i]); 01497 cn = new long [gndofe]; 01498 gtop->give_code_numbers (domproc[i],cn); 01499 globloc (gv,buff,cn,gndofe); 01500 delete [] cn; 01501 01502 MPI_Send(buff,maxnbdof,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD); 01503 } 01504 01505 nullv (buff,nbdofmas[domproc[0]]); 01506 01507 gndofe = gtop->give_ndofe (domproc[0]); 01508 cn = new long [gndofe]; 01509 gtop->give_code_numbers (domproc[0],cn); 01510 globloc (gv,buff,cn,gndofe); 01511 delete [] cn; 01512 01513 } 01514 else{ 01515 MPI_Recv(buff,maxnbdof,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 01516 } 01517 MPI_Barrier (MPI_COMM_WORLD); 01518 01519 j=0; 01520 for (i=nidof;i<ndof;i++){ 01521 lv[i]=buff[j]/2.0; j++; 01522 } 01523 01524 delete [] buff; 01525 01526 } 01527 */ 01528 01529 /** 01530 function computes norm of vector of unbalanced values 01531 function gathers vector of unbalanced values on boundaries 01532 function scatters correct values into subdomains 01533 function returns norm of vector of unbalanced values 01534 01535 @param lv - local vector 01536 @param domproc - domain-processor correspondence 01537 01538 JK, 6.11.2004 01539 */ 01540 /* 01541 double seqschur::unbalanced_values (double *lv,long *domproc,FILE *out) 01542 { 01543 long i; 01544 double aux,norm,parlocnorm; 01545 double *gv=NULL; 01546 MPI_Status stat; 01547 01548 if (myrank==0){ 01549 gv = new double [ndofcp]; 01550 nullv (gv,ndofcp); 01551 } 01552 01553 // norm of vector of unbalanced values at internal nodes 01554 parlocnorm =0; 01555 for (i=0;i<nidof;i++){ 01556 parlocnorm+=lv[i]*lv[i]; 01557 } 01558 01559 // gathering of partial norms from subdomains 01560 if (myrank==0){ 01561 norm=parlocnorm; 01562 for (i=1;i<nproc;i++){ 01563 MPI_Recv (&aux,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 01564 norm+=aux; 01565 } 01566 } 01567 else{ 01568 MPI_Send(&parlocnorm,1,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD); 01569 } 01570 MPI_Barrier (MPI_COMM_WORLD); 01571 01572 // gathering of unbalanced values at boundary nodes 01573 gather_bound_vect (lv,gv,domproc); 01574 01575 // contribution to the norm from boundary nodes 01576 if (myrank==0){ 01577 for (i=0;i<ndofcp;i++){ 01578 norm+=gv[i]*gv[i]; 01579 //fprintf (out,"\n%6ld %e",i,gv[i]); 01580 } 01581 } 01582 01583 // scattering of unbalanced values at boundary nodes into subdomains 01584 scatter_bound_vect (lv,gv,domproc); 01585 01586 // scattering of norm of vector of unbalanced values 01587 if (myrank==0){ 01588 for (i=1;i<nproc;i++){ 01589 MPI_Send(&norm,1,MPI_DOUBLE,i,myrank,MPI_COMM_WORLD); 01590 } 01591 } 01592 else{ 01593 MPI_Recv (&norm,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat); 01594 } 01595 MPI_Barrier (MPI_COMM_WORLD); 01596 01597 if (myrank==0) 01598 delete [] gv; 01599 01600 return norm; 01601 } 01602 */ 01603 01604 01605 01606 01607 01608 /** 01609 Function search prescribed common code numbers on boundary nodes in the given subdomain. 01610 @param top - pointer to the sequential general topology 01611 @param ltg - local to global array 01612 @param idgnn - array with node numbers of boundary nodes 01613 @param nbn - number of boundary nodes 01614 @param gnnc - array with boundary node numbers at positions of nodes with prescribed common code numbers 01615 @param id - searched node id 01616 @param ccn - searched common code number id 01617 01618 */ 01619 /* 01620 void seqschur::search_comcn(gtopology *top,long *ltg, long *idgnn, long nbn, long **gnnc, long id, long ccn) 01621 { 01622 long i, j, k, ii, ndofn, flag; 01623 for (i=0; i<nbn; i++) 01624 { 01625 flag = 0; 01626 ii = idgnn[i]; 01627 ndofn=top->give_ndofn(ii); 01628 for (j=0;j<ndofn;j++){ 01629 k=top->give_dof (i,j); 01630 if (k == ccn) 01631 { 01632 i = nbn; 01633 flag = 1; 01634 break; 01635 } 01636 } 01637 } 01638 if (flag) 01639 { 01640 for (i=id;i<top->nn;i++){ 01641 for (j=0;j<ndofn;j++){ 01642 k=top->give_dof (i,j); 01643 if (k == ccn) 01644 { 01645 gnnc[i][k]=ltg[ii]; 01646 } 01647 } 01648 } 01649 } 01650 } 01651 */