home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / eigtest < prev    next >
Encoding:
Text File  |  1995-11-15  |  1.8 KB  |  83 lines
  1.  
  2. // This routine test the eigen sensitivity routines.
  3.  
  4. rfile eigsen
  5. rfile eigsen2
  6. rfile banner
  7. rfile eignm
  8.  
  9. banner("Eigen Sensitivity Test.");
  10. // Generate a random 4x4 matrix for A
  11. A=5.*rand(4,4);
  12. A=real(A);
  13. //A=(A+A')/2;
  14.  
  15. // Generate a random 4x4 matrix for d(A)/dp matrix
  16. DADP=2.*rand(4,4);
  17.  
  18. // Set B=I for eigsen (Junkins method)
  19. B=eye(4,4);
  20.  
  21. // Set d(B)/dp = 0
  22. DBDP = zeros(4,4);
  23.  
  24. // Compute eigenvalues, right and left eigenvectors of A
  25. E=eig(A);
  26. X=E.vec;
  27. Lambda=E.val;
  28. ET=eig(A');
  29. Y=ET.vec;
  30.  
  31. // Clear unnecessary variables
  32. clear(E,ET);
  33.  
  34. Adum=eignm(A,B);
  35. X=Adum.Phi;
  36. Y=Adum.Psi;
  37. Lambda=Adum.Lambda;
  38. // Compute Eigen sensitivity using Junkins' Method.
  39. Es=eigsen(A,B,DADP,DBDP);
  40. //clear(Adum);
  41.  
  42. // Compute Eigen sensitivity using Nelson's Method.
  43. Es2=eigsen2(A,DADP);
  44.  
  45. // Compute norms of the errors
  46. REnorm=norm(Es.REVTD-Es2.dX);
  47.  
  48. Lnorm=norm(Es.EVAD-Es2.de);
  49.  
  50. // Clean-up
  51. clear(DADP,B,DBDP,X,Y,Lambda);
  52.  
  53. if (!(REnorm < 1.0e-12) ) {
  54.     printf("%s","FAILED Right Eigen Sensitivity test.\n");
  55.     printf("%s","The norm of the Right Eigenvector Sensitivity Error is");
  56.     sprintf(sdum,"%12.5e",REnorm);
  57.     s="   "+sdum;
  58.     printf("%s\n",s);
  59. else
  60.     printf("%s","Passed Right Eigenvector Sensitivity Check.\n");
  61.     printf("%s","The norm of the Right Eigenvector Sensitivity Error is");
  62.     sprintf(sdum,"%12.5e",REnorm);
  63.     s="   "+sdum;
  64.     printf("%s\n",s);
  65. }
  66.  
  67. printf("%s"," \n");
  68. if (!(Lnorm < 1.0e-12) ) {
  69.     printf("%s","Failed Eigenvalue sensitivity test.\n");
  70.     printf("%s","The norm of the Eigenvalue Sensitivity Error is");
  71.     sprintf(sdum,"%12.5e",Lnorm);
  72.     s="   "+sdum;
  73.     printf("%s\n",s);
  74. else
  75.     printf("%s","Passed Eigenvalue Sensitivity Test.\n");
  76.     printf("%s","The norm of the Eigenvalue Sensitivity Error is");
  77.     sprintf(sdum,"%12.5e",Lnorm);
  78.     s="   "+sdum;
  79.     printf("%s\n",s);
  80. }
  81.  
  82. clear(REnorm,Lnorm,s,sdum);
  83.