home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / eigsen < prev    next >
Encoding:
Text File  |  1995-11-15  |  4.1 KB  |  160 lines
  1. //---------------------------------------------------------------------------
  2. //
  3. // eigsen
  4. //
  5. // Syntax: D=eigsen(A,DADP)
  6. //
  7. // This routine computes the gradient of the eigenvalues of A and the
  8. // gradient of the eigenvectors of A. It assumes that A has distinct
  9. // eigenvalues. It uses Nelson's algorithm. The routine also normalizes the
  10. // eigenvectors and normalizes the gradient of the eigenvectors.
  11. //
  12. // If the routine is called as,
  13. //
  14. //    D=eigsen(A,DADP)
  15. //
  16. // where A is the square matrix of interest and DADP is the gradient of
  17. // A with respect to some scalar parameter p, it returns the gradient
  18. // of the eigenvalues (D.de) and the gradient of the eigenvalues (D.dX).
  19. //
  20. // The results are returned in a list:
  21. //
  22. //     D.de = Gradient of the eigenvalues
  23. //     D.dX = Gradient of the eigenvectors
  24. //
  25. // Note: The routine makes an effort to check if there are no repeated
  26. //       eigenvalues. However, only a warning is printed and the routine
  27. //       continues. It is up to the user to make sure there are no repeated
  28. //       eigenvalues.
  29. //
  30. // Ref: Nelson, R.B. "Simplified Calculation of Eigenvector Derivatives,"
  31. //      AIAA Journal, Vol. 14, No. 9, 1979, pp. 1201-1205.
  32. //
  33. // Copyright (C), by Jeffrey B. Layton, 1994
  34. // Version JBL 940918
  35. //---------------------------------------------------------------------------
  36.  
  37. rfile eignm
  38.  
  39. eigsen = function(A,DADP)
  40. {
  41.    local(nargs,i,k,Adum,AT,X,Y,D,con,F,Ap,V,C,de,dXdp,iloop,arg)
  42.  
  43. // Count number of inpute arguments
  44.    nargs=0;
  45.    if (exist(A)) {nargs=nargs+1;}
  46.    if (exist(DADP)) {nargs=nargs+1;}
  47.  
  48.    if (nargs < 2) {
  49.        error("EIGSEN: Wrong number of input arguments.");
  50.    }
  51.  
  52. // Initial Error Traps
  53.    if ( A.nr <= 1) {
  54.        error("EIGSEN: A must be a matrix.");
  55.        if (A.nr != A.nc) {
  56.            error("EIGSEN: A is not square.");
  57.        }
  58.        if (DADP.nr != DADP.nc) {
  59.            error("EIGSEN: DADP is not square.");
  60.        }
  61.        if ( (DADP.nr != A.nr) && (DADP.nc != A.nc) ) {
  62.            error("EIGSEN: A and DADP are not compatible.");
  63.        }
  64.    }
  65.  
  66. // Initialize the eigenvector gradient matrix
  67.    dXdp = zeros(A.nr,A.nc);
  68.  
  69. // Find eigenvalues and eigenvectors (both right and left) of A
  70.    Adum=eig(A);
  71.    D=conj(Adum.val');
  72.    X=Adum.vec;
  73.  
  74.    AT=A';
  75.    Adum=eig(AT);
  76.    Y=Adum.vec;
  77.  
  78. // Check for repeated eigenvalues
  79.    for (iloop in 1:(A.nr-1)) {
  80.         arg=A.val[iloop]-A.val[iloop+1];
  81.         if (real(arg) < (epsilon()*10) ) {
  82.             if (imag(arg) < (epsilon()*10) ) {
  83.                 error("eigsen: Repeated eigenvalue in A");
  84.             }
  85.         }
  86.    }
  87.  
  88. // Normalizing the right eigenvectors.
  89.    for (iloop in 1:A.nr) {
  90.        con=X[;iloop]'*X[;iloop];
  91.        if (con > 0.0) {
  92.            X[;iloop]=X[;iloop]/sqrt(con);
  93.        }
  94.    }
  95.  
  96. // Normalizing the left eigenvectors.
  97.    for (iloop in 1:A.nr) {
  98.        con=Y[;iloop]'*Y[;iloop];
  99.        if (con > 0.0) {
  100.            Y[;iloop]=Y[;iloop]/sqrt(con);
  101.        }
  102.    }
  103.  
  104. // Could just call eign to do the dirty work
  105.    Adum=eign(A,eye(A.nr,A.nc));
  106.    D=Adum.Lambda;
  107.    X=Adum.Phi;
  108.    Y=Adum.Psi;
  109.  
  110. // Initializing the right hand side matrix
  111.    F=zeros(A.nr,1);
  112.  
  113. // Initialize the eigenvalue gradient vector.
  114.    de=zeros(A.nr,1);
  115.  
  116. // Looping to find the particular solution in Nelson's
  117. // algorithm
  118.  
  119.    for (i in 1:A.nr) {
  120.  
  121. // Compute the gradient of the eigenvalues.
  122.        de[i]=Y[;i].'*DADP*X[;i];
  123.  
  124. // Finding the RHS of the modified problem of (A-lambda*I)
  125.        if (abs(D[i]) > epsilon() ) {
  126.            F=X[;i]*(de[i])-DADP*X[;i];
  127.  
  128. // Now Finding the largest value of an element of the eigenvector
  129.            abs(X[;i]).*abs(Y[;i])
  130.            k=maxi(abs(X[;i]).*abs(Y[;i]))
  131.  
  132. // Setting up the modified set of equations (A-Lambda*I)
  133.            Ap=A-(D[i]*eye(A.nr,A.nr));
  134.  
  135.            Ap[k;1:A.nr]=zeros(1,A.nr);
  136.            Ap[1:A.nr;k]=zeros(A.nr,1);
  137.  
  138.            Ap[k;k]=1;
  139.  
  140. // Setting up the RHS of the problem to find V.
  141.            F[k]=0;
  142.  
  143. // Solving for the vector V.
  144.            V=solve(Ap,F);
  145.            V[k]=0.0;
  146.  
  147. // Solving for the constant c.
  148.            C=-real((X[;i]'*V));
  149.  
  150. // Calculating dX/dp.
  151.            dXdp[;i]=V+C*X[;i];
  152.        else
  153.            D[i]
  154.            error("eigsen: Bug in Looping.");
  155.        }
  156.    }
  157.  
  158.    return << de=de; dX=dXdp >>
  159. };
  160.