00001 #include "aggregator.h"
00002 #include <stdlib.h>
00003 #include <limits.h>
00004
00005 aggregator::aggregator ()
00006 {
00007
00008 na=0;
00009
00010 nn=0;
00011
00012 n=0;
00013
00014 cms=0;
00015
00016
00017 impl=0;
00018
00019
00020 trbm=0;
00021
00022 nrbm=0;
00023
00024
00025 specrad=0.0;
00026
00027
00028 degree_k=0;
00029 degree_r=0;
00030
00031
00032 exinex=1;
00033
00034
00035 lnntagr=NULL;
00036 lntagr=NULL;
00037
00038 lnnagr=NULL;
00039 lnagr=NULL;
00040
00041
00042 lnuaggr=NULL;
00043
00044 luaggr=NULL;
00045
00046
00047 lmdm=NULL;
00048 lmcr=NULL;
00049 lmsky=NULL;
00050 sdirect=NULL;
00051
00052
00053 ssle = new slesolv();
00054
00055 }
00056
00057 aggregator::~aggregator ()
00058 {
00059 delete [] lnuaggr;
00060 delete [] luaggr;
00061
00062 delete [] lmdm;
00063 delete [] lmcr;
00064 delete [] lmsky;
00065 delete [] sdirect;
00066 }
00067
00068
00069
00070
00071
00072
00073
00074
00075 void aggregator::read (gtopology *gt,XFILE *in,long mespr)
00076 {
00077
00078
00079
00080 xfscanf (in,"%ld",&impl);
00081 if (impl!=1 && impl!=2)
00082 print_err("wrong type of BOSS implemenation is required",__FILE__,__LINE__,__func__);
00083
00084
00085
00086
00087 xfscanf (in,"%ld",&exinex);
00088 if (exinex!=1 && exinex!=2)
00089 print_err("wrong type of solver (exact/inexact) is required",__FILE__,__LINE__,__func__);
00090
00091
00092
00093 xfscanf (in,"%ld",°ree_k);
00094 if (mespr==1)
00095 fprintf (stdout,"\n degree of recursion in BOSS is %ld",degree_k);
00096
00097
00098
00099
00100 xfscanf (in,"%ld",&trbm);
00101
00102
00103
00104
00105 ssle->read (gt,in,mespr);
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 }
00147
00148
00149
00150
00151
00152
00153
00154
00155 void aggregator::print (FILE *out)
00156 {
00157
00158 fprintf (out,"\n%ld",impl);
00159
00160
00161 fprintf (out,"\n%ld",exinex);
00162
00163
00164 fprintf (out,"\n%ld",degree_k);
00165
00166
00167 fprintf (out,"\n%ld",trbm);
00168
00169
00170
00171
00172 fprintf (out,"\n");
00173 ssle->print (out);
00174
00175 }
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 void aggregator::define_na (long i)
00188 {
00189 na=i;
00190 }
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 void aggregator::define_graph (long nnod, long *nneigh, long **listneigh)
00204 {
00205
00206 nn = nnod;
00207
00208 nadjnodnod = nneigh;
00209
00210 adjnodnod = listneigh;
00211
00212 }
00213
00214
00215
00216
00217
00218
00219
00220 void aggregator::prepare_tlnagr(void)
00221 {
00222
00223
00224 long i,j,k;
00225 long local_n = 10*nn/na/11;
00226 long cnt;
00227 long *lnodes;
00228 long *vrstva;
00229 long vrstva_start, vrstva_konec, node, neig;
00230
00231 lnodes=NULL;
00232 vrstva=NULL;
00233
00234
00235 lnodes = new long [nn];
00236 memset (lnodes,0,sizeof(*lnodes)*nn);
00237 if ( lnodes == NULL)
00238 {
00239 printf(" nedostatek pameti\n");
00240 exit(-1);
00241 }
00242 vrstva = new long [nn+1];
00243 memset (vrstva,0,sizeof(*vrstva)*(nn+1));
00244
00245 if ( vrstva == NULL)
00246 {
00247 printf(" nedostatek pameti\n");
00248 exit(-1);
00249 }
00250
00251 for (i=1;i<nn;i++) lnodes[i]=-1;
00252 for (j=0;j<na;j++)
00253 {
00254
00255
00256 for( i=1;i<nn;i++)
00257 if (lnodes[i]==-1) break;
00258 cnt = 1;
00259 lnodes[i] = j;
00260 vrstva[1] = i;
00261 vrstva_start = 1;
00262 vrstva_konec = 1;
00263
00264 for(;cnt<local_n;)
00265 {
00266
00267
00268 for(i=vrstva_start; i<=vrstva_konec;i++)
00269 {
00270 node = vrstva[i];
00271
00272 for(k=0;k<nadjnodnod[node];k++)
00273 {
00274 neig = adjnodnod[node][k];
00275 if (lnodes[neig]==-1)
00276 {
00277 cnt++;
00278 lnodes[neig]=j;
00279 vrstva[cnt]=neig;
00280 if (cnt==local_n) break;
00281 }
00282 }
00283 if(cnt==local_n) break;
00284 }
00285 vrstva_start = vrstva_konec+1;
00286 vrstva_konec = cnt;
00287
00288 }
00289
00290 }
00291
00292 for(;;)
00293 {
00294
00295 cnt=1;
00296 for(node=0;node<nn;node++)
00297 {
00298 for(i=0;i<nadjnodnod[node];i++)
00299 {
00300 neig = adjnodnod[node][i];
00301 if (lnodes[neig]==-1)
00302 {
00303 lnodes[neig]=lnodes[i];
00304 cnt = 0;
00305 }
00306 }
00307 }
00308 if(cnt==1) break;
00309 }
00310
00311
00312
00313 lnntagr = new long [na];
00314 memset (lnntagr,0,sizeof(*lnntagr)*na);
00315
00316
00317 lntagr = new long * [na];
00318 lntagr[0] = new long [nn];
00319 for (i=0;i<nn;i++){
00320 lntagr[0][i]=0;
00321 }
00322
00323
00324
00325
00326 if ((lnntagr==NULL) ||(lntagr==NULL) ||(lntagr[0]==NULL))
00327 {
00328 printf("nedostatek pameti\n");
00329 exit(-1);
00330 }
00331
00332 for(i=0;i<na;i++) lnntagr[i]=0;
00333 for(i=0;i<nn;i++) lnntagr[lnodes[i]]++;
00334
00335
00336 cnt=0;
00337 for(i=0;i<na;i++)
00338 {
00339 for(j=0;j<nn;j++)
00340 if(lnodes[j]==i)
00341 {
00342 lntagr[0][cnt]=j;
00343 cnt++;
00344 }
00345 }
00346
00347 for(i=1;i<na;i++)
00348 lntagr[i]=lntagr[i-1]+lnntagr[i-1];
00349
00350 delete [] lnodes;
00351 delete [] vrstva;
00352
00353
00354
00355 }
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 void aggregator::metis_aggr (gtopology *gt)
00368 {
00369 long i,j;
00370
00371
00372 lnntagr = new long [na];
00373 for (i=0;i<na;i++){
00374 lnntagr[i]=gt->stop->nnsd[i];
00375 }
00376
00377
00378 lntagr = new long* [na];
00379 for (i=0;i<na;i++){
00380 lntagr[i] = new long [lnntagr[i]];
00381 }
00382 for (i=0;i<na;i++){
00383 for (j=0;j<lnntagr[i];j++){
00384 lntagr[i][j]=gt->stop->ltg[i][j];
00385 }
00386 }
00387
00388 }
00389
00390
00391
00392
00393
00394
00395
00396
00397 void aggregator::prepare_lnagr ()
00398 {
00399 long i,j,k,l,vrstva_zacatek,vrstva_konec,cnt;
00400 long node, neig;
00401 long *vrstva,*lnodes;
00402 vrstva=NULL;
00403 lnodes=NULL;
00404
00405
00406 vrstva = new long [nn];
00407 memset (vrstva,0,sizeof(*vrstva)*nn);
00408
00409 lnodes = new long [nn];
00410 memset (lnodes,0,sizeof(*lnodes)*nn);
00411
00412 lnnagr = new long [na];
00413 memset (lnnagr,0,sizeof(*lnnagr)*na);
00414
00415 lnagr = new long* [na];
00416
00417 for(i=0;i<na;i++)
00418 {
00419
00420 for(j=0; j<nn; j++) lnodes[j]=-1;
00421 for(j=0;j<lnntagr[i];j++)
00422 {
00423 vrstva[j]=lntagr[i][j];
00424 lnodes[vrstva[j]] = i;
00425 }
00426
00427 vrstva_zacatek = 0;
00428 vrstva_konec = lnntagr[i];
00429 for(k=0;k<degree_r;k++)
00430 {
00431 cnt = vrstva_konec;
00432 for(j=vrstva_zacatek;j<vrstva_konec;j++)
00433 {
00434
00435 node=vrstva[j];
00436 for(l=0;l<nadjnodnod[node];l++)
00437 {
00438 neig = adjnodnod[node][l];
00439 if (lnodes[neig]==-1)
00440 {
00441 lnodes[neig] = i;
00442 vrstva[cnt] = neig;
00443 cnt++;
00444 }
00445 }
00446 }
00447 vrstva_zacatek = vrstva_konec;
00448 vrstva_konec = cnt;
00449 }
00450
00451 lnnagr[i] = vrstva_konec;
00452
00453 lnagr[i] = new long [nn];
00454 for(l=0;l<vrstva_konec;l++)
00455 lnagr[i][l]=vrstva[l];
00456 }
00457
00458 delete [] vrstva;
00459 delete [] lnodes;
00460 }
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471 void aggregator::define_degree()
00472 {
00473 degree_r = 1;
00474 for(long i=0;i<degree_k;i++) degree_r *= 3;
00475 degree_r = (degree_r-1)/2;
00476 }
00477
00478
00479 void aggregator::assemble_smoothed_prol2 (long agrid,long cid,long *i,compvect *cvi,compvect *cvo,gmatrix *mtx)
00480 {
00481 compvect *cva;
00482
00483 cva = new compvect;
00484
00485
00486 mtx->cr->crxcv_cv (cvi,cvo);
00487
00488
00489
00490
00491
00492
00493
00494
00495 if (cvi->nz != cvo->nz){
00496 fprintf (stderr,"\n\n ruzne pocty nenulovych prvku ve vektorech cvi a cvo ve funkci assemble_smoothed_prol2");
00497 fprintf (stderr,"\n (file %s, line %d)\n",__FILE__,__LINE__);
00498 }
00499
00500 delete [] cva;
00501 cva=NULL;
00502 }
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515 void aggregator::assemble_smoothed_prol (gtopology *gt,gmatrix *mtx)
00516 {
00517 long i,j;
00518 double *v;
00519
00520 v = new double[n];
00521
00522 p = new double [na*nrbm*n];
00523
00524 switch (trbm){
00525 case 1:{
00526
00527
00528 for (i=0;i<na;i++){
00529 gener_rbm_1 (gt,i,v);
00530
00531
00532 mulS(degree_k,v,mtx);
00533
00534 copyv (v,p+i*n,n);
00535 }
00536 break;
00537 }
00538 case 2:{
00539
00540 for (i=0;i<na;i++){
00541 for (j=0;j<nrbm;j++){
00542
00543 gener_rbm_2 (gt,i,j,v);
00544
00545
00546 mulS(degree_k,v,mtx);
00547
00548 copyv (v,p+(i*nrbm+j)*n,n);
00549 }
00550 }
00551 break;
00552 }
00553 default:{
00554 print_err("unknown type of rigid body modes is required",__FILE__,__LINE__,__func__);
00555 }
00556 }
00557
00558 delete [] v;
00559 }
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569 void aggregator::mulS(long st,double *x, gmatrix *mtx)
00570 {
00571 double locro = specrad;
00572 double *wrk1;
00573 double *wrk2;
00574
00575 wrk1 = new double[n];
00576 wrk2 = new double[n];
00577
00578 for(long j=1;j<st;j++)
00579 {
00580 copyv(x,wrk1,n);
00581 copyv(x,wrk2,n);
00582 mulA(j-1,wrk2,mtx);
00583 for (long k=0; k<n;k++) x[k] = wrk1[k]-(4*wrk2[k])/(3*locro);
00584
00585 locro /= 9;
00586 }
00587 delete [] wrk1;
00588 delete [] wrk2;
00589 }
00590
00591
00592
00593
00594
00595
00596
00597 void aggregator::mulA(long st, double *x, gmatrix *mtx)
00598 {
00599 double *wrk1;
00600 double *wrk2;
00601 double locro = specrad;
00602
00603 wrk1 = new double[n];
00604 wrk2 = new double[n];
00605
00606 for(long k=1; k<st;k++) locro /= 9;
00607 copyv(x,wrk1,n);
00608 if ( st == 0 )
00609 mtx->gmxv(wrk1,x);
00610 else
00611 {
00612 mulA(st-1,x,mtx);
00613 copyv(x,wrk1,n);
00614 copyv(x,wrk2,n);
00615 mulA(st-1,wrk2,mtx);
00616 for (long k=0; k<n;k++) x[k] = wrk1[k]-(4*wrk2[k])/(3*locro);
00617
00618 copyv(x,wrk1,n);
00619 copyv(x,wrk2,n);
00620 mulA(st-1,wrk2,mtx);
00621 for (long k=0; k<n;k++) x[k] = wrk1[k]-(4*wrk2[k])/(3*locro);
00622
00623 }
00624 delete [] wrk1;
00625 delete [] wrk2;
00626 }
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638 void aggregator::gener_rbm_1 (gtopology *gt,long agid,double *v)
00639 {
00640 long i,cn,nid;
00641
00642 for (i=0;i<lnnagr[agid];i++){
00643 nid=lnagr[agid][i];
00644 cn=gt->give_dof (nid,0);
00645 if (cn>0)
00646 v[cn-1]=1.0;
00647 }
00648 }
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660 void aggregator::gener_rbm_2 (gtopology *gt,long agid,long cnid,double *v)
00661 {
00662 long i,cn,nid;
00663 double xp,yp,xk,yk;
00664
00665 if (cnid==2){
00666 nid=lnagr[agid][0];
00667 xp=gt->gnodes[nid].x;
00668 yp=gt->gnodes[nid].y;
00669 }
00670
00671 for (i=0;i<lnnagr[agid];i++){
00672 nid=lnagr[agid][i];
00673 if (cnid==0){
00674 cn=gt->give_dof (nid,cnid);
00675 if (cn>0)
00676 v[cn-1]=1.0;
00677 }
00678 if (cnid==1){
00679 cn=gt->give_dof (nid,cnid);
00680 if (cn>0)
00681 v[cn-1]=1.0;
00682 }
00683 if (cnid==2){
00684 xk=gt->gnodes[nid].x;
00685 yk=gt->gnodes[nid].y;
00686
00687 cn=gt->give_dof (nid,0);
00688 if (cn>0)
00689 v[cn-1]=yp-yk;
00690
00691 cn=gt->give_dof (nid,1);
00692 if (cn>0)
00693 v[cn-1]=xk-xp;
00694 }
00695 }
00696 }
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708 void aggregator::coarse_matrix (gmatrix *gm)
00709 {
00710 long i,j;
00711 double s,*am;
00712
00713
00714 cm = new densemat [1];
00715 cm->alloc (cms);
00716
00717
00718 am = new double [n*cms];
00719
00720 for (i=0;i<cms;i++){
00721 gm->gmxv (p+i*n,am+i*n);
00722 }
00723
00724 for (i=0;i<cms;i++){
00725 for (j=0;j<cms;j++){
00726 s=ss (p+i*n,am+j*n,n);
00727 cm->add_entry (s,i,j);
00728 }
00729 }
00730
00731 delete [] am;
00732 }
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745 void aggregator::assemble_aggr_unknowns (gtopology *gt)
00746 {
00747 long i,j,k,l,ann,nu,ndofn;
00748
00749 if (lnuaggr!=NULL){
00750 delete [] lnuaggr;
00751 }
00752 lnuaggr = new long [na];
00753
00754 maxnu=0;
00755 for (i=0;i<na;i++){
00756 nu=0;
00757 for (j=0;j<lnnagr[i];j++){
00758
00759 ann=lnagr[i][j];
00760
00761 ndofn = gt->give_ndofn (ann);
00762 for (k=0;k<ndofn;k++){
00763 l=gt->give_dof (ann,k);
00764 if (l>0)
00765 nu++;
00766 }
00767
00768 }
00769 lnuaggr[i]=nu;
00770 if (nu>maxnu)
00771 maxnu=nu;
00772 }
00773
00774 if (luaggr!=NULL){
00775 for (i=0;i<na;i++){
00776 delete [] luaggr[i];
00777 }
00778 delete [] luaggr;
00779 }
00780 luaggr = new long* [na];
00781 for (i=0;i<na;i++){
00782 luaggr[i] = new long [lnuaggr[i]];
00783 }
00784
00785
00786 for (i=0;i<na;i++){
00787 lnuaggr[i]=0;
00788 for (j=0;j<lnnagr[i];j++){
00789
00790 ann=lnagr[i][j];
00791
00792 ndofn = gt->give_ndofn (ann);
00793 for (k=0;k<ndofn;k++){
00794 l=gt->give_dof (ann,k);
00795 if (l>0){
00796 luaggr[i][lnuaggr[i]]=l;
00797 lnuaggr[i]++;
00798 }
00799 }
00800
00801 }
00802 }
00803
00804 }
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814 void aggregator::local_matrices (long bsize,gmatrix *gm)
00815 {
00816 long i,nu,pam,min,max;
00817
00818
00819
00820
00821
00822
00823 lmcr = new comprow [na];
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839 pam=0;
00840 min=LONG_MAX;
00841 max=0;
00842 for (i=0;i<na;i++){
00843 gm->cr->select_submatrix (luaggr[i],lnuaggr[i],&lmcr[i]);
00844
00845
00846
00847 nu=lmcr[i].n;
00848 if (lmcr[i].adr[nu]<min)
00849 min=lmcr[i].adr[nu];
00850 if (lmcr[i].adr[nu]>max)
00851 max=lmcr[i].adr[nu];
00852
00853 pam+=lmcr[i].adr[nu];
00854 fprintf (stdout,"\n number of nonzero matrix entries in compressed rows storage %ld",lmcr[i].adr[nu]);
00855 }
00856
00857 fprintf (stdout,"\n pocet prvku v maticich comp rows %10ld",pam);
00858 fprintf (stdout,"\n minimalni pocet prvku %10ld",min);
00859 fprintf (stdout,"\n maximalni pocet prvku %10ld",max);
00860
00861
00862 if (exinex==1){
00863
00864 lmsky = new skyline [na];
00865 for (i=0;i<na;i++){
00866 lmsky[i].assemble_from_cr (lmcr[i].n,lmcr[i].adr,lmcr[i].ci,lmcr[i].a);
00867 }
00868
00869 pam=0;
00870 min=LONG_MAX;
00871 max=0;
00872 for (i=0;i<na;i++){
00873 pam+=lmsky[i].negm;
00874 if (min>lmsky[i].negm)
00875 min=lmsky[i].negm;
00876 if (max<lmsky[i].negm)
00877 max=lmsky[i].negm;
00878 }
00879 fprintf (stdout,"\n pocet prvku v maticich skyline %10ld",pam);
00880 fprintf (stdout,"\n minimalni pocet prvku %10ld",min);
00881 fprintf (stdout,"\n maximalni pocet prvku %10ld",max);
00882 }
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933 min=LONG_MAX;
00934 max=0;
00935 for (i=0;i<na;i++){
00936 if (min>lnuaggr[i])
00937 min=lnuaggr[i];
00938 if (max<lnuaggr[i])
00939 max=lnuaggr[i];
00940 }
00941 fprintf (stdout,"\n minimalni pocet neznamych v agregatu %10ld",min);
00942 fprintf (stdout,"\n maximalni pocet neznamych v agregatu %10ld",max);
00943
00944 }
00945
00946 void aggregator::clean_memory ()
00947 {
00948 long i;
00949
00950 for (i=0;i<na;i++){
00951 delete [] lntagr[i];
00952 }
00953 delete [] lntagr;
00954
00955 delete [] lnntagr;
00956
00957
00958 for (i=0;i<na;i++){
00959 delete [] lnagr[i];
00960 }
00961 delete [] lnagr;
00962
00963 delete [] lnnagr;
00964 }
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976 void aggregator::prepare_boss (gtopology *gt,gmatrix *gm,FILE *out)
00977 {
00978 long i,j,k,l,v,min,bsize;
00979 time_t t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11;
00980 double *lhs,*rhs;
00981 lhs=NULL;
00982 rhs=NULL;
00983
00984
00985 na = gt->ns;
00986
00987
00988 nn = gt->nn;
00989
00990
00991 n = gm->n;
00992
00993
00994 switch (trbm){
00995 case 1:{
00996
00997 nrbm=1;
00998 break;
00999 }
01000 case 2:{
01001
01002 nrbm=3;
01003 break;
01004 }
01005 default:{
01006 print_err("unknown type of rigid body modes is required",__FILE__,__LINE__,__func__);
01007 }
01008 }
01009
01010
01011 cms=na*nrbm;
01012
01013
01014 specrad = gm->estim_spect_radius ();
01015
01016 t1 = time (NULL);
01017
01018 fprintf (stdout,"\n searching for adjacent nodes");
01019 gt->adjacnodes (out);
01020
01021
01022 nadjnodnod = gt->nadjnodnod;
01023
01024
01025 adjnodnod = gt->adjnodnod;
01026
01027 t2 = time (NULL);
01028
01029
01030
01031
01032
01033 if (impl==1){
01034 fprintf (stdout,"\n prepare_tlnagr");
01035
01036
01037
01038 prepare_tlnagr ();
01039 }
01040 if (impl==2){
01041
01042 metis_aggr (gt);
01043 }
01044
01045 t3 = time (NULL);
01046
01047
01048 define_degree ();
01049
01050 fprintf (stdout,"\n prepare_lnagr");
01051
01052
01053
01054 prepare_lnagr ();
01055
01056 t4 = time (NULL);
01057
01058 fprintf (stdout,"\n trideni");
01059
01060
01061
01062
01063 for (i=0;i<na;i++){
01064 for (j=0;j<lnnagr[i];j++){
01065 min=LONG_MAX;
01066 for (k=j;k<lnnagr[i];k++){
01067 if (min>lnagr[i][k]){
01068 min=lnagr[i][k];
01069 l=k;
01070 }
01071 }
01072 v=lnagr[i][j];
01073 lnagr[i][j]=lnagr[i][l];
01074 lnagr[i][l]=v;
01075 }
01076 }
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
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
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183 fprintf (stdout,"\n assemble_aggr_unknowns");
01184
01185 t5 = time (NULL);
01186
01187
01188
01189 assemble_aggr_unknowns (gt);
01190
01191
01192 t6 = time (NULL);
01193
01194
01195 bsize=gt->gnodes[0].ndofn;
01196
01197
01198 local_matrices (bsize,gm);
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216 t7 = time (NULL);
01217
01218 if (exinex==1){
01219
01220
01221
01222 for (i=0;i<na;i++){
01223
01224 lmsky[i].ldl_sky (lhs,rhs,1.0e-10,2);
01225 }
01226 }
01227
01228 t8 = time (NULL);
01229
01230
01231
01232
01233 assemble_smoothed_prol (gt,gm);
01234
01235 t9 = time (NULL);
01236
01237 coarse_matrix (gm);
01238
01239 t10 = time (NULL);
01240
01241 double zero=1.0e-15;
01242 cm->ll(lhs,rhs,zero,2);
01243
01244 t11 = time (NULL);
01245
01246 clean_memory ();
01247
01248
01249 fprintf (stdout,"\n\n Casy v BOSS");
01250 fprintf (stdout,"\n sestaveni sousednich uzlu %ld",t2-t1);
01251 fprintf (stdout,"\n sestaveni agregatu %ld",t3-t2);
01252 fprintf (stdout,"\n sestaveni agregatu %ld",t4-t3);
01253 fprintf (stdout,"\n trideni %ld",t5-t4);
01254 fprintf (stdout,"\n sestaveni neznamych %ld",t6-t5);
01255 fprintf (stdout,"\n lokalni matice %ld",t7-t6);
01256 fprintf (stdout,"\n rozklad lokalnich matic %ld",t8-t7);
01257 fprintf (stdout,"\n zhlazeny prolongator %ld",t9-t8);
01258 fprintf (stdout,"\n sestaveni coarse matice %ld",t10-t9);
01259 fprintf (stdout,"\n rozklad coarse matice %ld",t11-t10);
01260 fprintf (stdout,"\n");
01261
01262 }
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275 void aggregator::boss (gmatrix *gm,double *u,double *v)
01276 {
01277 long i;
01278 double *d,*dd,*pp,*z,*w,*ww;
01279
01280
01281 d = new double [n];
01282 dd = new double [maxnu];
01283 pp = new double [maxnu];
01284 z = new double [n];
01285 w = new double [cms];
01286 ww = new double [cms];
01287
01288
01289
01290
01291
01292 fillv (0.0,z,n);
01293
01294
01295
01296
01297 for (i=0;i<na;i++){
01298
01299 gm->gmxv (z,v);
01300
01301 subv (u,v,d,n);
01302
01303
01304 fillv (0.0,dd,maxnu);
01305
01306 globloc (d,dd,luaggr[i],lnuaggr[i]);
01307
01308
01309
01310
01311
01312 if (exinex==1){
01313
01314 lmsky[i].ldl_sky (pp,dd,1.0e-10,3);
01315 }
01316 if (exinex==2){
01317
01318 lmcr[i].cg (pp,dd,ssle->ni,ssle->err,ssle->ani,ssle->aerr,ssle->zero,0);
01319 }
01320
01321
01322 locglob (z,pp,luaggr[i],lnuaggr[i]);
01323 }
01324
01325
01326
01327
01328
01329
01330 gm->gmxv (z,v);
01331
01332 subv (u,v,d,n);
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358 for (i=na-1;i>-1;i--){
01359
01360 gm->gmxv (z,v);
01361
01362 subv (u,v,d,n);
01363
01364
01365 fillv (0.0,dd,maxnu);
01366
01367 globloc (d,dd,luaggr[i],lnuaggr[i]);
01368
01369
01370
01371
01372
01373 if (exinex==1){
01374
01375 lmsky[i].ldl_sky (pp,dd,1.0e-10,3);
01376 }
01377 if (exinex==2){
01378
01379 lmcr[i].cg (pp,dd,ssle->ni,ssle->err,ssle->ani,ssle->aerr,ssle->zero,0);
01380 }
01381
01382
01383 locglob (z,pp,luaggr[i],lnuaggr[i]);
01384 }
01385
01386
01387
01388
01389
01390 copyv (z,v,n);
01391
01392
01393 delete [] z;
01394 delete [] pp;
01395 delete [] dd;
01396 delete [] d;
01397 delete [] ww;
01398 delete [] w;
01399 }