00001 #include "seqfeti.h"
00002 #include "gmatrix.h"
00003 #include <limits.h>
00004
00005 seqfeti::seqfeti ()
00006 {
00007
00008
00009 fetiimpl=no_impl;
00010
00011
00012
00013 discont=1;
00014
00015 smst=skyline_matrix;
00016
00017
00018
00019 matassem=0;
00020
00021 vecassem=0;
00022
00023
00024
00025 ns=0;
00026
00027 ense=0;
00028
00029
00030
00031 ndofmax=0;
00032
00033 ndofcp=0;
00034
00035 thresh=0.0;
00036
00037 nicg=0;
00038
00039 anicg=0;
00040
00041 errcg=0.0;
00042
00043 aerrcg=0.0;
00044
00045 zero=1.0e-100;
00046
00047
00048 prec=noprecond;
00049
00050 gsize=0;
00051
00052 ndofmax=0;
00053
00054 ndofcp=0;
00055
00056
00057
00058 ndofmas=NULL;
00059
00060
00061 nsid=NULL;
00062
00063
00064
00065 cndom=NULL;
00066
00067
00068 ncdofd=NULL;
00069
00070
00071
00072 edofs=NULL;
00073
00074
00075 nrbmdom=NULL;
00076
00077
00078 rbmdom=NULL;
00079
00080
00081 rbmadr=NULL;
00082
00083
00084 se=NULL;
00085
00086
00087 g=NULL;
00088
00089
00090 invgg = NULL;
00091
00092
00093 ff = NULL;
00094
00095
00096 e = NULL;
00097
00098
00099 d = NULL;
00100
00101
00102 lambda = NULL;
00103
00104
00105
00106 b = NULL;
00107
00108
00109 h = NULL;
00110
00111
00112 ncbool=NULL;
00113
00114
00115 booldatar=NULL;
00116
00117
00118 booldatac=NULL;
00119
00120
00121 booldata=NULL;
00122
00123
00124 ccn=NULL;
00125
00126
00127
00128 smsky=NULL;
00129
00130 smdm=NULL;
00131
00132 }
00133
00134 seqfeti::~seqfeti ()
00135 {
00136 long i;
00137
00138
00139
00140
00141 for (i=0;i<ns;i++){
00142
00143 if (ccn!=NULL)
00144 delete [] ccn[i];
00145
00146 if (booldata!=NULL)
00147 delete [] booldata[i];
00148
00149 if (booldatac!=NULL)
00150 delete [] booldatac[i];
00151
00152 if (booldatar!=NULL)
00153 delete [] booldatar[i];
00154
00155 if (se!=NULL)
00156 delete [] se[i];
00157
00158 if (rbmdom!=NULL)
00159 delete [] rbmdom[i];
00160
00161 if (edofs!=NULL)
00162 delete [] edofs[i];
00163
00164 if (cndom!=NULL)
00165 delete [] cndom[i];
00166
00167 if (ff!=NULL)
00168 delete [] ff[i];
00169
00170 if (d!=NULL)
00171 delete [] d[i];
00172 }
00173
00174
00175 delete [] ccn;
00176
00177 delete [] booldata;
00178
00179 delete [] booldatac;
00180
00181 delete [] booldatar;
00182
00183 delete [] se;
00184
00185 delete [] rbmdom;
00186
00187 delete [] edofs;
00188
00189 delete [] cndom;
00190
00191 delete [] ff;
00192
00193 delete [] d;
00194
00195
00196
00197
00198
00199
00200 delete [] ncbool;
00201
00202 delete [] rbmadr;
00203
00204 delete [] nrbmdom;
00205
00206 delete [] ncdofd;
00207
00208 delete [] nsid;
00209
00210 delete [] ndofmas;
00211
00212
00213 delete [] g;
00214
00215 delete [] invgg;
00216
00217 delete [] e;
00218
00219 delete [] lambda;
00220
00221 delete [] b;
00222
00223 delete [] h;
00224
00225
00226 delete [] smsky;
00227
00228 delete [] smdm;
00229
00230 }
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240 void seqfeti::read (gtopology *top,XFILE *in)
00241 {
00242
00243
00244
00245
00246 xfscanf (in,"%m",&fetiimplem_kwdset,&fetiimpl);
00247
00248 if (fetiimpl!=boolean_matrices && fetiimpl!=nonredundant && fetiimpl!=redundant){
00249 print_err("unknown type of sequential FETI implementation",__FILE__,__LINE__,__func__);
00250 }
00251
00252
00253
00254
00255
00256
00257 xfscanf (in,"%ld %ld %le %ld %le",&ns,&ense,&thresh,&nicg,&errcg);
00258
00259 xfscanf (in,"%d",(int*)&prec);
00260
00261
00262
00263
00264
00265 if (fetiimpl==nonredundant || fetiimpl==redundant){
00266
00267 top->rst=1;
00268 }
00269 if (fetiimpl==boolean_matrices){
00270
00271
00272 top->rst=2;
00273 }
00274 if (prec==feti_dirichlet){
00275
00276
00277 top->cngen=2;
00278 }
00279
00280
00281 top->cngen=0;
00282
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292 void seqfeti::print (FILE *out)
00293 {
00294
00295 fprintf (out,"%d",fetiimpl);
00296
00297
00298
00299
00300
00301
00302 fprintf (out,"%ld %ld %le %ld %le\n",ns,ense,thresh,nicg,errcg);
00303
00304 fprintf (out,"%d\n",prec);
00305 }
00306
00307
00308
00309
00310
00311
00312
00313
00314 void seqfeti::read_booldata (XFILE *in)
00315 {
00316 long i,j;
00317
00318
00319 ncbool = new long [ns];
00320 booldatar = new long* [ns];
00321 booldatac = new long* [ns];
00322 booldata = new double* [ns];
00323
00324 for (i=0;i<ns;i++){
00325
00326
00327
00328 xfscanf (in,"%ld",ncbool+i);
00329
00330
00331 booldatar[i] = new long [ncbool[i]];
00332 booldatac[i] = new long [ncbool[i]];
00333 booldata[i] = new double [ncbool[i]];
00334
00335 for (j=0;j<ncbool[i];j++){
00336 xfscanf (in,"%ld %ld %lf",&booldatar[i][j],&booldatac[i][j],&booldata[i][j]);
00337 }
00338 }
00339
00340 }
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 void seqfeti::initiate (seqselnodes *selnodfeti,gtopology *top,FILE *out)
00356 {
00357 long i,j,k,ndofn,dof,nid;
00358 long *red;
00359 red = new long [ns+1];
00360
00361
00362 ndofcp = selnodfeti->tndofsn;
00363
00364
00365 ncdofd = new long [ns];
00366 for (i=0;i<ns;i++){
00367 ncdofd[i] = selnodfeti->snndofmas[i];
00368 }
00369
00370
00371 edofs = new long* [ns];
00372 for (i=0;i<ns;i++){
00373 edofs[i] = new long [ncdofd[i]];
00374 }
00375
00376 for (i=0;i<ns;i++){
00377 for (j=0;j<ncdofd[i];j++){
00378 edofs[i][j] = selnodfeti->lndofmas[i][j];
00379 }
00380 }
00381
00382
00383
00384
00385
00386
00387 red[0]=0;
00388 nid=0;
00389
00390 for (i=0;i<ns;i++){
00391 red[i+1]=0;
00392
00393 for (j=0;j<top->stop->nnsd[i];j++){
00394 ndofn = top->give_ndofn (nid);
00395 for (k=0;k<ndofn;k++){
00396 dof = top->give_dof (nid,k);
00397 if (dof>red[i+1]){
00398 red[i+1]=dof;
00399 }
00400 }
00401 nid++;
00402 }
00403 }
00404
00405
00406
00407
00408 for (i=0;i<ns;i++){
00409 for (j=0;j<ncdofd[i];j++){
00410 if (edofs[i][j]>0)
00411 edofs[i][j] -= red[i];
00412 else
00413 edofs[i][j] += red[i];
00414 }
00415 }
00416
00417
00418 fprintf (out,"\n\n\n check of the array edofs \n");
00419 for (i=0;i<ns;i++){
00420 fprintf (out,"\n domain %ld",i);
00421 for (j=0;j<ncdofd[i];j++){
00422 fprintf (out,"\n edofs %6ld %6ld",j,edofs[i][j]);
00423 }
00424 }
00425
00426
00427
00428
00429
00430
00431 wscalmat = new double [ndofcp];
00432 for (i = 0; i < ndofcp; i++){
00433
00434 }
00435
00436
00437
00438
00439
00440 ccn = new long* [ns];
00441 for (i=0;i<ns;i++){
00442 ccn[i] = new long [ncdofd[i]];
00443 }
00444
00445 for (i=0;i<ns;i++){
00446 for (j=0;j<ncdofd[i];j++){
00447 ccn[i][j] = selnodfeti->cndofmas[i][j];
00448 }
00449 }
00450
00451 delete [] red;
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465 }
00466
00467
00468
00469
00470
00471
00472
00473
00474 void seqfeti::det_ndofmax ()
00475 {
00476 long i,j;
00477
00478 if (fetiimpl==nonredundant || fetiimpl==redundant){
00479 ndofmax=0;
00480 for (i=0;i<ns;i++){
00481 if (ndofmax<ndofmas[i])
00482 ndofmax=ndofmas[i];
00483 }
00484 }
00485
00486 if (fetiimpl==boolean_matrices){
00487 ndofmax=0;
00488 ndofcp=0;
00489 for (i=0;i<ns;i++){
00490 for (j=0;j<ncbool[i];j++){
00491 if (ndofmax<booldatac[i][j])
00492 ndofmax=booldatac[i][j];
00493 if (ndofcp<booldatar[i][j])
00494 ndofcp=booldatar[i][j];
00495 }
00496 }
00497 }
00498 }
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513 void seqfeti::assemble_subdom_unknowns (gtopology *top,FILE *out)
00514 {
00515 long i,j,k,l,m,ndofn;
00516
00517
00518 nsid = new long [top->nn];
00519 l=0;
00520 for (i=0;i<ns;i++){
00521 for (j=0;j<top->stop->nnsd[i];j++){
00522 nsid[l]=i;
00523 l++;
00524 }
00525 }
00526
00527
00528
00529 ndofmas = new long [ns];
00530 for (i=0;i<ns;i++){
00531 ndofmas[i]=0;
00532 }
00533
00534 for (i=0;i<top->nn;i++){
00535
00536 ndofn=top->give_ndofn (i);
00537
00538 l=nsid[i];
00539 for (j=0;j<ndofn;j++){
00540
00541 k=top->give_dof (i,j);
00542 if (k>0)
00543 ndofmas[l]++;
00544 }
00545 }
00546
00547
00548 cndom = new long* [ns];
00549 for (i=0;i<ns;i++){
00550 cndom[i] = new long [ndofmas[i]];
00551 ndofmas[i]=0;
00552 }
00553
00554 for (i=0;i<top->nn;i++){
00555
00556 ndofn=top->give_ndofn (i);
00557
00558 l=nsid[i];
00559 for (j=0;j<ndofn;j++){
00560
00561 k=top->give_dof (i,j);
00562 if (k>0){
00563 cndom[l][ndofmas[l]]=k;
00564 ndofmas[l]++;
00565 }
00566 }
00567 }
00568
00569
00570
00571
00572
00573 if (prec == feti_dirichlet){
00574 for (i=0;i<ns;i++){
00575 for (j=0;j<ndofmas[i];j++){
00576 l=cndom[i][j];
00577 m=j;
00578 for (k=j;k<ndofmas[i];k++){
00579 if (l>cndom[i][k]){
00580 l=cndom[i][k];
00581 m=k;
00582 }
00583 }
00584 cndom[i][m]=cndom[i][j];
00585 cndom[i][j]=l;
00586 }
00587 }
00588 }
00589
00590
00591 fprintf (out,"\n\n kontrola pole ndofmas");
00592 for (i=0;i<ns;i++){
00593 fprintf (out,"\n ndofmas %6ld %6ld",i,ndofmas[i]);
00594 }
00595 fprintf (out,"\n\n kontrola pole cndom");
00596 for (i=0;i<ns;i++){
00597 fprintf (out,"\n domain %ld",i);
00598 for (j=0;j<ndofmas[i];j++){
00599 fprintf (out,"\n dom %6ld dof %6ld %6ld",i,j,cndom[i][j]);
00600 }
00601 }
00602
00603 }
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614 void seqfeti::subdomain_matrices (gmatrix *gm,FILE *out)
00615 {
00616 long i,j,k;
00617
00618 if (smst==dense_matrix){
00619 smdm = new densemat [ns];
00620 for (i=0;i<ns;i++){
00621
00622 }
00623 }
00624 if (smst==skyline_matrix){
00625 smsky = new skyline [ns];
00626 for (i=0;i<ns;i++){
00627 smsky[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofmas[i],cndom[i]);
00628 }
00629 }
00630
00631
00632 if (fetiimpl==nonredundant || fetiimpl==redundant){
00633 long *aux;
00634 ndofprec = new long [ns];
00635 cnprec = new long* [ns];
00636 cpreccn = new long* [ns];
00637
00638 for (i=0;i<ns;i++){
00639 aux = new long [ndofmas[i]];
00640 for (j=0;j<ndofmas[i];j++){
00641 aux[j]=0;
00642 }
00643
00644 for (j=0;j<ncdofd[i];j++){
00645 k=edofs[i][j];
00646 if (k<0)
00647 k=0-k-1;
00648 else
00649 k--;
00650 if (k>-1)
00651 aux[k]++;
00652 }
00653
00654 ndofprec[i]=0;
00655 for (j=0;j<ndofmas[i];j++){
00656 if (aux[j]>0)
00657 ndofprec[i]++;
00658 }
00659
00660 fprintf (out,"\n\n aux domena %ld\n",i);
00661 for (j=0;j<ndofmas[i];j++){
00662 fprintf (out," %ld",aux[j]);
00663 }
00664
00665
00666 cnprec[i] = new long [ndofprec[i]];
00667 cpreccn[i] = new long [ndofprec[i]];
00668 k=0;
00669 for (j=0;j<ndofmas[i];j++){
00670 if (aux[j]>0){
00671 cnprec[i][k]=cndom[i][j];
00672 cpreccn[i][k]=j;
00673 k++;
00674 }
00675 }
00676
00677 delete [] aux;
00678 }
00679
00680 fprintf (out,"\n\n\n\n kontrola ve funkci subdomain_matrices (file %s, line %d)\n\n",__FILE__,__LINE__);
00681 for (i=0;i<ns;i++){
00682 fprintf (out,"\n domena %5ld ndofprec %ld\n",i,ndofprec[i]);
00683 for (j=0;j<ndofprec[i];j++){
00684 fprintf (out," %ld",cnprec[i][j]);
00685 }
00686 fprintf (out,"\n");
00687 for (j=0;j<ndofprec[i];j++){
00688 fprintf (out," %ld",cpreccn[i][j]);
00689 }
00690 }
00691 fprintf (out,"\n\n\n\n\n");
00692 }
00693
00694
00695 if (fetiimpl==nonredundant || fetiimpl==redundant){
00696 if (prec==feti_lumped){
00697 smscr = new symcomprow [ns];
00698 for (i=0;i<ns;i++){
00699
00700 gm->scr->select_submatrix (&smscr[i],ndofprec[i],cnprec[i]);
00701 }
00702 }
00703
00704 if (prec==feti_dirichlet){
00705 psmdm = new densemat [ns];
00706
00707 double *x,*y,*condvect;
00708 for (i=0;i<ns;i++){
00709 x = new double [ndofmas[i]];
00710 y = new double [ndofmas[i]];
00711 psmdm[i].a = new double [ndofprec[i]*ndofprec[i]];
00712 psmdm[i].negm=ndofprec[i]*ndofprec[i];
00713 psmdm[i].n=ndofprec[i];
00714 condvect=new double [ndofprec[i]];
00715
00716 smsky[i].ldlkon_sky (psmdm[i].a,condvect,x,y,ndofprec[i],1);
00717 smsky[i].~skyline ();
00718
00719 delete [] x;
00720 delete [] y;
00721 delete [] condvect;
00722 }
00723
00724 smsky = new skyline [ns];
00725 for (i=0;i<ns;i++){
00726 smsky[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofmas[i],cndom[i]);
00727 }
00728
00729 }
00730 }
00731
00732 if (fetiimpl==boolean_matrices){
00733 if (prec==feti_lumped){
00734 smskyprec = new skyline [ns];
00735 for (i=0;i<ns;i++){
00736 smskyprec[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofmas[i],cndom[i]);
00737 }
00738 }
00739 if (prec==feti_dirichlet){
00740
00741 ndofprec = new long[ns];
00742 for (i=0;i<ns;i++){
00743 ndofprec[i] = 1;
00744 for(j = 1; j < ncbool[i]; j++){
00745 if(booldatac[i][j] != booldatac[i][j-1]){
00746 ndofprec[i]++;
00747 }
00748 }
00749 }
00750
00751 psmdm = new densemat [ns];
00752
00753 double *x,*y,*condvect;
00754 for (i=0;i<ns;i++){
00755 x = new double [ndofmas[i]];
00756 y = new double [ndofmas[i]];
00757 psmdm[i].a = new double [ndofprec[i]*ndofprec[i]];
00758 psmdm[i].negm=ndofprec[i]*ndofprec[i];
00759 psmdm[i].n=ndofprec[i];
00760 condvect=new double [ndofprec[i]];
00761
00762 smsky[i].ldlkon_sky (psmdm[i].a,condvect,x,y,ndofprec[i],1);
00763 smsky[i].~skyline ();
00764
00765 delete [] x;
00766 delete [] y;
00767 delete [] condvect;
00768 }
00769
00770 smsky = new skyline [ns];
00771 for (i=0;i<ns;i++){
00772 smsky[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofmas[i],cndom[i]);
00773 }
00774 }
00775 }
00776
00777 }
00778
00779
00780
00781
00782
00783
00784 void seqfeti::kernel (FILE *out)
00785 {
00786 long i,j;
00787
00788
00789 if (nrbmdom!=NULL){
00790 delete [] nrbmdom;
00791 }
00792 nrbmdom = new long [ns];
00793
00794
00795 if (se!=NULL){
00796 for (i=0;i<ns;i++){
00797 delete [] se[i];
00798 }
00799 delete [] se;
00800 }
00801 se = new long* [ns];
00802 for (i=0;i<ns;i++){
00803 se[i] = new long [ense];
00804 for (j=0;j<ense;j++){
00805 se[i][j]=-1;
00806 }
00807 }
00808
00809
00810 if (rbmdom!=NULL){
00811 for (i=0;i<ns;i++){
00812 delete [] rbmdom[i];
00813 }
00814 delete [] rbmdom;
00815 }
00816 rbmdom = new double* [ns];
00817 for (i=0;i<ns;i++){
00818 rbmdom[i] = new double [ense*ndofmas[i]];
00819 for (j=0;j<ense*ndofmas[i];j++){
00820 rbmdom[i][j]=0.0;
00821 }
00822 }
00823
00824
00825 if (smst==dense_matrix){
00826
00827 for (i=0;i<ns;i++){
00828 smdm[i].ker (rbmdom[i],nrbmdom[i],se[i],ense,thresh);
00829 }
00830 }
00831 if (smst==skyline_matrix){
00832
00833 for (i=0;i<ns;i++){
00834 smsky[i].ker (rbmdom[i],nrbmdom[i],se[i],ense,thresh,3);
00835
00836 fprintf (stdout,"\n the number of rigid body motions %ld",nrbmdom[i]);
00837
00838
00839 }
00840 }
00841
00842 }
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865 void seqfeti::boolean_matrix (long sdid,double *a)
00866 {
00867 long i,ri,ci;
00868
00869 for (i=0;i<ncbool[sdid];i++){
00870 ri=booldatar[sdid][i]-1;
00871 ci=booldatac[sdid][i]-1;
00872 a[ri*ndofmas[sdid]+ci]=booldata[sdid][i];
00873 }
00874 }
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887 void seqfeti::local_coarse (long nd,double *lv,double *cv)
00888 {
00889 if (fetiimpl==nonredundant || fetiimpl==redundant){
00890
00891 long i,j,k;
00892
00893 for (i=0;i<ncdofd[nd];i++){
00894 j=edofs[nd][i];
00895 k=ccn[nd][i]-1;
00896 if (j<0){
00897 j=0-j-1;
00898 cv[k]-=lv[j];
00899 }
00900 else{
00901 j=j-1;
00902 cv[k]+=lv[j];
00903 }
00904 }
00905 }
00906
00907 if (fetiimpl==boolean_matrices){
00908
00909 double *a,*av;
00910
00911
00912 a = new double [ndofcp*ndofmas[nd]];
00913 nullv (a,ndofcp*ndofmas[nd]);
00914
00915
00916 av = new double [ndofcp];
00917 nullv (av,ndofcp);
00918
00919
00920 boolean_matrix (nd,a);
00921
00922 mxv (a,lv,av,ndofcp,ndofmas[nd]);
00923
00924 addv (cv,av,ndofcp);
00925
00926 delete [] a;
00927 delete [] av;
00928 }
00929
00930 }
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942 void seqfeti::coarse_local (long nd,double *lv,double *cv)
00943 {
00944 if (fetiimpl==nonredundant || fetiimpl==redundant){
00945 long i,j,k;
00946
00947 for (i=0;i<ncdofd[nd];i++){
00948 k=ccn[nd][i]-1;
00949 j=edofs[nd][i];
00950 if (j<0){
00951 j=0-j-1;
00952 lv[j]-=cv[k];
00953 }
00954 else{
00955 j=j-1;
00956 lv[j]+=cv[k];
00957 }
00958 }
00959 }
00960
00961 if (fetiimpl==boolean_matrices){
00962
00963 double *a,*av;
00964
00965
00966 a = new double [ndofcp*ndofmas[nd]];
00967 nullv (a,ndofcp*ndofmas[nd]);
00968
00969
00970 av = new double [ndofmas[nd]];
00971 nullv (av,ndofmas[nd]);
00972
00973
00974 boolean_matrix (nd,a);
00975
00976 mtxv (a,cv,av,ndofcp,ndofmas[nd]);
00977
00978 addv (lv,av,ndofmas[nd]);
00979
00980 delete [] a;
00981 delete [] av;
00982 }
00983
00984 }
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999 void seqfeti::g_matrixsize (FILE *out)
01000 {
01001 long i;
01002
01003
01004 rbmadr = new long [ns+1];
01005 for (i=0;i<ns+1;i++){
01006 rbmadr[i]=0;
01007 }
01008
01009
01010 rbmadr[0]=0;
01011 for (i=0;i<ns;i++){
01012 rbmadr[i+1]=rbmadr[i]+nrbmdom[i];
01013 }
01014
01015
01016 gsize=rbmadr[ns];
01017
01018 }
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028 void seqfeti::g_matrix (FILE *out)
01029 {
01030 long i,k,l;
01031 double *av;
01032
01033
01034 if (g==NULL)
01035 g = new double [ndofcp*gsize];
01036
01037
01038 av = new double [ndofcp];
01039
01040 for (i=0;i<ns;i++){
01041 for (k=0;k<nrbmdom[i];k++){
01042 nullv (av,ndofcp);
01043 local_coarse (i,rbmdom[i]+k*ndofmas[i],av);
01044 for (l=0;l<ndofcp;l++){
01045 g[l*gsize+rbmadr[i]+k]=0.0-av[l];
01046 }
01047 }
01048 }
01049
01050 delete [] av;
01051 }
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061 void seqfeti::evector (FILE *out)
01062 {
01063 long i,j,k;
01064
01065 if (e==NULL)
01066 e = new double [gsize];
01067
01068 k=0;
01069 for (i=0;i<ns;i++){
01070 for (j=0;j<nrbmdom[i];j++){
01071 e[k]=0.0-ss(&rbmdom[i][j*ndofmas[i]],ff[i],ndofmas[i]);
01072 k++;
01073 }
01074 }
01075
01076 fprintf (out,"\n\n\n kontrola vektoru e \n");
01077 for (j=0;j<gsize;j++){
01078 fprintf (out,"\n %4ld %f",j,e[j]);
01079 }
01080
01081 }
01082
01083
01084
01085
01086
01087
01088 void seqfeti::inverse_matrix_GG (FILE *out)
01089 {
01090 long i,j;
01091 double *ee,*gg;
01092
01093 if (invgg == NULL){
01094 invgg = new double [gsize*gsize];
01095 }
01096 gg = new double [gsize*gsize];
01097 ee = new double [gsize*gsize];
01098
01099 for (i=0;i<gsize;i++){
01100 for (j=0;j<gsize;j++){
01101 ee[i*gsize+j]=0.0;
01102 }
01103 ee[i*gsize+i]=1.0;
01104 }
01105
01106 mtxm (g,g,gg,ndofcp,gsize,gsize);
01107
01108 fprintf (out,"\n\n\n kontrola matice G^T G \n");
01109 for (i=0;i<gsize;i++){
01110 fprintf (out,"\n %4ld ",i);
01111 for (j=0;j<gsize;j++){
01112 fprintf (out," %lf",gg[i*gsize+j]);
01113 }
01114 }
01115
01116 gemp (gg,invgg,ee,gsize,gsize,zero,1);
01117
01118 fprintf (out,"\n\n\n kontrola matice (G^T G)^{-1} \n");
01119 for (i=0;i<gsize;i++){
01120 fprintf (out,"\n %4ld ",i);
01121 for (j=0;j<gsize;j++){
01122 fprintf (out," %lf",invgg[i*gsize+j]);
01123 }
01124 }
01125
01126 delete [] gg;
01127 delete [] ee;
01128 }
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139 void seqfeti::feti_projection (double *v)
01140 {
01141 long i;
01142 double *p,*q,*w;
01143
01144 p = new double [gsize];
01145 q = new double [gsize];
01146 w = new double [ndofcp];
01147
01148
01149 mtxv (g,v,p,ndofcp,gsize);
01150
01151
01152 mxv (invgg,p,q,gsize,gsize);
01153
01154
01155 mxv (g,q,w,ndofcp,gsize);
01156
01157
01158 for (i=0;i<ndofcp;i++){
01159 v[i]-=w[i];
01160 }
01161
01162 delete [] w; delete [] q; delete [] p;
01163 }
01164
01165
01166
01167
01168
01169
01170
01171
01172 void seqfeti::define_h (FILE *out)
01173 {
01174 long i;
01175
01176 if (h==NULL){
01177 h = new double [ndofcp];
01178 }
01179
01180 for (i=0;i<ndofcp;i+=2){
01181 h[i]=1.0e-5;
01182 h[i+1]=0.0;
01183 }
01184 }
01185
01186
01187
01188
01189
01190
01191
01192
01193 void seqfeti::define_b (FILE *out)
01194 {
01195 long i;
01196
01197 if (b==NULL){
01198 b = new double [ndofcp];
01199 }
01200
01201 for (i=0;i<ndofcp;i++){
01202 b[i]=0.0;
01203 }
01204 }
01205
01206
01207 void seqfeti::scaling(double *invect,double *outvect,long n,FILE *out)
01208 {
01209 double *aux;
01210 long i;
01211 aux = new double[n];
01212
01213 for(i = 0; i < n; i++){
01214 aux[i] = invect[i]/wscalmat[i];
01215 }
01216
01217
01218
01219
01220
01221
01222
01223
01224 for(i = 0; i < n; i++){
01225 outvect[i] = aux[i];
01226 }
01227
01228 delete []aux;
01229 }
01230
01231 void seqfeti::lscaling(double *invect,double *outvect,long ndom,FILE *out)
01232 {
01233 double *aux;
01234 long i;
01235
01236 aux = new double[nlwscalmat[ndom]];
01237
01238 for(i = 0; i < nlwscalmat[ndom]; i++){
01239 aux[i] = invect[i]/lwscalmat[ndom][i];
01240 }
01241
01242
01243 fprintf (out,"\n\n\n kontrola skalovani \n");
01244 for(i = 0; i < nlwscalmat[ndom]; i++){
01245 fprintf (out,"%ld pred %le po %le wscalmat %lf ccn %ld\n",i+1,invect[i],aux[i],lwscalmat[ndom][i],ccn[ndom][i]);
01246 }
01247
01248
01249 for(i = 0; i < nlwscalmat[ndom]; i++){
01250 outvect[i] = aux[i];
01251 }
01252
01253 delete []aux;
01254 }
01255
01256
01257
01258
01259
01260
01261 void seqfeti::lumpedprec (long nd,double *dd,double *pp)
01262 {
01263 long i;
01264 double *d,*p;
01265
01266 if (fetiimpl == nonredundant || fetiimpl==redundant){
01267 d = new double [ndofprec[nd]];
01268 p = new double [ndofprec[nd]];
01269
01270 for (i=0;i<ndofprec[nd];i++){
01271 d[i]=0.0;
01272 p[i]=0.0;
01273 }
01274
01275 for (i=0;i<ndofprec[nd];i++){
01276 d[i]=dd[cpreccn[nd][i]];
01277 }
01278
01279 smscr[nd].mxv_scr (d,p);
01280
01281 for (i=0;i<ndofprec[nd];i++){
01282 pp[cpreccn[nd][i]]=p[i];
01283 }
01284 delete [] p;
01285 delete [] d;
01286 }
01287 if(fetiimpl == boolean_matrices){
01288
01289
01290
01291
01292
01293 smskyprec[nd].mxv_sky (dd,pp);
01294
01295
01296 }
01297 }
01298
01299
01300
01301
01302 void seqfeti::dirichletprec (long nd,double *dd,double *pp)
01303 {
01304 long i,j;
01305 double *d,*p;
01306
01307 if (fetiimpl == nonredundant || fetiimpl==redundant){
01308 d = new double [ndofprec[nd]];
01309 p = new double [ndofprec[nd]];
01310
01311 for (i=0;i<ndofprec[nd];i++){
01312 d[i]=0.0;
01313 p[i]=0.0;
01314 }
01315
01316 for (i=0;i<ndofprec[nd];i++){
01317 d[i]=dd[cpreccn[nd][i]];
01318 }
01319
01320 psmdm[nd].mxv_dm (d,p);
01321
01322
01323 for (i=0;i<ndofprec[nd];i++){
01324 pp[cpreccn[nd][i]]=p[i];
01325 }
01326
01327 delete [] p;
01328 delete [] d;
01329 }
01330 if (fetiimpl == boolean_matrices){
01331 d = new double [ndofprec[nd]];
01332 p = new double [ndofprec[nd]];
01333 j = 0;
01334 for (i=0;i<ndofmas[nd];i++){
01335 if(dd[i] != 0){
01336 d[j]=dd[i];
01337 j++;
01338 }
01339 }
01340 psmdm[nd].mxv_dm (d,p);
01341
01342 j = 0;
01343 for (i=0;i<ndofmas[nd];i++){
01344 if(dd[i] != 0){
01345 pp[i]=p[j];
01346 j++;
01347 }
01348 else{
01349 pp[i] = 0;
01350 }
01351 }
01352 }
01353 }
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367 void seqfeti::mpcg (FILE *out)
01368 {
01369 long i,j;
01370 double nom,denom,alpha,beta;
01371 double *s,*sd,*r,*p,*pp,*invgge;
01372
01373 if (lambda==NULL)
01374 lambda = new double [ndofcp];
01375
01376 if (ndofmax==0){
01377 print_err("the attribute ndofmax is not determined",__FILE__,__LINE__,__func__);
01378 }
01379
01380
01381 sd = new double [ndofmax];
01382 pp = new double [ndofmax];
01383
01384
01385 s = new double [ndofcp];
01386
01387 r = new double [ndofcp];
01388
01389 p = new double [ndofcp];
01390
01391
01392
01393 invgge = new double [gsize];
01394
01395 mxv (invgg,e,invgge,gsize,gsize);
01396
01397 mxv (g,invgge,lambda,ndofcp,gsize);
01398 delete [] invgge;
01399
01400 nullv (r,ndofcp);
01401 for (j=0;j<ns;j++){
01402 nullv (sd,ndofmas[j]);
01403 coarse_local (j,sd,lambda);
01404
01405
01406
01407 subv (sd,ff[j],ndofmas[j]);
01408
01409
01410 if (smst==dense_matrix)
01411 break;
01412 if (smst==skyline_matrix)
01413 smsky[j].ldl_feti_sky (pp,sd,nrbmdom[j],se[j],zero);
01414
01415
01416
01417
01418 local_coarse (j,pp,r);
01419 }
01420
01421 cmulv(-1.0,r,ndofcp);
01422
01423
01424 if (discont==1){
01425
01426
01427
01428 for (i=0;i<ndofcp;i++){
01429 r[i]-=b[i]+h[i]*lambda[i];
01430 }
01431 }
01432
01433
01434 feti_projection (r);
01435
01436
01437 for (i=0;i<ndofcp;i++){
01438 s[i]=r[i];
01439 }
01440
01441
01442 nom = ss (r,r,ndofcp);
01443
01444
01445
01446
01447
01448 for (i=0;i<nicg;i++){
01449
01450
01451
01452
01453 nullv (p,ndofcp);
01454 for (j=0;j<ns;j++){
01455 nullv (sd,ndofmax);
01456 coarse_local (j,sd,s);
01457
01458
01459 if (smst==dense_matrix)
01460 break;
01461 if (smst==skyline_matrix)
01462 smsky[j].ldl_feti_sky (pp,sd,nrbmdom[j],se[j],zero);
01463
01464 local_coarse (j,pp,p);
01465 }
01466
01467 if (discont==1){
01468
01469
01470
01471 for (j=0;j<ndofcp;j++){
01472
01473 p[j]+=h[j]*s[j];
01474 }
01475 }
01476
01477
01478 denom = ss (s,p,ndofcp);
01479
01480 fprintf (out,"\n MPCG in FETI, iterace %ld, kontrola citatele a jmenovatele alpha %e / %e",i,nom,denom);
01481
01482 if (fabs(denom)<zero){
01483 print_err("zero denominator of alpha in the modified conjugate gradient method",__FILE__,__LINE__,__func__);
01484 break;
01485 }
01486
01487
01488
01489
01490 alpha = nom/denom;
01491
01492
01493
01494
01495 for (j=0;j<ndofcp;j++){
01496 r[j]-=alpha*p[j];
01497 lambda[j]+=alpha*s[j];
01498 }
01499
01500 feti_projection (r);
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581 denom = nom;
01582 if (fabs(denom)<zero){
01583 print_err("zero denominator of beta in the modified conjugate gradient method",__FILE__,__LINE__,__func__);
01584 break;
01585 }
01586
01587 nom=ss(r,r,ndofcp);
01588
01589 fprintf (stdout,"\n iteration %4ld norm r %e",i,sqrt(nom));
01590 fprintf (out,"\n iteration %4ld norm r %e",i,sqrt(nom));
01591
01592 fprintf (out,"\n kontrola citatele a jmenovatele pred betou %e / %e",nom,denom);
01593
01594 if (sqrt(nom)<errcg){
01595 break;
01596 }
01597
01598
01599
01600 beta = nom/denom;
01601
01602
01603 for (j=0;j<ndofcp;j++){
01604 s[j]=beta*s[j]+r[j];
01605 }
01606
01607 }
01608
01609
01610 anicg=i;
01611
01612 aerrcg=nom;
01613
01614 delete [] sd;
01615 delete [] pp;
01616 delete [] s;
01617 delete [] r;
01618 delete [] p;
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
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871 void seqfeti::nodalunknowns (double *lhs, FILE *out)
01872 {
01873 long i,j;
01874
01875 if (d == NULL){
01876 d = new double* [ns];
01877 for (i=0;i<ns;i++){
01878 d[i] = new double [ndofmas[i]];
01879 for (j=0;j<ndofmas[i];j++){
01880 d[i][j]=0.0;
01881 }
01882 }
01883 }
01884 else{
01885 for (i=0;i<ns;i++){
01886 for (j=0;j<ndofmas[i];j++){
01887 d[i][j]=0.0;
01888 }
01889 }
01890 }
01891
01892 lagrmultnodunknowns (out);
01893
01894
01895 for (i=0;i<ns;i++){
01896 for (j=0;j<ndofmas[i];j++){
01897 lhs[cndom[i][j]-1]=d[i][j];
01898 }
01899 }
01900
01901 }
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913 void seqfeti::lagrmultnodunknowns (FILE *out)
01914 {
01915 long i,j,k,l;
01916 double *r,*pp,*alpha,*u,*av;
01917
01918 av = new double [6];
01919 pp = new double [ndofmax];
01920 r = new double [ndofcp];
01921 alpha = new double [gsize];
01922 u = new double [gsize];
01923
01924
01925 nullv (r,ndofcp);
01926
01927 for (j=0;j<ns;j++){
01928 nullv (pp,ndofmax);
01929
01930 coarse_local (j,pp,lambda);
01931
01932
01933 subv (pp,ff[j],ndofmas[j]);
01934
01935
01936 if (smst==dense_matrix)
01937 break;
01938 if (smst==skyline_matrix)
01939 smsky[j].ldl_feti_sky (d[j],pp,nrbmdom[j],se[j],zero);
01940
01941
01942 local_coarse (j,d[j],r);
01943 }
01944
01945 cmulv (-1.0,r,ndofcp);
01946
01947 fprintf (out,"\n");
01948 for (i=0;i<ndofcp;i++){
01949 fprintf (out,"\n g %20.15le",r[i]);
01950 }
01951
01952 if (discont==1){
01953
01954
01955
01956 for (i=0;i<ndofcp;i++){
01957 r[i]-=b[i]+h[i]*lambda[i];
01958 }
01959 }
01960
01961
01962 mtxv (g,r,u,ndofcp,gsize);
01963
01964 mxv (invgg,u,alpha,gsize,gsize);
01965
01966
01967 for (i=0;i<ns;i++){
01968
01969 cmulv (-1.0,d[i],ndofmas[i]);
01970
01971 nullv (av,6);
01972 l=rbmadr[i];
01973 for (j=0;j<nrbmdom[i];j++){
01974 av[j] = alpha[l]; l++;
01975 }
01976
01977 for (j=0;j<ndofmas[i];j++){
01978 for (k=0;k<nrbmdom[i];k++){
01979 d[i][j]+=rbmdom[i][k*ndofmas[i]+j]*av[k];
01980 }
01981 }
01982 }
01983
01984 delete [] av;
01985 delete [] pp;
01986 delete [] r;
01987 delete [] alpha;
01988 delete [] u;
01989 }
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001 void seqfeti::assemble_ff (double *rhs)
02002 {
02003 long i,j;
02004
02005 if (ff==NULL){
02006 ff = new double* [ns];
02007 for (i=0;i<ns;i++){
02008 ff[i] = new double [ndofmas[i]];
02009 for (j=0;j<ndofmas[i];j++){
02010 ff[i][j]=0.0;
02011 }
02012 }
02013 }else{
02014 for (i=0;i<ns;i++){
02015 for (j=0;j<ndofmas[i];j++){
02016 ff[i][j]=0.0;
02017 }
02018 }
02019 }
02020
02021
02022 for (i=0;i<ns;i++){
02023 globloc (rhs,ff[i],cndom[i],ndofmas[i]);
02024 }
02025
02026 }
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037 void seqfeti::matrices_assembl (gmatrix *gm,FILE *out)
02038 {
02039 long i,j;
02040 time_t t1,t2,t3,t4,t5;
02041
02042
02043 det_ndofmax ();
02044
02045 t1 = time (NULL);
02046
02047
02048 subdomain_matrices (gm,out);
02049
02050 t2 = time (NULL);
02051
02052
02053 kernel (out);
02054
02055
02056 g_matrixsize (out);
02057
02058 if (gsize==0){
02059 print_err("gsize is 0",__FILE__,__LINE__,__func__);
02060 abort ();
02061 }
02062
02063 t3 = time (NULL);
02064
02065
02066 g_matrix (out);
02067
02068 t4 = time (NULL);
02069
02070 fprintf (out,"\n\n\n kontrola matice G \n");
02071 for (i=0;i<ndofcp;i++){
02072 fprintf (out,"\n %4ld",i);
02073 for (j=0;j<gsize;j++){
02074 fprintf (out," %f",g[i*gsize+j]);
02075 }
02076 }
02077
02078
02079 inverse_matrix_GG (out);
02080
02081 t5 = time (NULL);
02082
02083 fprintf (stdout,"\n time of subdomain matrices assembling %ld",t2-t1);
02084 fprintf (stdout,"\n time of kernel computation %ld",t3-t2);
02085 fprintf (stdout,"\n time of G matrix assembling %ld",t4-t3);
02086 fprintf (stdout,"\n time of inverse matrix G^T G computation %ld",t5-t4);
02087
02088
02089 matassem=1;
02090 }
02091
02092
02093
02094
02095
02096
02097
02098
02099
02100 void seqfeti::vectors_assembl (double *rhs,FILE *out)
02101 {
02102 time_t t1,t2,t3;
02103
02104 t1 = time (NULL);
02105
02106
02107 assemble_ff (rhs);
02108
02109 t2 = time (NULL);
02110
02111
02112 evector (out);
02113
02114 t3 = time (NULL);
02115
02116 fprintf (stdout,"\n time of f assembling %ld",t2-t1);
02117 fprintf (stdout,"\n time of e assembling %ld",t3-t2);
02118 }
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131 void seqfeti::solve_system (gtopology *top,gmatrix *gm,
02132 double *lhs,double *rhs,FILE *out)
02133 {
02134 long i,j;
02135
02136 if (matassem==0)
02137 matrices_assembl (gm,out);
02138
02139 vectors_assembl (rhs,out);
02140
02141 fprintf (stdout,"\n\n gsize %ld \n",gsize);
02142 fprintf (out,"\n ndofcp %ld",ndofcp);
02143 fprintf (stdout,"\n ndofcp %ld",ndofcp);
02144
02145 define_b (out);
02146 define_h (out);
02147
02148
02149 mpcg (out);
02150
02151 fprintf (out,"\n\n\n kontrola multiplikatoru \n");
02152 for (j=0;j<ndofcp;j++){
02153 fprintf (out,"\n %4ld %20.15le",j,lambda[j]);
02154 }
02155
02156
02157 nodalunknowns (lhs,out);
02158
02159 fprintf (out,"\n\n kontrola");
02160 for (i=0;i<ns;i++){
02161 fprintf (out,"\n");
02162 for (j=0;j<ndofmas[i];j++){
02163 fprintf (out,"\n displ %4ld %4ld %20.15le",i,j,d[i][j]);
02164 }
02165 }
02166
02167 }