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 smst=skyline_matrix;
00012
00013
00014
00015 ns=0;
00016
00017 ense=0;
00018
00019
00020
00021 ndofmax=0;
00022
00023 ndofcp=0;
00024
00025 thresh=0.0;
00026
00027 nicg=0;
00028
00029 anicg=0;
00030
00031 errcg=0.0;
00032
00033 aerrcg=0.0;
00034
00035 zero=1.0e-100;
00036
00037
00038 prec=noprecond;
00039
00040 gsize=0;
00041
00042 ndofmax=0;
00043
00044 ndofcp=0;
00045
00046
00047
00048 ndofmas=NULL;
00049
00050
00051 nsid=NULL;
00052
00053
00054
00055 cndom=NULL;
00056
00057
00058 ncdofd=NULL;
00059
00060
00061
00062 edofs=NULL;
00063
00064
00065 nrbmdom=NULL;
00066
00067
00068 rbmdom=NULL;
00069
00070
00071 rbmadr=NULL;
00072
00073
00074 se=NULL;
00075
00076
00077 g=NULL;
00078
00079
00080 invgg = NULL;
00081
00082
00083 ff = NULL;
00084
00085
00086 e = NULL;
00087
00088
00089 d = NULL;
00090
00091
00092 lambda = NULL;
00093
00094
00095
00096 ncbool=NULL;
00097
00098
00099 booldatar=NULL;
00100
00101
00102 booldatac=NULL;
00103
00104
00105 booldata=NULL;
00106
00107
00108 ccn=NULL;
00109
00110
00111
00112 smsky=NULL;
00113
00114 smdm=NULL;
00115
00116 }
00117
00118 seqfeti::~seqfeti ()
00119 {
00120 long i;
00121
00122
00123
00124
00125 for (i=0;i<ns;i++){
00126
00127 delete [] ccn[i];
00128
00129 delete [] booldata[i];
00130
00131 delete [] booldatac[i];
00132
00133 delete [] booldatar[i];
00134
00135 delete [] se[i];
00136
00137 delete [] rbmdom[i];
00138
00139 delete [] edofs[i];
00140
00141 delete [] cndom[i];
00142 }
00143
00144
00145 delete [] ccn;
00146
00147 delete [] booldata;
00148
00149 delete [] booldatac;
00150
00151 delete [] booldatar;
00152
00153 delete [] se;
00154
00155 delete [] rbmdom;
00156
00157 delete [] edofs;
00158
00159 delete [] cndom;
00160
00161
00162
00163
00164
00165
00166 delete [] ncbool;
00167
00168 delete [] g;
00169
00170 delete [] rbmadr;
00171
00172 delete [] nrbmdom;
00173
00174 delete [] ncdofd;
00175
00176 delete [] nsid;
00177
00178 delete [] ndofmas;
00179
00180
00181 for (i=0;i<ns;i++){
00182 smsky[i].~skyline ();
00183 smdm[i].~densemat ();
00184 }
00185 }
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 void seqfeti::read (gtopology *top,XFILE *in)
00196 {
00197
00198
00199
00200
00201 xfscanf (in,"%m",&fetiimplem_kwdset,&fetiimpl);
00202
00203 if (fetiimpl!=boolean_matrices && fetiimpl!=nonredundant && fetiimpl!=redundant){
00204 print_err("unknown type of sequential FETI implementation",__FILE__,__LINE__,__func__);
00205 }
00206
00207
00208
00209
00210
00211
00212 xfscanf (in,"%ld %ld %le %ld %le",&ns,&ense,&thresh,&nicg,&errcg);
00213
00214 xfscanf (in,"%d",(int*)&prec);
00215
00216
00217
00218
00219
00220 if (fetiimpl==nonredundant || fetiimpl==redundant){
00221
00222 top->rst=1;
00223 }
00224 if (fetiimpl==boolean_matrices){
00225
00226
00227 top->rst=2;
00228 }
00229 if (prec==feti_dirichlet){
00230
00231
00232 top->cngen=2;
00233 }
00234
00235
00236 top->cngen=0;
00237
00238 }
00239
00240
00241
00242
00243
00244
00245
00246
00247 void seqfeti::print (FILE *out)
00248 {
00249
00250 fprintf (out,"%d",fetiimpl);
00251
00252
00253
00254
00255
00256
00257 fprintf (out,"%ld %ld %le %ld %le\n",ns,ense,thresh,nicg,errcg);
00258
00259 fprintf (out,"%d\n",prec);
00260 }
00261
00262
00263
00264
00265
00266
00267
00268
00269 void seqfeti::read_booldata (XFILE *in)
00270 {
00271 long i,j;
00272
00273
00274 ncbool = new long [ns];
00275 booldatar = new long* [ns];
00276 booldatac = new long* [ns];
00277 booldata = new double* [ns];
00278
00279 for (i=0;i<ns;i++){
00280
00281
00282
00283 xfscanf (in,"%ld",ncbool+i);
00284
00285
00286 booldatar[i] = new long [ncbool[i]];
00287 booldatac[i] = new long [ncbool[i]];
00288 booldata[i] = new double [ncbool[i]];
00289
00290 for (j=0;j<ncbool[i];j++){
00291 xfscanf (in,"%ld %ld %lf",&booldatar[i][j],&booldatac[i][j],&booldata[i][j]);
00292 }
00293 }
00294
00295 }
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310 void seqfeti::initiate (seqselnodes *selnodfeti,gtopology *top,FILE *out)
00311 {
00312 long i,j,k,ndofn,dof,nid;
00313 long *red;
00314 red = new long [ns+1];
00315
00316
00317 ndofcp = selnodfeti->tndofsn;
00318
00319
00320 ncdofd = new long [ns];
00321 for (i=0;i<ns;i++){
00322 ncdofd[i] = selnodfeti->snndofmas[i];
00323 }
00324
00325
00326 edofs = new long* [ns];
00327 for (i=0;i<ns;i++){
00328 edofs[i] = new long [ncdofd[i]];
00329 }
00330
00331 for (i=0;i<ns;i++){
00332 for (j=0;j<ncdofd[i];j++){
00333 edofs[i][j] = selnodfeti->lndofmas[i][j];
00334 }
00335 }
00336
00337
00338
00339
00340
00341
00342 red[0]=0;
00343 nid=0;
00344
00345 for (i=0;i<ns;i++){
00346 red[i+1]=0;
00347
00348 for (j=0;j<top->stop->nnsd[i];j++){
00349 ndofn = top->give_ndofn (nid);
00350 for (k=0;k<ndofn;k++){
00351 dof = top->give_dof (nid,k);
00352 if (dof>red[i+1]){
00353 red[i+1]=dof;
00354 }
00355 }
00356 nid++;
00357 }
00358 }
00359
00360
00361
00362
00363 for (i=0;i<ns;i++){
00364 for (j=0;j<ncdofd[i];j++){
00365 if (edofs[i][j]>0)
00366 edofs[i][j] -= red[i];
00367 else
00368 edofs[i][j] += red[i];
00369 }
00370 }
00371
00372
00373 fprintf (out,"\n\n\n check of the array edofs \n");
00374 for (i=0;i<ns;i++){
00375 fprintf (out,"\n domain %ld",i);
00376 for (j=0;j<ncdofd[i];j++){
00377 fprintf (out,"\n edofs %6ld %6ld",j,edofs[i][j]);
00378 }
00379 }
00380
00381
00382
00383
00384
00385
00386 wscalmat = new double [ndofcp];
00387 for (i = 0; i < ndofcp; i++){
00388
00389 }
00390
00391
00392
00393
00394
00395 ccn = new long* [ns];
00396 for (i=0;i<ns;i++){
00397 ccn[i] = new long [ncdofd[i]];
00398 }
00399
00400 for (i=0;i<ns;i++){
00401 for (j=0;j<ncdofd[i];j++){
00402 ccn[i][j] = selnodfeti->cndofmas[i][j];
00403 }
00404 }
00405
00406 delete [] red;
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 }
00421
00422
00423
00424
00425
00426
00427
00428
00429 void seqfeti::det_ndofmax ()
00430 {
00431 long i,j;
00432
00433 if (fetiimpl==nonredundant || fetiimpl==redundant){
00434 ndofmax=0;
00435 for (i=0;i<ns;i++){
00436 if (ndofmax<ndofmas[i])
00437 ndofmax=ndofmas[i];
00438 }
00439 }
00440
00441 if (fetiimpl==boolean_matrices){
00442 ndofmax=0;
00443 ndofcp=0;
00444 for (i=0;i<ns;i++){
00445 for (j=0;j<ncbool[i];j++){
00446 if (ndofmax<booldatac[i][j])
00447 ndofmax=booldatac[i][j];
00448 if (ndofcp<booldatar[i][j])
00449 ndofcp=booldatar[i][j];
00450 }
00451 }
00452 }
00453 }
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468 void seqfeti::assemble_subdom_unknowns (gtopology *top,FILE *out)
00469 {
00470 long i,j,k,l,m,ndofn;
00471
00472
00473 nsid = new long [top->nn];
00474 l=0;
00475 for (i=0;i<ns;i++){
00476 for (j=0;j<top->stop->nnsd[i];j++){
00477 nsid[l]=i;
00478 l++;
00479 }
00480 }
00481
00482
00483
00484 ndofmas = new long [ns];
00485 for (i=0;i<ns;i++){
00486 ndofmas[i]=0;
00487 }
00488
00489 for (i=0;i<top->nn;i++){
00490
00491 ndofn=top->give_ndofn (i);
00492
00493 l=nsid[i];
00494 for (j=0;j<ndofn;j++){
00495
00496 k=top->give_dof (i,j);
00497 if (k>0)
00498 ndofmas[l]++;
00499 }
00500 }
00501
00502
00503 cndom = new long* [ns];
00504 for (i=0;i<ns;i++){
00505 cndom[i] = new long [ndofmas[i]];
00506 ndofmas[i]=0;
00507 }
00508
00509 for (i=0;i<top->nn;i++){
00510
00511 ndofn=top->give_ndofn (i);
00512
00513 l=nsid[i];
00514 for (j=0;j<ndofn;j++){
00515
00516 k=top->give_dof (i,j);
00517 if (k>0){
00518 cndom[l][ndofmas[l]]=k;
00519 ndofmas[l]++;
00520 }
00521 }
00522 }
00523
00524
00525
00526
00527
00528 if (prec == feti_dirichlet){
00529 for (i=0;i<ns;i++){
00530 for (j=0;j<ndofmas[i];j++){
00531 l=cndom[i][j];
00532 m=j;
00533 for (k=j;k<ndofmas[i];k++){
00534 if (l>cndom[i][k]){
00535 l=cndom[i][k];
00536 m=k;
00537 }
00538 }
00539 cndom[i][m]=cndom[i][j];
00540 cndom[i][j]=l;
00541 }
00542 }
00543 }
00544
00545
00546 fprintf (out,"\n\n kontrola pole ndofmas");
00547 for (i=0;i<ns;i++){
00548 fprintf (out,"\n ndofmas %6ld %6ld",i,ndofmas[i]);
00549 }
00550 fprintf (out,"\n\n kontrola pole cndom");
00551 for (i=0;i<ns;i++){
00552 fprintf (out,"\n domain %ld",i);
00553 for (j=0;j<ndofmas[i];j++){
00554 fprintf (out,"\n dom %6ld dof %6ld %6ld",i,j,cndom[i][j]);
00555 }
00556 }
00557
00558
00559
00560 jum = new double* [ns];
00561 for (i=0;i<ns;i++){
00562 jum[i] = new double [ndofmas[i]];
00563 }
00564
00565 for (i=0;i<ns;i++){
00566 for (j=0;j<ndofmas[i];j++){
00567 jum[i][j]=0.0;
00568 }
00569 }
00570
00571 }
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582 void seqfeti::subdomain_matrices (gmatrix *gm,FILE *out)
00583 {
00584 long i,j,k;
00585
00586 if (smst==dense_matrix){
00587 smdm = new densemat [ns];
00588 for (i=0;i<ns;i++){
00589
00590 }
00591 }
00592 if (smst==skyline_matrix){
00593 smsky = new skyline [ns];
00594 for (i=0;i<ns;i++){
00595 smsky[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofmas[i],cndom[i]);
00596 }
00597 }
00598
00599
00600 if (fetiimpl==nonredundant || fetiimpl==redundant){
00601 long *aux;
00602 ndofprec = new long [ns];
00603 cnprec = new long* [ns];
00604 cpreccn = new long* [ns];
00605
00606 for (i=0;i<ns;i++){
00607 aux = new long [ndofmas[i]];
00608 for (j=0;j<ndofmas[i];j++){
00609 aux[j]=0;
00610 }
00611
00612 for (j=0;j<ncdofd[i];j++){
00613 k=edofs[i][j];
00614 if (k<0)
00615 k=0-k-1;
00616 else
00617 k--;
00618 if (k>-1)
00619 aux[k]++;
00620 }
00621
00622 ndofprec[i]=0;
00623 for (j=0;j<ndofmas[i];j++){
00624 if (aux[j]>0)
00625 ndofprec[i]++;
00626 }
00627
00628 fprintf (out,"\n\n aux domena %ld\n",i);
00629 for (j=0;j<ndofmas[i];j++){
00630 fprintf (out," %ld",aux[j]);
00631 }
00632
00633
00634 cnprec[i] = new long [ndofprec[i]];
00635 cpreccn[i] = new long [ndofprec[i]];
00636 k=0;
00637 for (j=0;j<ndofmas[i];j++){
00638 if (aux[j]>0){
00639 cnprec[i][k]=cndom[i][j];
00640 cpreccn[i][k]=j;
00641 k++;
00642 }
00643 }
00644
00645 delete [] aux;
00646 }
00647
00648 fprintf (out,"\n\n\n\n kontrola ve funkci subdomain_matrices (file %s, line %d)\n\n",__FILE__,__LINE__);
00649 for (i=0;i<ns;i++){
00650 fprintf (out,"\n domena %5ld ndofprec %ld\n",i,ndofprec[i]);
00651 for (j=0;j<ndofprec[i];j++){
00652 fprintf (out," %ld",cnprec[i][j]);
00653 }
00654 fprintf (out,"\n");
00655 for (j=0;j<ndofprec[i];j++){
00656 fprintf (out," %ld",cpreccn[i][j]);
00657 }
00658 }
00659 fprintf (out,"\n\n\n\n\n");
00660 }
00661
00662
00663 if (fetiimpl==nonredundant || fetiimpl==redundant){
00664 if (prec==feti_lumped){
00665 smscr = new symcomprow [ns];
00666 for (i=0;i<ns;i++){
00667
00668 gm->scr->select_submatrix (&smscr[i],ndofprec[i],cnprec[i]);
00669 }
00670 }
00671
00672 if (prec==feti_dirichlet){
00673 psmdm = new densemat [ns];
00674
00675 double *x,*y,*condvect;
00676 for (i=0;i<ns;i++){
00677 x = new double [ndofmas[i]];
00678 y = new double [ndofmas[i]];
00679 psmdm[i].a = new double [ndofprec[i]*ndofprec[i]];
00680 psmdm[i].negm=ndofprec[i]*ndofprec[i];
00681 psmdm[i].n=ndofprec[i];
00682 condvect=new double [ndofprec[i]];
00683
00684 smsky[i].ldlkon_sky (psmdm[i].a,condvect,x,y,ndofprec[i],1);
00685 smsky[i].~skyline ();
00686
00687 delete [] x;
00688 delete [] y;
00689 delete [] condvect;
00690 }
00691
00692 smsky = new skyline [ns];
00693 for (i=0;i<ns;i++){
00694 smsky[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofmas[i],cndom[i]);
00695 }
00696
00697 }
00698 }
00699
00700 if (fetiimpl==boolean_matrices){
00701 if (prec==feti_lumped){
00702 smskyprec = new skyline [ns];
00703 for (i=0;i<ns;i++){
00704 smskyprec[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofmas[i],cndom[i]);
00705 }
00706 }
00707 if (prec==feti_dirichlet){
00708
00709 ndofprec = new long[ns];
00710 for (i=0;i<ns;i++){
00711 ndofprec[i] = 1;
00712 for(j = 1; j < ncbool[i]; j++){
00713 if(booldatac[i][j] != booldatac[i][j-1]){
00714 ndofprec[i]++;
00715 }
00716 }
00717 }
00718
00719 psmdm = new densemat [ns];
00720
00721 double *x,*y,*condvect;
00722 for (i=0;i<ns;i++){
00723 x = new double [ndofmas[i]];
00724 y = new double [ndofmas[i]];
00725 psmdm[i].a = new double [ndofprec[i]*ndofprec[i]];
00726 psmdm[i].negm=ndofprec[i]*ndofprec[i];
00727 psmdm[i].n=ndofprec[i];
00728 condvect=new double [ndofprec[i]];
00729
00730 smsky[i].ldlkon_sky (psmdm[i].a,condvect,x,y,ndofprec[i],1);
00731 smsky[i].~skyline ();
00732
00733 delete [] x;
00734 delete [] y;
00735 delete [] condvect;
00736 }
00737
00738 smsky = new skyline [ns];
00739 for (i=0;i<ns;i++){
00740 smsky[i].assemble_from_scr (gm->scr->adr,gm->scr->ci,gm->scr->a,ndofmas[i],cndom[i]);
00741 }
00742 }
00743 }
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755 }
00756
00757
00758
00759
00760
00761
00762 void seqfeti::kernel (FILE *out)
00763 {
00764 long i,j;
00765
00766
00767 if (nrbmdom!=NULL){
00768 delete [] nrbmdom;
00769 }
00770 nrbmdom = new long [ns];
00771
00772
00773 if (se!=NULL){
00774 for (i=0;i<ns;i++){
00775 delete [] se[i];
00776 }
00777 delete [] se;
00778 }
00779 se = new long* [ns];
00780 for (i=0;i<ns;i++){
00781 se[i] = new long [ense];
00782 for (j=0;j<ense;j++){
00783 se[i][j]=-1;
00784 }
00785 }
00786
00787
00788 if (rbmdom!=NULL){
00789 for (i=0;i<ns;i++){
00790 delete [] rbmdom[i];
00791 }
00792 delete [] rbmdom;
00793 }
00794 rbmdom = new double* [ns];
00795 for (i=0;i<ns;i++){
00796 rbmdom[i] = new double [ense*ndofmas[i]];
00797 for (j=0;j<ense*ndofmas[i];j++){
00798 rbmdom[i][j]=0.0;
00799 }
00800 }
00801
00802
00803 if (smst==dense_matrix){
00804
00805 for (i=0;i<ns;i++){
00806 smdm[i].ker (rbmdom[i],nrbmdom[i],se[i],ense,thresh);
00807 }
00808 }
00809 if (smst==skyline_matrix){
00810
00811 for (i=0;i<ns;i++){
00812 smsky[i].ker (rbmdom[i],nrbmdom[i],se[i],ense,thresh,3);
00813
00814
00815 }
00816 }
00817
00818 }
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841 void seqfeti::boolean_matrix (long sdid,double *a)
00842 {
00843 long i,ri,ci;
00844
00845 for (i=0;i<ncbool[sdid];i++){
00846 ri=booldatar[sdid][i]-1;
00847 ci=booldatac[sdid][i]-1;
00848 a[ri*ndofmas[sdid]+ci]=booldata[sdid][i];
00849 }
00850 }
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863 void seqfeti::local_coarse (long nd,double *lv,double *cv)
00864 {
00865 if (fetiimpl==nonredundant || fetiimpl==redundant){
00866
00867 long i,j,k;
00868
00869 for (i=0;i<ncdofd[nd];i++){
00870 j=edofs[nd][i];
00871 k=ccn[nd][i]-1;
00872 if (j<0){
00873 j=0-j-1;
00874 cv[k]-=lv[j];
00875 }
00876 else{
00877 j=j-1;
00878 cv[k]+=lv[j];
00879 }
00880 }
00881 }
00882
00883 if (fetiimpl==boolean_matrices){
00884
00885 double *a,*av;
00886
00887
00888 a = new double [ndofcp*ndofmas[nd]];
00889 nullv (a,ndofcp*ndofmas[nd]);
00890
00891
00892 av = new double [ndofcp];
00893 nullv (av,ndofcp);
00894
00895
00896 boolean_matrix (nd,a);
00897
00898 mxv (a,lv,av,ndofcp,ndofmas[nd]);
00899
00900 addv (cv,av,ndofcp);
00901
00902 delete [] a;
00903 delete [] av;
00904 }
00905
00906 }
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918 void seqfeti::coarse_local (long nd,double *lv,double *cv)
00919 {
00920 if (fetiimpl==nonredundant || fetiimpl==redundant){
00921 long i,j,k;
00922
00923 for (i=0;i<ncdofd[nd];i++){
00924 k=ccn[nd][i]-1;
00925 j=edofs[nd][i];
00926 if (j<0){
00927 j=0-j-1;
00928 lv[j]-=cv[k];
00929 }
00930 else{
00931 j=j-1;
00932 lv[j]+=cv[k];
00933 }
00934 }
00935 }
00936
00937 if (fetiimpl==boolean_matrices){
00938
00939 double *a,*av;
00940
00941
00942 a = new double [ndofcp*ndofmas[nd]];
00943 nullv (a,ndofcp*ndofmas[nd]);
00944
00945
00946 av = new double [ndofmas[nd]];
00947 nullv (av,ndofmas[nd]);
00948
00949
00950 boolean_matrix (nd,a);
00951
00952 mtxv (a,cv,av,ndofcp,ndofmas[nd]);
00953
00954 addv (lv,av,ndofmas[nd]);
00955
00956 delete [] a;
00957 delete [] av;
00958 }
00959
00960 }
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975 void seqfeti::g_matrixsize (FILE *out)
00976 {
00977 long i;
00978
00979
00980 rbmadr = new long [ns+1];
00981 for (i=0;i<ns+1;i++){
00982 rbmadr[i]=0;
00983 }
00984
00985
00986 rbmadr[0]=0;
00987 for (i=0;i<ns;i++){
00988 rbmadr[i+1]=rbmadr[i]+nrbmdom[i];
00989 }
00990
00991
00992 gsize=rbmadr[ns];
00993
00994 }
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004 void seqfeti::g_matrix (FILE *out)
01005 {
01006 long i,k,l;
01007 double *av;
01008
01009
01010 if (g!=NULL)
01011 delete [] g;
01012 g = new double [ndofcp*gsize];
01013
01014
01015 av = new double [ndofcp];
01016
01017 for (i=0;i<ns;i++){
01018 for (k=0;k<nrbmdom[i];k++){
01019 nullv (av,ndofcp);
01020 local_coarse (i,rbmdom[i]+k*ndofmas[i],av);
01021 for (l=0;l<ndofcp;l++){
01022 g[l*gsize+rbmadr[i]+k]=0.0-av[l];
01023 }
01024 }
01025 }
01026
01027 delete [] av;
01028 }
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038 void seqfeti::evector (FILE *out)
01039 {
01040 long i,j,k;
01041
01042 if (e!=NULL)
01043 delete [] e;
01044 e = new double [gsize];
01045
01046 k=0;
01047 for (i=0;i<ns;i++){
01048 for (j=0;j<nrbmdom[i];j++){
01049 e[k]=0.0-ss(&rbmdom[i][j*ndofmas[i]],ff[i],ndofmas[i]);
01050 k++;
01051 }
01052 }
01053
01054 }
01055
01056
01057
01058
01059
01060
01061 void seqfeti::inverse_matrix_GG (FILE *out)
01062 {
01063 long i,j;
01064 double *ee,*gg;
01065
01066 if (invgg != NULL){
01067 delete [] invgg;
01068 }
01069 invgg = new double [gsize*gsize];
01070
01071 gg = new double [gsize*gsize];
01072 ee = new double [gsize*gsize];
01073
01074 for (i=0;i<gsize;i++){
01075 for (j=0;j<gsize;j++){
01076 ee[i*gsize+j]=0.0;
01077 }
01078 ee[i*gsize+i]=1.0;
01079 }
01080
01081 mtxm (g,g,gg,ndofcp,gsize,gsize);
01082
01083 fprintf (out,"\n\n\n kontrola matice HH \n");
01084 for (i=0;i<gsize;i++){
01085 fprintf (out,"\n %4ld ",i);
01086 for (j=0;j<gsize;j++){
01087 fprintf (out," %lf",gg[i*gsize+j]);
01088 }
01089 }
01090
01091 gemp (gg,invgg,ee,gsize,gsize,zero,1);
01092
01093
01094
01095
01096
01097 fprintf (out,"\n\n\n kontrola matice HHinv \n");
01098 for (i=0;i<gsize;i++){
01099 fprintf (out,"\n %4ld ",i);
01100 for (j=0;j<gsize;j++){
01101 fprintf (out," %lf",invgg[i*gsize+j]);
01102 }
01103 }
01104
01105 delete [] gg;
01106 delete [] ee;
01107 }
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118 void seqfeti::feti_projection (double *v)
01119 {
01120 long i;
01121 double *p,*q,*w;
01122
01123 p = new double [gsize];
01124 q = new double [gsize];
01125 w = new double [ndofcp];
01126
01127
01128 mtxv (g,v,p,ndofcp,gsize);
01129
01130
01131 mxv (invgg,p,q,gsize,gsize);
01132
01133
01134 mxv (g,q,w,ndofcp,gsize);
01135
01136
01137 for (i=0;i<ndofcp;i++){
01138 v[i]-=w[i];
01139 }
01140
01141 delete [] w; delete [] q; delete [] p;
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 void seqfeti::scaling(double *invect,double *outvect,long n,FILE *out)
01182 {
01183 double *aux;
01184 long i;
01185 aux = new double[n];
01186
01187 for(i = 0; i < n; i++){
01188 aux[i] = invect[i]/wscalmat[i];
01189 }
01190
01191
01192
01193
01194
01195
01196
01197
01198 for(i = 0; i < n; i++){
01199 outvect[i] = aux[i];
01200 }
01201
01202 delete []aux;
01203 }
01204
01205 void seqfeti::lscaling(double *invect,double *outvect,long ndom,FILE *out)
01206 {
01207 double *aux;
01208 long i;
01209
01210 aux = new double[nlwscalmat[ndom]];
01211
01212 for(i = 0; i < nlwscalmat[ndom]; i++){
01213 aux[i] = invect[i]/lwscalmat[ndom][i];
01214 }
01215
01216
01217 fprintf (out,"\n\n\n kontrola skalovani \n");
01218 for(i = 0; i < nlwscalmat[ndom]; i++){
01219 fprintf (out,"%ld pred %le po %le wscalmat %lf ccn %ld\n",i+1,invect[i],aux[i],lwscalmat[ndom][i],ccn[ndom][i]);
01220 }
01221
01222
01223 for(i = 0; i < nlwscalmat[ndom]; i++){
01224 outvect[i] = aux[i];
01225 }
01226
01227 delete []aux;
01228 }
01229
01230
01231
01232
01233
01234
01235 void seqfeti::lumpedprec (long nd,double *dd,double *pp)
01236 {
01237 long i;
01238 double *d,*p;
01239
01240 if (fetiimpl == nonredundant || fetiimpl==redundant){
01241 d = new double [ndofprec[nd]];
01242 p = new double [ndofprec[nd]];
01243
01244 for (i=0;i<ndofprec[nd];i++){
01245 d[i]=0.0;
01246 p[i]=0.0;
01247 }
01248
01249 for (i=0;i<ndofprec[nd];i++){
01250 d[i]=dd[cpreccn[nd][i]];
01251 }
01252
01253 smscr[nd].mxv_scr (d,p);
01254
01255 for (i=0;i<ndofprec[nd];i++){
01256 pp[cpreccn[nd][i]]=p[i];
01257 }
01258 delete [] p;
01259 delete [] d;
01260 }
01261 if(fetiimpl == boolean_matrices){
01262
01263
01264
01265
01266
01267 smskyprec[nd].mxv_sky (dd,pp);
01268
01269
01270 }
01271 }
01272
01273
01274
01275
01276 void seqfeti::dirichletprec (long nd,double *dd,double *pp)
01277 {
01278 long i,j;
01279 double *d,*p;
01280
01281 if (fetiimpl == nonredundant || fetiimpl==redundant){
01282 d = new double [ndofprec[nd]];
01283 p = new double [ndofprec[nd]];
01284
01285 for (i=0;i<ndofprec[nd];i++){
01286 d[i]=0.0;
01287 p[i]=0.0;
01288 }
01289
01290 for (i=0;i<ndofprec[nd];i++){
01291 d[i]=dd[cpreccn[nd][i]];
01292 }
01293
01294 psmdm[nd].mxv_dm (d,p);
01295
01296
01297 for (i=0;i<ndofprec[nd];i++){
01298 pp[cpreccn[nd][i]]=p[i];
01299 }
01300
01301 delete [] p;
01302 delete [] d;
01303 }
01304 if (fetiimpl == boolean_matrices){
01305 d = new double [ndofprec[nd]];
01306 p = new double [ndofprec[nd]];
01307 j = 0;
01308 for (i=0;i<ndofmas[nd];i++){
01309 if(dd[i] != 0){
01310 d[j]=dd[i];
01311 j++;
01312 }
01313 }
01314 psmdm[nd].mxv_dm (d,p);
01315
01316 j = 0;
01317 for (i=0;i<ndofmas[nd];i++){
01318 if(dd[i] != 0){
01319 pp[i]=p[j];
01320 j++;
01321 }
01322 else{
01323 pp[i] = 0;
01324 }
01325 }
01326 }
01327 }
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341 void seqfeti::mpcg (FILE *out)
01342 {
01343 long i,j;
01344 double nom,denom,alpha,beta;
01345 double *s,*sd,*r,*p,*pp,*invgge;
01346
01347 if (lambda==NULL)
01348 lambda = new double [ndofcp];
01349
01350 if (ndofmax==0){
01351 print_err("the attribute ndofmax is not determined",__FILE__,__LINE__,__func__);
01352 }
01353
01354
01355 sd = new double [ndofmax];
01356 pp = new double [ndofmax];
01357
01358
01359 s = new double [ndofcp];
01360
01361 r = new double [ndofcp];
01362
01363 p = new double [ndofcp];
01364
01365
01366
01367 invgge = new double [gsize];
01368
01369 mxv (invgg,e,invgge,gsize,gsize);
01370
01371 mxv (g,invgge,lambda,ndofcp,gsize);
01372 delete [] invgge;
01373
01374 nullv (r,ndofcp);
01375 for (j=0;j<ns;j++){
01376 nullv (sd,ndofmas[j]);
01377 coarse_local (j,sd,lambda);
01378
01379
01380
01381 subv (sd,ff[j],ndofmas[j]);
01382
01383
01384 if (smst==dense_matrix)
01385 break;
01386 if (smst==skyline_matrix)
01387 smsky[j].ldl_feti_sky (pp,sd,nrbmdom[j],se[j],zero);
01388
01389
01390 subv (pp,jum[j],ndofmas[j]);
01391
01392 local_coarse (j,pp,r);
01393 }
01394
01395 cmulv(-1.0,r,ndofcp);
01396
01397
01398
01399
01400
01401
01402
01403 feti_projection (r);
01404
01405
01406 for (i=0;i<ndofcp;i++){
01407 s[i]=r[i];
01408 }
01409
01410
01411 nom = ss (r,r,ndofcp);
01412
01413
01414
01415
01416
01417 for (i=0;i<nicg;i++){
01418
01419
01420
01421
01422 nullv (p,ndofcp);
01423 for (j=0;j<ns;j++){
01424 nullv (sd,ndofmax);
01425 coarse_local (j,sd,s);
01426
01427
01428 if (smst==dense_matrix)
01429 break;
01430 if (smst==skyline_matrix)
01431 smsky[j].ldl_feti_sky (pp,sd,nrbmdom[j],se[j],zero);
01432
01433 local_coarse (j,pp,p);
01434 }
01435
01436
01437 denom = ss (s,p,ndofcp);
01438
01439 fprintf (out,"\n MPCG in FETI, iterace %ld, kontrola citatele a jmenovatele alpha %e / %e",i,nom,denom);
01440
01441 if (fabs(denom)<zero){
01442 print_err("zero denominator of alpha in the modified conjugate gradient method",__FILE__,__LINE__,__func__);
01443 break;
01444 }
01445
01446
01447
01448
01449 alpha = nom/denom;
01450
01451
01452
01453
01454 for (j=0;j<ndofcp;j++){
01455 r[j]-=alpha*p[j];
01456 lambda[j]+=alpha*s[j];
01457 }
01458
01459 feti_projection (r);
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
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 denom = nom;
01541 if (fabs(denom)<zero){
01542 print_err("zero denominator of beta in the modified conjugate gradient method",__FILE__,__LINE__,__func__);
01543 break;
01544 }
01545
01546 nom=ss(r,r,ndofcp);
01547
01548 fprintf (stdout,"\n iteration %4ld norm r %e",i,sqrt(nom));
01549 fprintf (out,"\n iteration %4ld norm r %e",i,sqrt(nom));
01550
01551 fprintf (out,"\n kontrola citatele a jmenovatele pred betou %e / %e",nom,denom);
01552
01553 if (sqrt(nom)<errcg){
01554 break;
01555 }
01556
01557
01558
01559 beta = nom/denom;
01560
01561
01562 for (j=0;j<ndofcp;j++){
01563 s[j]=beta*s[j]+r[j];
01564 }
01565
01566 }
01567
01568
01569 anicg=i;
01570
01571 aerrcg=nom;
01572
01573 delete [] d;
01574 delete [] p;
01575 delete [] g;
01576
01577 delete [] pp;
01578 }
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
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
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
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098
02099
02100
02101
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114 void seqfeti::nodalunknowns (double *lhs, FILE *out)
02115 {
02116 long i,j;
02117
02118 if (d == NULL){
02119 d = new double* [ns];
02120 for (i=0;i<ns;i++){
02121 d[i] = new double [ndofmas[i]];
02122 for (j=0;j<ndofmas[i];j++){
02123 d[i][j]=0.0;
02124 }
02125 }
02126 }
02127
02128 lagrmultdispl (out);
02129
02130
02131 for (i=0;i<ns;i++){
02132 for (j=0;j<ndofmas[i];j++){
02133 lhs[cndom[i][j]-1]=d[i][j];
02134 }
02135 }
02136
02137 }
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149 void seqfeti::lagrmultdispl (FILE *out)
02150 {
02151 long i,j,k,l;
02152 double *r,*pp,*alpha,*u,*av;
02153
02154
02155 av = new double [6];
02156 pp = new double [ndofmax];
02157 r = new double [ndofcp];
02158 alpha = new double [gsize];
02159 u = new double [gsize];
02160
02161
02162 nullv (r,ndofcp);
02163
02164 for (j=0;j<ns;j++){
02165 nullv (pp,ndofmax);
02166
02167 coarse_local (j,pp,lambda);
02168
02169
02170 subv (pp,ff[j],ndofmas[j]);
02171
02172
02173 if (smst==dense_matrix)
02174 break;
02175 if (smst==skyline_matrix)
02176 smsky[j].ldl_feti_sky (d[j],pp,nrbmdom[j],se[j],zero);
02177
02178
02179 subv (d[j],jum[j],ndofmas[j]);
02180
02181
02182 local_coarse (j,d[j],r);
02183 }
02184
02185 fprintf (out,"\n");
02186 for (i=0;i<ndofcp;i++){
02187 fprintf (out,"\n g %20.15le",r[i]);
02188 }
02189
02190
02191
02192
02193
02194 mtxv (g,r,u,ndofcp,gsize);
02195 fprintf (out,"\n");
02196 for (i=0;i<gsize;i++){
02197 fprintf (out,"\n u %20.15le",u[i]);
02198 }
02199
02200
02201 mxv (invgg,u,alpha,gsize,gsize);
02202 fprintf (out,"\n");
02203 for (i=0;i<gsize;i++){
02204 fprintf (out,"\n alpha %20.15le",alpha[i]);
02205 }
02206
02207 cmulv (-1.0,alpha,gsize);
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221 for (i=0;i<ns;i++){
02222 cmulv (-1.0,d[i],ndofmas[i]);
02223
02224 nullv (av,6);
02225 l=rbmadr[i];
02226 for (j=0;j<nrbmdom[i];j++){
02227 av[j] = alpha[l]; l++;
02228 }
02229
02230 fprintf (out,"\n");
02231 for (j=0;j<nrbmdom[i];j++){
02232 fprintf (out,"\n test av %20.15le",av[j]);
02233 }
02234
02235 for (j=0;j<ndofmas[i];j++){
02236 for (k=0;k<nrbmdom[i];k++){
02237 d[i][j]+=rbmdom[i][k*ndofmas[i]+j]*av[k];
02238 }
02239 }
02240 }
02241
02242 delete [] av;
02243 delete [] pp;
02244 delete [] r;
02245 delete [] alpha;
02246 delete [] u;
02247 }
02248
02249
02250
02251
02252
02253
02254
02255
02256
02257
02258
02259 void seqfeti::assemble_ff (double *rhs)
02260 {
02261 long i,j;
02262
02263
02264 if (ff != NULL){
02265
02266 for (i=0;i<ns;i++){
02267 delete [] ff[i];
02268 }
02269 delete [] ff;
02270 }
02271
02272 ff = new double* [ns];
02273 for (i=0;i<ns;i++){
02274 ff[i] = new double [ndofmas[i]];
02275 for (j=0;j<ndofmas[i];j++){
02276 ff[i][j]=0.0;
02277 }
02278 }
02279
02280
02281 for (i=0;i<ns;i++){
02282 globloc (rhs,ff[i],cndom[i],ndofmas[i]);
02283 }
02284
02285 }
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302 void seqfeti::solve_system (gtopology *top,gmatrix *gm,
02303 double *lhs,double *rhs,FILE *out)
02304 {
02305 long i,j;
02306 time_t t1,t2,t3,t4,t5,t6;
02307
02308
02309 det_ndofmax ();
02310
02311
02312 assemble_ff (rhs);
02313
02314
02315 subdomain_matrices (gm,out);
02316
02317 t1 = time (NULL);
02318
02319
02320 kernel (out);
02321
02322 t2 = time (NULL);
02323
02324
02325 g_matrixsize (out);
02326
02327 fprintf (stdout,"\n\n gsize %ld \n",gsize);
02328
02329 if (gsize==0){
02330 fprintf (stdout,"\n\n gsize = 0 \n");
02331 fprintf (stderr,"\n\n gsize = 0 \n");
02332 abort ();
02333 }
02334
02335 fprintf (out,"\n ndofcp %ld",ndofcp);
02336
02337
02338 g_matrix (out);
02339
02340 fprintf (out,"\n\n\n kontrola matice H \n");
02341 for (i=0;i<ndofcp;i++){
02342 fprintf (out,"\n %4ld",i);
02343 for (j=0;j<gsize;j++){
02344 fprintf (out," %f",g[i*gsize+j]);
02345 }
02346 }
02347
02348
02349 evector (out);
02350
02351 fprintf (out,"\n\n\n kontrola vektoru e \n");
02352 for (j=0;j<gsize;j++){
02353 fprintf (out,"\n %4ld %f",j,e[j]);
02354 }
02355
02356 t3 = time (NULL);
02357
02358
02359 inverse_matrix_GG (out);
02360
02361
02362
02363 t4 = time (NULL);
02364
02365
02366 mpcg (out);
02367
02368
02369 fprintf (out,"\n\n\n kontrola multiplikatoru \n");
02370 for (j=0;j<ndofcp;j++){
02371 fprintf (out,"\n %4ld %20.15le",j,lambda[j]);
02372 }
02373
02374 t5 = time (NULL);
02375
02376 nodalunknowns (lhs,out);
02377
02378 fprintf (out,"\n\n kontrola");
02379 for (i=0;i<ns;i++){
02380 fprintf (out,"\n");
02381 for (j=0;j<ndofmas[i];j++){
02382 fprintf (out,"\n displ %4ld %4ld %20.15le",i,j,d[i][j]);
02383 }
02384 }
02385
02386
02387 t6 = time (NULL);
02388
02389
02390 fprintf (out,"\n\n evaluation of kernel %ld",t2-t1);
02391 fprintf (out,"\n\n matrix G computation %ld",t3-t2);
02392 fprintf (out,"\n\n decomposition of G %ld",t4-t3);
02393 fprintf (out,"\n\n modified conjug. gradient m. %ld",t5-t4);
02394 fprintf (out,"\n\n displacements computation %ld",t6-t5);
02395
02396 fprintf (stdout,"\n\n evaluation of kernel %ld",t2-t1);
02397 fprintf (stdout,"\n\n matrix G computation %ld",t3-t2);
02398 fprintf (stdout,"\n\n decomposition of G %ld",t4-t3);
02399 fprintf (stdout,"\n\n modified conjug. gradient m. %ld",t5-t4);
02400 fprintf (stdout,"\n\n displacements computation %ld",t6-t5);
02401
02402
02403 }