00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef CHOLESKY_H
00027 #define CHOLESKY_H
00028
00029 #include <cmath>
00030
00031
00032
00033 namespace SPUC
00034 {
00035
00036
00037
00038
00039
00040
00041 template <class SPDMatrix, class SymmMatrix>
00042 int Cholesky_upper_factorization(SPDMatrix &A, SymmMatrix &L)
00043 {
00044 int M = A.dim(1);
00045 int N = A.dim(2);
00046
00047 assert(M == N);
00048
00049
00050
00051 if (M != L.dim(1) || N != L.dim(2))
00052 L = SymmMatrix(N,N);
00053
00054 int i,j,k;
00055
00056
00057 typename SPDMatrix::element_type dot=0;
00058
00059
00060 for (j=1; j<=N; j++)
00061 {
00062 dot= 0;
00063
00064 for (i=1; i<j; i++)
00065 dot = dot + L(j,i)*L(j,i);
00066
00067 L(j,j) = A(j,j) - dot;
00068
00069 for (i=j+1; i<=N; i++)
00070 {
00071 dot = 0;
00072 for (k=1; k<j; k++)
00073 dot = dot + L(i,k)*L(j,k);
00074 L(i,j) = A(j,i) - dot;
00075 }
00076
00077 if (L(j,j) <= 0.0) return 1;
00078
00079 L(j,j) = sqrt( L(j,j) );
00080
00081 for (i=j+1; i<=N; i++)
00082 L(i,j) = L(i,j) / L(j,j);
00083
00084 }
00085
00086 return 0;
00087 }
00088
00089
00090
00091
00092 }
00093
00094
00095 #endif
00096