00001 #include "eigsol.h"
00002 #include "vector.h"
00003 #include <math.h>
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 void jacobi (double *a,double *evec,double *eval,long n,long ni,long &ani,double limit, long normalize)
00029 {
00030 long i,j,k,l,ii,jj,en,nn;
00031 double c,s,q,r,t,ai,aj,maxa;
00032
00033 if (normalize)
00034 {
00035 maxa = 0.0;
00036 for (i=0; i<n*n; i++)
00037 {
00038 if (a[0] != 0.0)
00039 {
00040 maxa = a[0];
00041 break;;
00042 }
00043 }
00044 if (maxa != 0.0)
00045 {
00046 for (i=0; i<n*n; i++)
00047 {
00048 if ((fabs(a[i]) > fabs(maxa)) && (a[i] != 0.0))
00049 maxa = a[i];
00050 }
00051 for (i=0; i<n*n; i++)
00052 a[i] /= maxa;
00053 }
00054 }
00055 en=(n*n-n)/2;
00056
00057 k=0;
00058 for (i=0;i<n;i++){
00059 for (j=0;j<n;j++){
00060 evec[k]=0.0; k++;
00061 }
00062 }
00063 for (i=0;i<n;i++){
00064 evec[i*n+i]=1.0;
00065 }
00066
00067
00068 for (k=0;k<ni;k++){
00069
00070 nn=0;
00071 for (i=0;i<n;i++){
00072 for (j=i+1;j<n;j++){
00073 if (fabs(a[i*n+j])<limit){
00074 nn++; continue;
00075 }
00076 else{
00077 q=(a[j*n+j]-a[i*n+i])/2.0/a[i*n+j];
00078 r=sqrt(1.0+q*q);
00079 if (q>=0.0) t=-q+r;
00080 else t=-q-r;
00081 r=sqrt(1.0+t*t);
00082 c=1.0/r;
00083 s=t*c;
00084
00085
00086 ii=i; jj=j;
00087 for (l=0;l<n;l++){
00088 ai=a[ii];
00089 aj=a[jj];
00090 a[ii]=ai*c-aj*s;
00091 a[jj]=ai*s+aj*c;
00092 ii+=n; jj+=n;
00093 }
00094
00095
00096 ii=i*n; jj=j*n;
00097 for (l=0;l<n;l++){
00098 ai=a[ii];
00099 aj=a[jj];
00100 a[ii]=ai*c-aj*s;
00101 a[jj]=ai*s+aj*c;
00102 ii++; jj++;
00103 }
00104
00105
00106 ii=i; jj=j;
00107 for (l=0;l<n;l++){
00108 ai=evec[ii];
00109 aj=evec[jj];
00110 evec[ii]=ai*c-aj*s;
00111 evec[jj]=ai*s+aj*c;
00112 ii+=n; jj+=n;
00113 }
00114
00115 }
00116 }
00117 }
00118 if (nn==en){
00119 break;
00120 }
00121 }
00122
00123
00124 ani=k;
00125
00126 for (i=0;i<n;i++){
00127 eval[i]=a[i*n+i];
00128 }
00129 if ((normalize) && (maxa != 0.0))
00130 {
00131 for (i=0;i<n;i++)
00132 eval[i]*=maxa;
00133 }
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 }
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 void gen_jacobi (double *a,double *b,double *x,double *w,long n,long ni,
00180 double *thresholds,long nthr,double zero)
00181 {
00182 long i,j,k,l,nt,numel,nze,stop;
00183 double discr,sol1,sol2,alpha,beta,limit,c0,c1,c2,denom;
00184
00185 numel=(n*n-n)/2;
00186 nt=0;
00187 limit=thresholds[nt];
00188 stop=0;
00189
00190
00191 for (i=0;i<n*n;i++){
00192 x[i]=0.0;
00193 }
00194 for (i=0;i<n;i++){
00195 x[i*n+i]=1.0;
00196 }
00197
00198
00199 for (i=0;i<ni;i++){
00200 fprintf (stdout,"\n\n JM iterace cislo %ld",i);
00201
00202
00203 nze=0;
00204 for (j=0;j<n;j++){
00205 for (k=j+1;k<n;k++){
00206 if (fabs(a[j*n+k])<limit && fabs(b[j*n+k])<limit) nze++;
00207 else{
00208
00209 c2=a[j*n+k]*b[j*n+j]-a[j*n+j]*b[j*n+k];
00210 c1=a[k*n+k]*b[j*n+j]-a[j*n+j]*b[k*n+k];
00211 c0=a[k*n+k]*b[j*n+k]-a[j*n+k]*b[k*n+k];
00212
00213 if (fabs(c2)<zero){
00214
00215 if (fabs(c1)<zero){
00216 fprintf (stderr,"\n\n non-solvable equation in function jacobi.\n");
00217 }
00218 alpha=(0.0-c0)/c1;
00219 }
00220 else{
00221 discr=c1*c1-4.0*c2*c0;
00222 if (discr<0.0){
00223 fprintf (stderr,"\n\n negative discriminant in quadratic equation in function jacobi.\n");
00224 }
00225 discr=sqrt(discr);
00226 sol1=(-1.0*c1+discr)/2.0/c2;
00227 sol2=(-1.0*c1-discr)/2.0/c2;
00228
00229
00230 if (fabs(sol1)<fabs(sol2)) alpha=sol2;
00231 else alpha=sol1;
00232
00233 }
00234 denom=a[k*n+k]+alpha*a[j*n+k];
00235 if (fabs(denom)<zero){
00236 fprintf (stderr,"\n\n zero denominator in expression for beta in function jacobi.\n");
00237 }
00238 beta=(0.0-a[j*n+k]-a[j*n+j]*alpha)/denom;
00239
00240
00241 for (l=0;l<n;l++){
00242 sol1=a[j*n+l];
00243 a[j*n+l]=a[j*n+l]+beta*a[k*n+l];
00244 a[k*n+l]=a[k*n+l]+alpha*sol1;
00245 sol1=b[j*n+l];
00246 b[j*n+l]=b[j*n+l]+beta*b[k*n+l];
00247 b[k*n+l]=b[k*n+l]+alpha*sol1;
00248 }
00249
00250
00251 for (l=0;l<n;l++){
00252 sol1=a[l*n+j];
00253 a[l*n+j]=a[l*n+j]+beta*a[l*n+k];
00254 a[l*n+k]=a[l*n+k]+alpha*sol1;
00255 sol1=b[l*n+j];
00256 b[l*n+j]=b[l*n+j]+beta*b[l*n+k];
00257 b[l*n+k]=b[l*n+k]+alpha*sol1;
00258 }
00259
00260
00261 for (l=0;l<n;l++){
00262 sol1=x[l*n+j];
00263 x[l*n+j]=x[l*n+j]+beta*x[l*n+k];
00264 x[l*n+k]=x[l*n+k]+alpha*sol1;
00265 }
00266
00267 }
00268 }
00269 }
00270 if (nze==numel){
00271 fprintf (stdout,"\n\n zmena zavory");
00272 nt++;
00273 if (nt==nthr) stop=1;
00274 else limit=thresholds[nt];
00275 }
00276
00277 if (stop==1) break;
00278 }
00279
00280
00281 for (i=0;i<n;i++){
00282 w[i]=a[i*n+i]/b[i*n+i];
00283 }
00284
00285
00286 for (i=0;i<n;i++){
00287 sol1=fabs(w[i]); k=i;
00288 for (j=i+1;j<n;j++){
00289 if (sol1>fabs(w[j])){
00290 sol1=fabs(w[j]); k=j;
00291 }
00292 }
00293 sol2=w[i];
00294 w[i]=w[k];
00295 w[k]=sol2;
00296 for (j=0;j<n;j++){
00297 sol2=x[j*n+i];
00298 x[j*n+i]=x[j*n+k];
00299 x[j*n+k]=sol2;
00300 }
00301 }
00302 }
00303
00304
00305
00306 void gen_jacobi2 (double *k,double *m,double *x,double *w,long n,long ni,
00307 double *gate,long ng,double limit)
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327 {
00328 long i,j,l,ii,jj,ll,rr,ss,zero,numel,nz;
00329 double kii,kij,kjj,mii,mij,mjj,s,k0,k1,k2,d,min,alpha,beta,err;
00330
00331 numel=(n*n-n)/2;
00332
00333
00334
00335
00336
00337 nullv (x,n*n);
00338
00339 for (i=0;i<n;i++){
00340 x[i*n+i]=1.0;
00341 }
00342
00343
00344
00345
00346 nz=0; err=gate[nz];
00347 for (l=0;l<ni;l++){
00348 zero=0;
00349 for (i=0;i<n;i++){
00350 ii=i*n+i; jj=ii; kii=k[ii]; mii=m[ii];
00351 for (j=i+1;j<n;j++){
00352 ii++; kij=k[ii]; mij=m[ii];
00353 jj+=n+1; kjj=k[jj]; mjj=m[jj];
00354 if (fabs(kii)<limit || fabs(kjj)<limit ||
00355 fabs(mii)<limit || fabs(mjj)<limit){
00356 if (fabs(kij)<err && fabs(mij)<err){
00357 zero++; continue;
00358 }
00359 }
00360 else{
00361 if (fabs(kij*kij/kii/kjj)<err && fabs(mij*mij/mii/mjj)<err){
00362 zero++; continue;
00363 }
00364 }
00365 k2=kij*mjj-kjj*mij;
00366 k1=kii*mjj-kjj*mii;
00367 k0=kii*mij-kij*mii;
00368 if (fabs(k2)<limit){
00369 if (fabs(k1)<limit){
00370 fprintf (stderr,"\n V Jacobiove metode je zdegenerovana kvadraticka rovnice");
00371 fprintf (stderr,"\n Vypocet konci. \n\n");
00372 }
00373 else{
00374 beta=-1.0*k0/k1;
00375 }
00376 }
00377 else{
00378 d=k1*k1-4.0*k2*k0;
00379 if (d<-1.0*limit){
00380 fprintf (stderr,"\n V Jacobiove metode je zaporny diskriminant.");
00381 fprintf (stderr,"\n Vypocet konci. \n\n");
00382 }
00383 s=-k1+sqrt(d);
00384 min=-k1-sqrt(d);
00385 if (fabs(s)>fabs(min)) beta=s/2.0/k2;
00386 else beta=min/2.0/k2;
00387
00388 }
00389 if (fabs(kii+beta*kij)<limit){
00390 beta=-1.0*kij/kjj;
00391 if (fabs(mii+beta*mij)<limit){
00392 fprintf (stderr,"\n Zmatky v Jacobiove metode rotaci.");
00393 fprintf (stderr,"\n Vypocet konci. \n\n");
00394 }
00395 else{
00396 alpha=(mij+beta*mjj)/(mii+beta*mij)*(-1.0);
00397 }
00398 }
00399 else{
00400 alpha=(kij+beta*kjj)/(kii+beta*kij)*(-1.0);
00401 }
00402
00403
00404
00405
00406 k[i*n+i]+=2.0*beta*kij+beta*beta*kjj;
00407 m[i*n+i]+=2.0*beta*mij+beta*beta*mjj;
00408 k[j*n+j]+=2.0*alpha*kij+alpha*alpha*kii;
00409 m[j*n+j]+=2.0*alpha*mij+alpha*alpha*mii;
00410 k[i*n+j]=0.0; m[i*n+j]=0.0;
00411 rr=i; ss=j;
00412 for (ll=0;ll<i;ll++){
00413 s=k[rr];
00414 k[rr]+=beta*k[ss];
00415 k[ss]+=alpha*s;
00416 s=m[rr];
00417 m[rr]+=beta*m[ss];
00418 m[ss]+=alpha*s;
00419 rr+=n; ss+=n;
00420 }
00421 rr=i*n+i+1; ss=(i+1)*n+j;
00422 for (ll=i+1;ll<j;ll++){
00423 s=k[rr];
00424 k[rr]+=beta*k[ss];
00425 k[ss]+=alpha*s;
00426 s=m[rr];
00427 m[rr]+=beta*m[ss];
00428 m[ss]+=alpha*s;
00429 rr++; ss+=n;
00430 }
00431 rr=i*n+j+1; ss=j*n+j+1;
00432 for (ll=j+1;ll<n;ll++){
00433 s=k[rr];
00434 k[rr]+=beta*k[ss];
00435 k[ss]+=alpha*s;
00436 s=m[rr];
00437 m[rr]+=beta*m[ss];
00438 m[ss]+=alpha*s;
00439 rr++; ss++;
00440 }
00441
00442
00443
00444
00445 rr=i; ss=j;
00446 for (ll=0;ll<n;ll++){
00447 s=x[rr];
00448 x[rr]+=beta*x[ss];
00449 x[ss]+=alpha*s;
00450 rr+=n; ss+=n;
00451 }
00452 }
00453 }
00454
00455 fprintf (stdout,"\n iterace %ld %ld",l,zero);
00456 if (zero==numel){
00457 nz++;
00458 if (nz<ng){
00459 err=gate[nz];
00460 fprintf (stdout,"\n Bude se pracovat s %ld. zavorou.",nz+1);
00461 fprintf (stdout,"\n Velikost zavory je %e",err);
00462 }
00463 else{
00464 break;
00465 }
00466 }
00467 }
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 for (i=0;i<n;i++){
00485 w[i]=k[i*n+i]/m[i*n+i];
00486 }
00487
00488
00489
00490
00491 for (i=0;i<n;i++){
00492 min=1.0e50;
00493 for (j=i;j<n;j++){
00494 if (w[j]<min){
00495 min=w[j]; ii=j;
00496 }
00497 }
00498 if (ii!=i){
00499 s=w[i];
00500 w[i]=w[ii];
00501 w[ii]=s;
00502 }
00503 }
00504 }