home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / jpinv < prev    next >
Encoding:
Text File  |  1995-11-15  |  1.2 KB  |  60 lines
  1. //---------------------------------------------------------------
  2. //
  3. // jpinv
  4. //
  5. // Syntax: Ap = jpinv(A)
  6. //
  7. // This routine computes the pseudo-inverse using Greville's
  8. // Iterative Method.
  9. //
  10. // A is the matrix to be inverted and the results is passed
  11. // back as Ap.
  12. //
  13. // Copyright (C), by Jeffrey B. Layton, 1994
  14. // Version JBL 940918
  15. //---------------------------------------------------------------
  16.  
  17. jpinv = function(A)
  18. {
  19.    local(nargs,b,mm,nn,k,Thresh,con,conn,dk,ck,bkT,Ap)
  20.  
  21. // Count number of input arguments
  22.    nargs=0;
  23.    if (exist(A)) {nargs=nargs+1;}
  24.  
  25.    if (nargs != 1) {
  26.        error("JPINV: Wrong number of input arguments.");
  27.    }
  28.  
  29.    Thresh = max(size(A))*epsilon()*norm(A,"2");
  30.    con = A[;1]'*A[;1];
  31.    conn = norm(con,"2");
  32.    Ap = zeros(A.nr,A.nc);
  33.  
  34.    if (conn >= Thresh) {
  35.       Ap=A[;1]'/[A[;1]'*A[;1]];
  36.    else
  37.       b=size(A[;1]');
  38.       mm=b[1];
  39.       nn=b[2];
  40.       Ap=zeros(mm,nn);
  41.    }
  42.  
  43.    if (A.nc != 1) {
  44.       for (k in 2:A.nc) {
  45.           dk=Ap*A[;k];
  46.           ck=A[;k]-A[;1:(k-1)]*dk;
  47.           con=ck'*ck;
  48.           conn=norm(con,"2");
  49.           if (conn >= Thresh) {
  50.              bkT=ck'/(ck'*ck);
  51.           else
  52.              bkT=dk'*Ap/(1+dk'*dk);
  53.           }
  54.           Ap=[Ap-dk*bkT;bkT];
  55.       }
  56.    }
  57.  
  58.    return Ap
  59. };
  60.