00001 #include <math.h>
00002 #include "edsolver.h"
00003 #include "global.h"
00004 #include "globmat.h"
00005 #include "mechprint.h"
00006 #include "elemswitch.h"
00007
00008 #include "gmatrix.h"
00009 #include "eigsol.h"
00010 #include "vector.h"
00011
00012
00013 void solve_eigen_dynamics (double *x,double *w)
00014 {
00015 long i;
00016
00017
00018 stiffness_matrix (0);
00019
00020 mass_matrix (0);
00021
00022 switch (Mp->eigsol.teigsol){
00023 case inv_iteration:{
00024 inverse_iteration (x,w);
00025 break;
00026 }
00027 case subspace_it_jacobi:{
00028 subspace_iter_jac (x,w);
00029 break;
00030 }
00031 case subspace_it_gsortho:{
00032 subspace_iter_ortho (x,w);
00033 break;
00034 }
00035 case shifted_subspace_it_gsortho:{
00036 subspace_shift_iter_ortho (x,w);
00037 break;
00038 } default:{
00039 print_err("unknown eigenvalue solver is required",__FILE__,__LINE__,__func__);
00040 }
00041 }
00042
00043 print_eigenvalues (w);
00044 print_eigenvectors ();
00045
00046
00047
00048 print_init(-1, "wt");
00049 for (i=0;i<Mp->eigsol.neigv;i++){
00050
00051 compute_req_val(i);
00052
00053
00054
00055 print_step(i, i+1, w[i], NULL);
00056 }
00057 print_close();
00058
00059
00060
00061
00062 }
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 void inverse_iteration (double *x,double *w)
00074 {
00075 long i,ni,n;
00076 double nom,denom,zero,rho,prevrho,err,error;
00077 double *p,*z;
00078
00079
00080 ni=Mp->eigsol.nies;
00081
00082 err=Mp->eigsol.erres;
00083
00084 n=Ndofm;
00085
00086 zero=Mp->zero;
00087
00088
00089
00090 p = new double [n];
00091 z = new double [n];
00092
00093
00094
00095
00096
00097
00098 for (i=0;i<n;i++){
00099 x[i]=1.0;
00100 }
00101
00102
00103 Mmat->gmxv (x,z);
00104
00105
00106 for (i=0;i<ni;i++){
00107
00108 copyv (z,p,n);
00109
00110
00111
00112 Mp->ssle->solve_system (Gtm,Smat,x,p,Out);
00113
00114
00115 nom=ss(x,z,n);
00116
00117
00118 Mmat->gmxv (x,z);
00119
00120
00121 denom=ss(x,z,n);
00122 rho=nom/denom;
00123 denom=1.0/sqrt(denom);
00124
00125 cmulv (denom,z,n);
00126
00127 if (Mespr==1) printf ("\n iterace %4ld rho %e",i,rho);
00128 if (i!=0){
00129 error=(prevrho-rho)/prevrho;
00130 if (error<err){
00131 cmulv (denom,x,n); break;
00132 }
00133 }
00134 prevrho=rho;
00135 }
00136 Mp->eigsol.anies=i;
00137 Mp->eigsol.aerres=error;
00138 w[0]=rho;
00139
00140 delete [] z; delete [] p;
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151 void subspace_iter_ortho (double *x,double *w)
00152 {
00153 long i,j,k,l,ii,ni,n,nv,nrv;
00154 double err,error,maxerror,zero,nom,denom,alpha;
00155 double *z,*p,*wo;
00156
00157
00158 n=Ndofm;
00159
00160 ni=Mp->eigsol.nies;
00161
00162 err=Mp->eigsol.erres;
00163
00164 nrv=Mp->eigsol.neigv;
00165
00166 nv=Mp->eigsol.nev;
00167
00168 zero=Mp->zero;
00169
00170
00171
00172
00173 z = new double [n];
00174 p = new double [n];
00175 wo = new double [nv];
00176
00177
00178
00179
00180
00181 for (i=0;i<nv*n;i++){
00182 x[i]=1.0;
00183 }
00184
00185
00186
00187 for (i=0;i<ni;i++){
00188 for (j=0;j<nv;j++){
00189
00190
00191 Mmat->gmxv (x+j*n,z);
00192
00193 copyv (z,p,n);
00194
00195
00196
00197 Mp->ssle->solve_system (Gtm,Smat,x+j*n,z,Out);
00198
00199
00200 nom=ss(x+j*n,p,n);
00201
00202 Mmat->gmxv (x+j*n,z);
00203
00204 denom=ss(x+j*n,z,n);
00205 w[j]=nom/denom;
00206 }
00207
00208
00209 for (j=0;j<nv;j++){
00210
00211 Mmat->gmxv (x+j*n,z);
00212
00213 for (k=0;k<j;k++){
00214 alpha=ss(x+k*n,z,n);
00215 for (l=0;l<n;l++){
00216 x[j*n+l]-=alpha*x[k*n+l];
00217 }
00218 }
00219
00220 Mmat->gmxv (x+j*n,z);
00221
00222 alpha=ss(x+j*n,z,n);
00223 alpha=1.0/sqrt(alpha);
00224 cmulv (alpha,x+j*n,n);
00225 }
00226
00227
00228 if (i!=0){
00229 ii=0; maxerror=0.0;
00230 for (j=0;j<nrv;j++){
00231 error=(wo[j]-w[j])/w[j];
00232 if (error>maxerror) maxerror=error;
00233 if (error<err) ii++;
00234 }
00235
00236 if (Mespr==1) printf ("\n iterace %4ld pocet doiterovanych tvaru %4ld max. odchylka %e",i,ii,maxerror);
00237
00238
00239 if (ii==nrv){
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265 break;
00266 }
00267 }
00268
00269 for (j=0;j<nv;j++){
00270 wo[j]=w[j];
00271 }
00272 }
00273
00274 Mp->eigsol.anies=i;
00275 Mp->eigsol.aerres=maxerror;
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 delete [] z;
00301 delete [] wo;
00302 delete [] p;
00303 }
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316 void subspace_iter_jac (double *x,double *w)
00317 {
00318 long i,j,k,l,ii,ni,n,nv,nrv,nij,njacthr,stop,*ind;
00319 double err,error,maxerror,zero,alpha;
00320 double *z,*p,*wo,*redstiff,*redmass,*redeigv,*jacthr,*aux,*invec;
00321
00322
00323 n=Ndofm;
00324
00325 ni=Mp->eigsol.nies;
00326
00327 err=Mp->eigsol.erres;
00328
00329 nrv=Mp->eigsol.neigv;
00330
00331 nv=Mp->eigsol.nev;
00332
00333 nij=Mp->eigsol.nijmr;
00334
00335 njacthr=Mp->eigsol.njacthr;
00336
00337 zero=Mp->zero;
00338
00339 stop=0;
00340
00341 jacthr=Mp->eigsol.jacthr;
00342
00343 z = new double [n*nv];
00344 p = new double [n*nv];
00345 wo = new double [nv];
00346 redstiff = new double [nv*nv];
00347 redmass = new double [nv*nv];
00348 redeigv = new double [nv*nv];
00349 aux = new double [nv];
00350 invec = new double [n];
00351 ind = new long [n];
00352
00353 for (i=0;i<n;i++){
00354
00355 }
00356 for (i=0;i<n;i++){
00357 alpha=invec[i]; k=i;
00358 for (j=i+1;j<n;j++){
00359 if (alpha<invec[j]){
00360 alpha=invec[j]; k=j;
00361 }
00362 }
00363 ind[i]=k;
00364 alpha=invec[i];
00365 invec[i]=invec[k];
00366 invec[k]=alpha;
00367 }
00368
00369
00370
00371
00372 Mp->ssle->solve_system (Gtm,Smat,z,z,Out);
00373
00374
00375 for (i=0;i<nv*n;i++){
00376 x[i]=0.0;
00377 }
00378 for (i=0;i<n;i++){
00379 x[i]=1.0;
00380 }
00381 for (i=1;i<nv;i++){
00382 x[i*n+ind[i]]=1.0;
00383 }
00384
00385
00386 for (j=0;j<nv;j++){
00387
00388 Mmat->gmxv (x+j*n,z+j*n);
00389 }
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 for (i=0;i<ni;i++){
00410
00411
00412 for (j=0;j<nv;j++){
00413
00414 copyv (z+j*n,p+j*n,n);
00415
00416
00417
00418 Mp->ssle->solve_system (Gtm,Smat,x+j*n,z+j*n,Out);
00419 }
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435 mtxmccr (x,p,redstiff,n,nv,nv);
00436
00437
00438 for (j=0;j<nv;j++){
00439
00440 Mmat->gmxv (x+j*n,z+j*n);
00441 }
00442
00443
00444
00445 mtxmccr (x,z,redmass,n,nv,nv);
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459 gen_jacobi (redstiff,redmass,redeigv,w,nv,nij,jacthr,njacthr,zero);
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471 for (j=0;j<n;j++){
00472 for (k=0;k<nv;k++){
00473 aux[k]=0.0;
00474 for (l=0;l<nv;l++){
00475 aux[k]+=z[l*n+j]*redeigv[l*nv+k];
00476 }
00477 }
00478 for (k=0;k<nv;k++){
00479 z[k*n+j]=aux[k];
00480 }
00481 }
00482
00483
00484
00485
00486 double norm = sqrt(ss (z,z,n*nv));
00487 cmulv (1.0/norm,z,n*nv);
00488
00489
00490
00491 if (i!=0){
00492 ii=0; maxerror=0.0;
00493 for (j=0;j<nrv;j++){
00494 error=(wo[j]-w[j])/w[j];
00495 if (error>maxerror) maxerror=error;
00496 if (error<err) ii++;
00497 }
00498
00499 if (Mespr==1) fprintf (stdout,"\n iterace %4ld pocet doiterovanych tvaru %4ld max. odchylka %e",i,ii,maxerror);
00500
00501 if (ii==nrv) stop=1;
00502 }
00503
00504 for (j=0;j<nv;j++){
00505 wo[j]=w[j];
00506 }
00507
00508 if (stop==1) break;
00509 }
00510
00511 Mp->eigsol.anies=i;
00512 Mp->eigsol.aerres=maxerror;
00513
00514 delete [] z;
00515 }
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527 void subspace_shift_iter_ortho (double *x,double *w)
00528 {
00529 long i,j,k,l,ii,ni,n,nv,nrv;
00530 double err,error,maxerror,zero,nom,denom,alpha,shift;
00531 double *z,*p,*wo;
00532
00533
00534 n=Ndofm;
00535
00536 ni=Mp->eigsol.nies;
00537
00538 err=Mp->eigsol.erres;
00539
00540 nrv=Mp->eigsol.neigv;
00541
00542 nv=Mp->eigsol.nev;
00543
00544 zero=Mp->zero;
00545
00546 shift = Mp->eigsol.shift;
00547
00548 z = new double [n];
00549 p = new double [n];
00550 wo = new double [nv];
00551
00552
00553 for (i=0;i<nv*n;i++){
00554 x[i]=1.0;
00555 }
00556
00557
00558 Smat->addgm (shift,*Mmat);
00559
00560
00561 for (i=0;i<ni;i++){
00562 for (j=0;j<nv;j++){
00563
00564 Mmat->gmxv (x+j*n,z);
00565
00566 copyv (z,p,n);
00567
00568
00569
00570 Mp->ssle->solve_system (Gtm,Smat,x+j*n,z,Out);
00571
00572
00573 nom=ss(x+j*n,p,n);
00574
00575 Mmat->gmxv (x+j*n,z);
00576
00577 denom=ss(x+j*n,z,n);
00578 w[j]=nom/denom;
00579 }
00580
00581
00582 for (j=0;j<nv;j++){
00583
00584 Mmat->gmxv (x+j*n,z);
00585
00586 for (k=0;k<j;k++){
00587 alpha=ss(x+k*n,z,n);
00588 for (l=0;l<n;l++){
00589 x[j*n+l]-=alpha*x[k*n+l];
00590 }
00591 }
00592
00593 Mmat->gmxv (x+j*n,z);
00594
00595 alpha=ss(x+j*n,z,n);
00596 alpha=1.0/sqrt(alpha);
00597 cmulv (alpha,x+j*n,n);
00598 }
00599
00600
00601 if (i>5){
00602 ii=0; maxerror=0.0;
00603 for (j=0;j<nrv;j++){
00604 error=fabs(wo[j]-w[j])/w[j];
00605
00606 if (error>maxerror)
00607 maxerror=error;
00608 if (error<err)
00609 ii++;
00610 }
00611
00612 if (Mespr==1) printf ("\n iterace %4ld pocet doiterovanych tvaru %4ld max. odchylka %e",i,ii,maxerror);
00613
00614
00615 if (ii==nrv){
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641 break;
00642 }
00643 }
00644
00645 for (j=0;j<nv;j++){
00646 wo[j]=w[j];
00647 }
00648 }
00649
00650 Mp->eigsol.anies=i;
00651 Mp->eigsol.aerres=maxerror;
00652
00653 for (i=0;i<nv;i++){
00654 w[i]-=shift;
00655 }
00656
00657
00658 delete [] z;
00659 delete [] wo;
00660 delete [] p;
00661 }