home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume34 / newmat06 / part06 / cholesky.cxx < prev    next >
C/C++ Source or Header  |  1992-12-06  |  2KB  |  72 lines
  1. //$$ cholesky.cxx                     cholesky decomposition
  2.  
  3. // Copyright (C) 1991,2: R B Davies
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.h"
  8.  
  9. #include "newmat.h"
  10.  
  11.  
  12. /********* Cholesky decomposition of a positive definite matrix *************/
  13.  
  14. // Suppose S is symmetrix and positive definite. Then there exists a unique
  15. // lower triangular matrix L such that L L.t() = S;
  16.  
  17. static Real square(Real x) { return x*x; }
  18.  
  19. ReturnMatrix Cholesky(const SymmetricMatrix& S)
  20. {
  21.    Tracer trace("Cholesky");
  22.    int nr = S.Nrows();
  23.    LowerTriangularMatrix T(nr);
  24.    Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
  25.    for (int i=0; i<nr; i++)
  26.    {
  27.       Real* tj = t; Real sum; int k;
  28.       for (int j=0; j<i; j++)
  29.       {
  30.      Real* tk = ti; sum = 0.0; k = j;
  31.      while (k--) { sum += *tj++ * *tk++; }
  32.      *tk = (*s++ - sum) / *tj++;
  33.       }
  34.       sum = 0.0; k = i;
  35.       while (k--) { sum += square(*ti++); }
  36.       Real d = *s++ - sum;
  37.       if (d<=0.0) Throw(NPDException(S));
  38.       *ti++ = sqrt(d);
  39.    }
  40.    T.Release(); return (ReturnMatrix)T;
  41. }
  42.  
  43. ReturnMatrix Cholesky(const SymmetricBandMatrix& S)
  44. {
  45.    Tracer trace("Band-Cholesky");
  46.    int nr = S.Nrows(); int m = S.lower;
  47.    LowerBandMatrix T(nr,m);
  48.    Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
  49.  
  50.    for (int i=0; i<nr; i++)
  51.    {
  52.       Real* tj = t; Real sum; int l;
  53.       if (i<m) { l = m-i; s += l; ti += l; l = i; }
  54.       else { t += (m+1); l = m; }
  55.  
  56.       for (int j=0; j<l; j++)
  57.       {
  58.      Real* tk = ti; sum = 0.0; int k = j; tj += (m-j);
  59.      while (k--) { sum += *tj++ * *tk++; }
  60.      *tk = (*s++ - sum) / *tj++;
  61.       }
  62.       sum = 0.0;
  63.       while (l--) { sum += square(*ti++); }
  64.       Real d = *s++ - sum;
  65.       if (d<=0.0)  Throw(NPDException(S));
  66.       *ti++ = sqrt(d);
  67.    }
  68.  
  69.    T.Release(); return (ReturnMatrix)T;
  70. }
  71.  
  72.