00001 #include <math.h>
00002 #include "pedsolver.h"
00003 #include "pglobal.h"
00004 #include "xfile.h"
00005 #include "mpi.h"
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 void parallel_solution_eigen_dynamics (double *x,double *w)
00017 {
00018 long i;
00019
00020
00021 stiffness_matrix (0);
00022
00023 mass_matrix (0);
00024
00025 switch (Mp->eigsol.teigsol){
00026 case inv_iteration:{
00027 inverse_iteration (x,w);
00028 break;
00029 }
00030 case subspace_it_jacobi:{
00031 subspace_iter_jac (x,w);
00032 break;
00033 }
00034 case subspace_it_gsortho:{
00035 subspace_iter_ortho (x,w);
00036 break;
00037 }
00038 default:{
00039 par_print_err(Myrank+1, proc_name, "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 void parallel_inverse_iteration (double *x,double *w)
00073 {
00074 long i,ni,n;
00075 double nom,denom,zero,rho,prevrho,err,error;
00076 double *p,*z;
00077
00078
00079 ni=Mp->eigsol.nies;
00080
00081 err=Mp->eigsol.erres;
00082
00083 n=Ndofm;
00084
00085 zero=Mp->zero;
00086
00087
00088 p = new double [n];
00089 z = new double [n];
00090
00091
00092 for (i=0;i<n;i++){
00093 x[i]=1.0;
00094 }
00095
00096 Mmat->gmxv (x,z);
00097
00098
00099 for (i=0;i<ni;i++){
00100
00101 copyv (z,p,n);
00102
00103 Psol->par_linear_solver (Gtm,Smat,x,p,Out,Mespr);
00104
00105
00106 nom=Psol->pss (x,z,Out);
00107
00108
00109 Mmat->gmxv (x,z);
00110
00111
00112 denom=Psol->pss(x,z,Out);
00113 rho=nom/denom;
00114 denom=1.0/sqrt(denom);
00115
00116 cmulv (denom,z,n);
00117
00118 if (Mespr==1) printf ("\n iterace %4ld rho %e",i,rho);
00119 if (i!=0){
00120 error=(prevrho-rho)/prevrho;
00121 if (error<err){
00122 cmulv (denom,x,n);
00123 break;
00124 }
00125 }
00126 prevrho=rho;
00127 }
00128 Mp->eigsol.anies=i;
00129 Mp->eigsol.aerres=error;
00130 w[0]=rho;
00131
00132 delete [] z; delete [] p;
00133 }
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145 void parallel_subspace_iter_ortho (double *x,double *w)
00146 {
00147 long i,j,k,l,ii,ni,n,nv,nrv;
00148 double err,error,maxerror,zero,nom,denom,alpha;
00149 double *z,*p,*wo;
00150
00151
00152 n=Ndofm;
00153
00154 ni=Mp->eigsol.nies;
00155
00156 err=Mp->eigsol.erres;
00157
00158 nrv=Mp->eigsol.neigv;
00159
00160 nv=Mp->eigsol.nev;
00161
00162 zero=Mp->zero;
00163
00164
00165 z = new double [n];
00166 p = new double [n];
00167 wo = new double [nv];
00168
00169
00170 for (i=0;i<nv*n;i++){
00171 x[i]=1.0;
00172 }
00173
00174
00175
00176 for (i=0;i<ni;i++){
00177 for (j=0;j<nv;j++){
00178
00179 Mmat->gmxv (x+j*n,z);
00180
00181 copyv (z,p,n);
00182
00183 Psol->par_linear_solver (Gtm,Smat,x+j*n,p,Out,Mespr);
00184
00185
00186 nom=Psol->pss(x+j*n,z,Out);
00187
00188 Mmat->gmxv (x+j*n,z);
00189
00190 denom=Psol->pss(x+j*n,z,Out);
00191 w[j]=nom/denom;
00192 }
00193
00194
00195 for (j=0;j<nv;j++){
00196
00197
00198
00199 Mmat->gmxv (x+j*n,z);
00200
00201 for (k=0;k<j;k++){
00202 alpha=Psol->pss(x+k*n,z,Out);
00203 for (l=0;l<n;l++){
00204 x[j*n+l]-=alpha*x[k*n+l];
00205 }
00206 }
00207
00208
00209 Mmat->gmxv (x+j*n,z);
00210
00211 alpha=Psol->pss(x+j*n,z,Out);
00212 alpha=1.0/sqrt(alpha);
00213 cmulv (alpha,x+j*n,n);
00214 }
00215
00216
00217 if (i!=0){
00218 ii=0; maxerror=0.0;
00219 for (j=0;j<nrv;j++){
00220 error=(wo[j]-w[j])/w[j];
00221 if (error>maxerror) maxerror=error;
00222 if (error<err) ii++;
00223 }
00224
00225 if (Mespr==1) printf ("\n iterace %4ld pocet doiterovanych tvaru %4ld max. odchylka %e",i,ii,maxerror);
00226
00227
00228 if (ii==nrv){
00229
00230 break;
00231 }
00232 }
00233
00234 for (j=0;j<nv;j++){
00235 wo[j]=w[j];
00236 }
00237 }
00238
00239 Mp->eigsol.anies=i;
00240 Mp->eigsol.aerres=maxerror;
00241
00242 delete [] z;
00243 delete [] wo;
00244 delete [] p;
00245 }
00246