home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2398 / dsweep.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1990-12-28  |  2.6 KB  |  115 lines
  1. /*
  2. Name          dsweep - Inverts a positive definite matrix in place by a 
  3.                        diagonal sweep.
  4.  
  5. Usage         #include "usual.h"
  6.               #include "tools.h"
  7.               INTEGER dsweep(REAL *a, INTEGER n, REAL eps)
  8.  
  9. Prototype in  tools.h
  10.  
  11. Description   a is a symmetric, positive definite, n by n matrix stored 
  12.               columnwise with no unused space and first element in a[0]; i.e.,
  13.               for (j=1; j<=n; j++) for (i=1; i<=n; i++) aij=a[n*(j-1)+(i-1)];
  14.               will traverse the matrix with aij being the element in the 
  15.               i-th row and j-th column.  On return a contains the inverse of a
  16.               stored columnwise.  eps is an input constant used as a relative 
  17.               tolerance in testing for degenerate rank.  A reasonable value 
  18.               for eps is 1.e-13. 
  19.  
  20. Remark        dsweep.cc is dsweep.f translated to C++.
  21.  
  22. Reference     Schatzoff, M. et al.  Efficient Calculation of all Possible 
  23.               Regressions. Technometrics, 10. 769-79 (November 1968)
  24.  
  25. Return value  The return value ier is an error parameter coded as follows:
  26.               ier=0, no error, rank of a is n;  ier > 0, a is singular, rank 
  27.               of a is n-ier.  If ier > 0 then dsweep returns a g-inverse. 
  28.  
  29. Functions     Library: (none)
  30. called        Sublib:  (none)
  31. */
  32.  
  33. #include "usual.h"
  34. #include "tools.h"
  35.  
  36.  
  37. INTEGER dsweep(REAL *a, INTEGER n, REAL eps)
  38. {
  39.   INTEGER i, j, k, kk, ier;
  40.   REAL    test, d, aik, akj, akk, tol;
  41.  
  42.   a--;   // Adjust offset to allow Fortran style indexing.
  43.  
  44.   tol=0.0;
  45.  
  46.   for (i=1; i<=n; i++) {
  47.     test=a[n*(i-1)+i];
  48.     if(test > tol) tol=test;
  49.   }
  50.  
  51.   tol=tol*eps;
  52.   ier=0;
  53.  
  54.   for (k=1; k<=n; k++) {
  55.  
  56.     kk=n*(k-1)+k;
  57.  
  58.     akk=a[kk];
  59.  
  60.     if (akk < tol) {
  61.  
  62.       ier=ier+1;
  63.  
  64.       for (j=k; j<=n; j++)     // zero out row k of upper triangle
  65.         a[n*(j-1)+k]=0.0;
  66.  
  67.       for (i=1; i<=k-1; i++)  // zero out column k of upper triangle
  68.         a[kk-i]=0.0;
  69.  
  70.     } 
  71.     else {  // sweep
  72.  
  73.       d=1.0/akk;
  74.  
  75.       for (i=1; i <= n; i++) {
  76.       for (j=i; j <= n; j++) {
  77.         if (i != k && j != k)  {
  78.  
  79.           if (i<k) 
  80.             aik=a[n*(k-1)+i];
  81.           else 
  82.             aik=a[n*(i-1)+k];
  83.  
  84.           if (k<j) 
  85.             akj=a[n*(j-1)+k];
  86.           else 
  87.             akj=-a[n*(k-1)+j];
  88.  
  89.           a[n*(j-1)+i] -= aik*akj*d;
  90.         }
  91.       }
  92.       }
  93.  
  94.       for (j=k; j<=n; j++)
  95.         a[n*(j-1)+k] *= d;
  96.  
  97.       for (i=1; i<=k-1; i++)
  98.         a[kk-i] = -a[kk-i]*d;
  99.  
  100.       a[kk]=d;
  101.  
  102.     }
  103.   }
  104.  
  105.   // copy the upper triangle to the lower
  106.  
  107.   for (i=1; i <= n; i++)
  108.   for (j=i; j <= n; j++)
  109.     a[n*(i-1)+j]=a[n*(j-1)+i];
  110.  
  111.   return ier;
  112. }
  113.  
  114.  
  115.