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