Blender V2.61 - r43446

cholesky.h

Go to the documentation of this file.
00001 /*
00002  */
00003 
00004 /*
00005 *
00006 * Template Numerical Toolkit (TNT): Linear Algebra Module
00007 *
00008 * Mathematical and Computational Sciences Division
00009 * National Institute of Technology,
00010 * Gaithersburg, MD USA
00011 *
00012 *
00013 * This software was developed at the National Institute of Standards and
00014 * Technology (NIST) by employees of the Federal Government in the course
00015 * of their official duties. Pursuant to title 17 Section 105 of the
00016 * United States Code, this software is not subject to copyright protection
00017 * and is in the public domain.  The Template Numerical Toolkit (TNT) is
00018 * an experimental system.  NIST assumes no responsibility whatsoever for
00019 * its use by other parties, and makes no guarantees, expressed or implied,
00020 * about its quality, reliability, or any other characteristic.
00021 *
00022 * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
00023 * see http://math.nist.gov/tnt for latest updates.
00024 *
00025 */
00026 
00027 #ifndef CHOLESKY_H
00028 #define CHOLESKY_H
00029 
00030 #include <cmath>
00031 
00032 // index method
00033 
00034 namespace TNT
00035 {
00036 
00037 
00038 //
00039 // Only upper part of A is used.  Cholesky factor is returned in
00040 // lower part of L.  Returns 0 if successful, 1 otherwise.
00041 //
00042 template <class SPDMatrix, class SymmMatrix>
00043 int Cholesky_upper_factorization(SPDMatrix &A, SymmMatrix &L)
00044 {
00045     Subscript M = A.dim(1);
00046     Subscript N = A.dim(2);
00047 
00048     assert(M == N);                 // make sure A is square
00049 
00050     // readjust size of L, if necessary
00051 
00052     if (M != L.dim(1) || N != L.dim(2))
00053         L = SymmMatrix(N,N);
00054 
00055     Subscript i,j,k;
00056 
00057 
00058     typename SPDMatrix::element_type dot=0;
00059 
00060 
00061     for (j=1; j<=N; j++)                // form column j of L
00062     {
00063         dot= 0;
00064 
00065         for (i=1; i<j; i++)             // for k= 1 TO j-1
00066             dot = dot +  L(j,i)*L(j,i);
00067 
00068         L(j,j) = A(j,j) - dot;
00069 
00070         for (i=j+1; i<=N; i++)
00071         {
00072             dot = 0;
00073             for (k=1; k<j; k++)
00074                 dot = dot +  L(i,k)*L(j,k);
00075             L(i,j) = A(j,i) - dot;
00076         }
00077 
00078         if (L(j,j) <= 0.0) return 1;
00079 
00080         L(j,j) = sqrt( L(j,j) );
00081 
00082         for (i=j+1; i<=N; i++)
00083             L(i,j) = L(i,j) / L(j,j);
00084 
00085     }
00086 
00087     return 0;
00088 }
00089 
00090 
00091 
00092 
00093 }  
00094 // namespace TNT
00095 
00096 #endif
00097 // CHOLESKY_H
00098