69 double A_00=O2SCL_IX2(A,0,0);
76 O2SCL_ERR2(
"Matrix not positive definite (A[0][0]<=0) in ",
80 double L_00=sqrt(A_00);
81 O2SCL_IX2(A,0,0)=L_00;
84 double A_10=O2SCL_IX2(A,1,0);
85 double A_11=O2SCL_IX2(A,1,1);
87 double L_10=A_10/L_00;
88 double diag=A_11-L_10*L_10;
91 O2SCL_ERR2(
"Matrix not positive definite (diag<=0 for 2x2) in ",
94 double L_11=sqrt(diag);
96 O2SCL_IX2(A,1,0)=L_10;
97 O2SCL_IX2(A,1,1)=L_11;
101 double A_kk=O2SCL_IX2(A,k,k);
106 double A_ki=O2SCL_IX2(A,k,i);
107 double A_ii=O2SCL_IX2(A,i,i);
113 sum+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,k,j);
117 A_ki=(A_ki-sum)/A_ii;
118 O2SCL_IX2(A,k,i)=A_ki;
123 double diag=A_kk-sum*sum;
126 O2SCL_ERR2(
"Matrix not positive definite (diag<=0) in ",
130 double L_kk=sqrt(diag);
132 O2SCL_IX2(A,k,k)=L_kk;
142 double A_ij=O2SCL_IX2(A,i,j);
143 O2SCL_IX2(A,j,i)=A_ij;
156 template<
class mat_t,
class vec_t,
class vec2_t>
158 const vec_t &b, vec2_t &x) {
165 o2scl_cblas::o2cblas_Lower,
166 o2scl_cblas::o2cblas_NoTrans,
167 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
171 o2scl_cblas::o2cblas_Upper,
172 o2scl_cblas::o2cblas_NoTrans,
173 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
180 template<
class mat_t,
class vec_t>
185 o2scl_cblas::o2cblas_Lower,
186 o2scl_cblas::o2cblas_NoTrans,
187 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
191 o2scl_cblas::o2cblas_Upper,
192 o2scl_cblas::o2cblas_NoTrans,
193 o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
214 O2SCL_IX2(LLT,j,j)=1.0/O2SCL_IX2(LLT,j,j);
215 double ajj=-O2SCL_IX2(LLT,j,j);
224 for (
size_t ii=N-j-1;ii>0 && ii--;) {
226 const size_t j_min=0;
227 const size_t j_max=ii;
229 for (
size_t jj=j_min;jj<j_max;jj++) {
230 temp+=O2SCL_IX2(LLT,jx+j+1,j)*O2SCL_IX2(LLT,ii+j+1,jj+j+1);
233 O2SCL_IX2(LLT,ix+j+1,j)=temp+O2SCL_IX2(LLT,ix+j+1,j)*
234 O2SCL_IX2(LLT,ii+j+1,ii+j+1);
255 for (j=i+1;j<N;++j) {
261 for(
size_t k=j;k<N;k++) {
262 sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,j);
266 O2SCL_IX2(LLT,i,j)=sum;
273 for(
size_t k=i;k<N;k++) {
274 sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,i);
277 O2SCL_IX2(LLT,i,i)=sum;
284 O2SCL_IX2(LLT,j,i)=O2SCL_IX2(LLT,i,j);
298 template<
class mat_t,
class vec_t>
308 const double C_ii=O2SCL_IX2(A,i,i);
309 O2SCL_IX(D,i)=C_ii*C_ii;
315 O2SCL_IX2(A,i,j)=O2SCL_IX2(A,i,j)/sqrt(O2SCL_IX(D,j));
326 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.
void cholesky_decomp(const size_t M, mat_t &A)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix...
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
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.