56 for(
size_t i=0;i<N;i++) {
57 if (O2SCL_IX2(A,i,i)==0.0)
return 1;
94 for (j = 0; j < N - 1; j++) {
97 double ajj, max = fabs(O2SCL_IX2(A,j,j));
100 for (i = j + 1; i < N; i++) {
101 double aij = fabs (O2SCL_IX2(A,i,j));
114 temp=O2SCL_IX2(A,j,k);
115 O2SCL_IX2(A,j,k)=O2SCL_IX2(A,i_pivot,k);
116 O2SCL_IX2(A,i_pivot,k)=temp;
122 ajj = O2SCL_IX2(A,j,j);
125 for (i = j + 1; i < N; i++) {
126 double aij = O2SCL_IX2(A,i,j) / ajj;
127 O2SCL_IX2(A,i,j)=aij;
128 for (k = j + 1; k < N; k++) {
129 double aik = O2SCL_IX2(A,i,k);
130 double ajk = O2SCL_IX2(A,j,k);
131 O2SCL_IX2(A,i,k)=aik - aij * ajk;
140 #ifdef O2SCL_NEVER_DEFINED
142 template<
class mat_t,
class vec_
size_t>
143 int LU_decomp_L3_sub(
const size_t M,
const size_t N, mat_t &A,
144 vec_size_t &ipiv,
size_t istart,
149 }
else if (N <= 24) {
151 return LU_decomp_L2(A, ipiv);
165 const size_t N1=((N >= 16) ? ((N + 8) / 16) * 8 : N / 2);
166 const size_t N2 = N - N1;
167 const size_t M2 = M - N1;
188 status=LU_decomp_L3(M,N1,AL,ipiv,0,0);
194 apply_pivots(&AR.matrix, &ipiv1.vector);
197 dtrsm_submat(CblasLeft,CblasLower,CblasNoTrans,CblasUnit,
203 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, &A21.matrix, &A12.matrix, 1.0, &A22.matrix);
206 status = LU_decomp_L3(&A22.matrix, &ipiv2.vector);
211 apply_pivots(&A21.matrix, &ipiv2.vector);
214 for (i = 0; i < N2; ++i)
216 unsigned int * ptr = gsl_vector_uint_ptr(&ipiv2.vector, i);
234 template<
class mat_t,
class vec_t>
235 int LU_svx(
const size_t N,
const mat_t &LU,
239 O2SCL_ERR(
"LU matrix is singular in LU_svx().",
248 o2scl_cblas::o2cblas_Lower,
249 o2scl_cblas::o2cblas_NoTrans,
250 o2scl_cblas::o2cblas_Unit,
255 o2scl_cblas::o2cblas_Upper,
256 o2scl_cblas::o2cblas_NoTrans,
257 o2scl_cblas::o2cblas_NonUnit,
271 template<
class mat_t,
class mat_row_t>
280 for (j = 0; j < N - 1; j++) {
283 double ajj, max = fabs(O2SCL_IX2(A,j,j));
286 for (i = j + 1; i < N; i++) {
287 double aij = fabs (O2SCL_IX2(A,i,j));
300 mat_row_t r2(A,i_pivot);
303 O2SCL_IX(r1,k)=O2SCL_IX(r2,k);
310 ajj = O2SCL_IX2(A,j,j);
313 for (i = j + 1; i < N; i++) {
314 double aij = O2SCL_IX2(A,i,j) / ajj;
315 O2SCL_IX2(A,i,j)=aij;
316 for (k = j + 1; k < N; k++) {
317 double aik = O2SCL_IX2(A,i,k);
318 double ajk = O2SCL_IX2(A,j,k);
319 O2SCL_IX2(A,i,k)=aik - aij * ajk;
334 template<
class mat_t,
class vec_t,
class vec2_t>
336 const vec_t &b, vec2_t &x) {
339 O2SCL_ERR(
"LU matrix is singular in LU_solve().",
357 template<
class mat_t,
class mat2_t,
class vec_t,
class vec2_t,
class vec3_t>
358 int LU_refine(
const size_t N,
const mat_t &A,
const mat2_t &LU,
363 O2SCL_ERR(
"LU matrix is singular in LU_refine().",
370 o2scl_cblas::o2cblas_NoTrans,
371 N,N,1.0,A,x,-1.0,residual);
375 int status=
LU_svx(N,LU,p,residual);
395 template<
class mat_t,
class mat2_t,
class mat_col_t>
400 O2SCL_ERR(
"LU matrix is singular in LU_invert().",
410 for(
size_t j=0;j<N;j++) {
411 if (i==j) O2SCL_IX2(inverse,i,j)=1.0;
412 else O2SCL_IX2(inverse,i,j)=0.0;
416 for (i = 0; i < N; i++) {
417 mat_col_t c(inverse,i);
418 int status_i=
LU_svx(N,LU,p,c);
434 template<
class mat_t>
435 double LU_det(
const size_t N,
const mat_t &LU,
int signum) {
439 double det=(double)signum;
442 det*=O2SCL_IX2(LU,i,i);
456 template<
class mat_t>
462 for (i = 0; i < N; i++) {
463 lndet+=log(fabs(O2SCL_IX2(LU,i,i)));
475 template<
class mat_t>
476 int LU_sgndet(
const size_t N,
const mat_t &LU,
int signum) {
482 double u=O2SCL_IX2(LU,i,i);