CHOLESKY.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
CHOLESKY.CPP

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ cholesky.cpp                     cholesky decomposition
00003 
00004 // Copyright (C) 1991,2,3,4: R B Davies
00005 
00006 #define WANT_MATH
00007 
00008 #include "include.h"
00009 
00010 #include "newmat.h"
00011 
00012 #ifdef use_namespace
00013 namespace NEWMAT {
00014 #endif
00015 
00016     /********* Cholesky decomposition of a positive definite matrix *************/
00017 
00018     // Suppose S is symmetrix and positive definite. Then there exists a unique
00019     // lower triangular matrix L such that L L.t() = S;
00020 
00021     inline Real square(Real x) { return x*x; }
00022 
00023     ReturnMatrix Cholesky(const SymmetricMatrix& S) {
00024         Tracer trace("Cholesky");
00025         int nr = S.Nrows();
00026         LowerTriangularMatrix T(nr);
00027         Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
00028         for (int i=0; i<nr; i++) {
00029             Real* tj = t; Real sum; int k;
00030             for (int j=0; j<i; j++) {
00031                 Real* tk = ti; sum = 0.0; k = j;
00032                 while (k--) { sum += *tj++ * *tk++; }
00033                 *tk = (*s++ - sum) / *tj++;
00034             }
00035             sum = 0.0; k = i;
00036             while (k--) { sum += square(*ti++); }
00037             Real d = *s++ - sum;
00038             if (d<=0.0) Throw(NPDException(S));
00039             *ti++ = sqrt(d);
00040         }
00041         T.Release(); return T.ForReturn();
00042     }
00043 
00044     ReturnMatrix Cholesky(const SymmetricBandMatrix& S) {
00045         Tracer trace("Band-Cholesky");
00046         int nr = S.Nrows(); int m = S.lower;
00047         LowerBandMatrix T(nr,m);
00048         Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
00049 
00050         for (int i=0; i<nr; i++) {
00051             Real* tj = t; Real sum; int l;
00052             if (i<m) { l = m-i; s += l; ti += l; l = i; }
00053             else { t += (m+1); l = m; }
00054 
00055             for (int j=0; j<l; j++) {
00056                 Real* tk = ti; sum = 0.0; int k = j; tj += (m-j);
00057                 while (k--) { sum += *tj++ * *tk++; }
00058                 *tk = (*s++ - sum) / *tj++;
00059             }
00060             sum = 0.0;
00061             while (l--) { sum += square(*ti++); }
00062             Real d = *s++ - sum;
00063             if (d<=0.0)  Throw(NPDException(S));
00064             *ti++ = sqrt(d);
00065         }
00066 
00067         T.Release(); return T.ForReturn();
00068     }
00069 
00070 #ifdef use_namespace
00071 }
00072 #endif
Generated on Fri Aug 23 01:37:01 2002 for VirtualU by
1.2.17