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