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;
170 template<
class mat_t,
class vec_t,
class vec2_t>
172 const vec_t &b, vec2_t &x) {
179 o2scl_cblas::o2cblas_Lower,
180 o2scl_cblas::o2cblas_NoTrans,
181 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
185 o2scl_cblas::o2cblas_Upper,
186 o2scl_cblas::o2cblas_NoTrans,
187 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
194 template<
class mat_t,
class vec_t>
199 o2scl_cblas::o2cblas_Lower,
200 o2scl_cblas::o2cblas_NoTrans,
201 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
205 o2scl_cblas::o2cblas_Upper,
206 o2scl_cblas::o2cblas_NoTrans,
207 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
228 O2SCL_IX2(LLT,j,j)=1.0/O2SCL_IX2(LLT,j,j);
229 double ajj=-O2SCL_IX2(LLT,j,j);
238 for (
size_t ii=N-j-1;ii>0 && ii--;) {
240 const size_t j_min=0;
241 const size_t j_max=ii;
243 for (
size_t jj=j_min;jj<j_max;jj++) {
244 temp+=O2SCL_IX2(LLT,jx+j+1,j)*O2SCL_IX2(LLT,ii+j+1,jj+j+1);
247 O2SCL_IX2(LLT,ix+j+1,j)=temp+O2SCL_IX2(LLT,ix+j+1,j)*
248 O2SCL_IX2(LLT,ii+j+1,ii+j+1);
269 for (j=i+1;j<N;++j) {
275 for(
size_t k=j;k<N;k++) {
276 sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,j);
280 O2SCL_IX2(LLT,i,j)=sum;
287 for(
size_t k=i;k<N;k++) {
288 sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,i);
291 O2SCL_IX2(LLT,i,i)=sum;
298 O2SCL_IX2(LLT,j,i)=O2SCL_IX2(LLT,i,j);
312 template<
class mat_t,
class vec_t>
322 const double C_ii=O2SCL_IX2(A,i,i);
323 O2SCL_IX(D,i)=C_ii*C_ii;
329 O2SCL_IX2(A,i,j)=O2SCL_IX2(A,i,j)/sqrt(O2SCL_IX(D,j));
340 O2SCL_IX2(A,i,j)=O2SCL_IX2(A,j,i);
void dtrsv(const enum o2cblas_order order, const enum o2cblas_uplo Uplo, const enum o2cblas_transpose TransA, const enum o2cblas_diag Diag, const size_t M, const size_t N, const mat_t &A, vec_t &X)
Compute .
invalid argument supplied by user
void cholesky_svx(const size_t N, const mat_t &LLT, vec_t &x)
Solve a linear system in place using a Cholesky decomposition.
int cholesky_decomp_unit(const size_t N, mat_t &A, vec_t &D)
Cholesky decomposition with unit-diagonal triangular parts.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
int cholesky_decomp(const size_t M, mat_t &A, bool err_on_fail=true)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix...
The namespace for linear algebra classes and functions.
double dnrm2_subrow(const mat_t &M, const size_t ir, const size_t ic, const size_t N)
Compute the norm of a subrow of a matrix.
void cholesky_solve(const size_t N, const mat_t &LLT, const vec_t &b, vec2_t &x)
Solve a symmetric positive-definite linear system after a Cholesky decomposition. ...
void cholesky_invert(const size_t N, mat_t &LLT)
Compute the inverse of a symmetric positive definite matrix given the Cholesky decomposition.
void dscal_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M, const double alpha)
Compute for a subcolumn of a matrix.