63 bool err_on_fail=
true) {
71 double A_00=O2SCL_IX2(A,0,0);
79 O2SCL_ERR2(
"Matrix not positive definite (A[0][0]<=0) in ",
86 double L_00=sqrt(A_00);
87 O2SCL_IX2(A,0,0)=L_00;
90 double A_10=O2SCL_IX2(A,1,0);
91 double A_11=O2SCL_IX2(A,1,1);
93 double L_10=A_10/L_00;
94 double diag=A_11-L_10*L_10;
98 O2SCL_ERR2(
"Matrix not positive definite (diag<=0 for 2x2) in ",
104 double L_11=sqrt(diag);
106 O2SCL_IX2(A,1,0)=L_10;
107 O2SCL_IX2(A,1,1)=L_11;
111 double A_kk=O2SCL_IX2(A,k,k);
116 double A_ki=O2SCL_IX2(A,k,i);
117 double A_ii=O2SCL_IX2(A,i,i);
123 sum+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,k,j);
127 A_ki=(A_ki-sum)/A_ii;
128 O2SCL_IX2(A,k,i)=A_ki;
133 double diag=A_kk-sum*sum;
137 O2SCL_ERR2(
"Matrix not positive definite (diag<=0) in ",
144 double L_kk=sqrt(diag);
146 O2SCL_IX2(A,k,k)=L_kk;
156 double A_ij=O2SCL_IX2(A,i,j);
157 O2SCL_IX2(A,j,i)=A_ij;
167 template<
class mat_t>
172 for (
size_t i=0;i<M;i++) {
173 det*=O2SCL_IX2(A,i,i);
182 template<
class mat_t>
187 for (
size_t i = 0; i < M; i++) {
188 lndet+=log(fabs(O2SCL_IX2(A,i,i)));
200 template<
class mat_t,
class vec_t,
class vec2_t>
202 const vec_t &b, vec2_t &x) {
209 o2scl_cblas::o2cblas_Lower,
210 o2scl_cblas::o2cblas_NoTrans,
211 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
215 o2scl_cblas::o2cblas_Upper,
216 o2scl_cblas::o2cblas_NoTrans,
217 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
224 template<
class mat_t,
class vec_t>
229 o2scl_cblas::o2cblas_Lower,
230 o2scl_cblas::o2cblas_NoTrans,
231 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
235 o2scl_cblas::o2cblas_Upper,
236 o2scl_cblas::o2cblas_NoTrans,
237 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
258 O2SCL_IX2(LLT,j,j)=1.0/O2SCL_IX2(LLT,j,j);
259 double ajj=-O2SCL_IX2(LLT,j,j);
268 for (
size_t ii=N-j-1;ii>0 && ii--;) {
270 const size_t j_min=0;
271 const size_t j_max=ii;
273 for (
size_t jj=j_min;jj<j_max;jj++) {
274 temp+=O2SCL_IX2(LLT,jx+j+1,j)*O2SCL_IX2(LLT,ii+j+1,jj+j+1);
277 O2SCL_IX2(LLT,ix+j+1,j)=temp+O2SCL_IX2(LLT,ix+j+1,j)*
278 O2SCL_IX2(LLT,ii+j+1,ii+j+1);
299 for (j=i+1;j<N;++j) {
305 for(
size_t k=j;k<N;k++) {
306 sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,j);
310 O2SCL_IX2(LLT,i,j)=sum;
317 for(
size_t k=i;k<N;k++) {
318 sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,i);
321 O2SCL_IX2(LLT,i,i)=sum;
328 O2SCL_IX2(LLT,j,i)=O2SCL_IX2(LLT,i,j);
342 template<
class mat_t,
class vec_t>
352 const double C_ii=O2SCL_IX2(A,i,i);
353 O2SCL_IX(D,i)=C_ii*C_ii;
359 O2SCL_IX2(A,i,j)=O2SCL_IX2(A,i,j)/sqrt(O2SCL_IX(D,j));
370 O2SCL_IX2(A,i,j)=O2SCL_IX2(A,j,i);