59 void resize(
size_t n) {
73 ubvector v1, v2, v3, v4;
74 void resize(
size_t n) {
89 ubvector v1, v2, v3, v4, v5;
90 void resize(
size_t n) {
116 template<
class vec_t,
class vec2_t,
class vec3_t,
117 class vec4_t,
class mem_t,
class mem_vec_t>
119 const vec3_t &b, vec4_t &x,
size_t N, mem_t &m) {
121 mem_vec_t &gamma=m.v1;
122 mem_vec_t &alpha=m.v2;
133 alpha[0] = O2SCL_IX(diag,0);
134 gamma[0] = O2SCL_IX(offdiag,0) / alpha[0];
136 for (i = 1; i < N - 1; i++) {
137 alpha[i] = O2SCL_IX(diag,i) - O2SCL_IX(offdiag,i-1) * gamma[i - 1];
138 gamma[i] = O2SCL_IX(offdiag,i) / alpha[i];
142 alpha[N-1]=O2SCL_IX(diag,N-1)-O2SCL_IX(offdiag,N-2)*gamma[N-2];
147 for (i = 1; i < N; i++) {
148 z[i] = O2SCL_IX(b,i) - gamma[i - 1] * z[i - 1];
149 }
for (i = 0; i < N; i++) {
150 c[i] = z[i] / alpha[i];
154 O2SCL_IX(x,N-1)=c[N-1];
156 for (i = N - 2, j = 0; j <= N - 2; j++, i--) {
157 O2SCL_IX(x,i) = c[i] - gamma[i] * O2SCL_IX(x,i+1);
181 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t,
182 class vec5_t,
class mem_t,
class mem_vec_t>
184 const vec3_t &belowdiag,
const vec4_t &rhs,
185 vec5_t &x,
size_t N, mem_t &m) {
186 mem_vec_t &alpha=m.v1;
196 alpha[0] = O2SCL_IX(diag,0);
197 z[0] = O2SCL_IX(rhs,0);
199 for (i = 1; i < N; i++) {
200 const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1];
201 alpha[i] = O2SCL_IX(diag,i) - t*O2SCL_IX(abovediag,i-1);
202 z[i] = O2SCL_IX(rhs,i) - t*z[i-1];
211 O2SCL_IX(x,N-1) = z[N - 1]/alpha[N - 1];
213 for (i = N - 2, j = 0; j <= N - 2; j++, i--) {
214 O2SCL_IX(x,i) = (z[i] - O2SCL_IX(abovediag,i) *
215 O2SCL_IX(x,i+1))/alpha[i];
239 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t,
240 class mem_t,
class mem_vec_t>
242 const vec3_t &b, vec4_t &x,
size_t N,
245 mem_vec_t &delta=m.v1;
246 mem_vec_t &gamma=m.v2;
247 mem_vec_t &alpha=m.v3;
257 x[0] = b[0] / O2SCL_IX(diag,0);
261 alpha[0] = O2SCL_IX(diag,0);
262 gamma[0] = O2SCL_IX(offdiag,0) / alpha[0];
263 delta[0] = O2SCL_IX(offdiag,N-1) / alpha[0];
265 for (i = 1; i < N - 2; i++) {
266 alpha[i] = O2SCL_IX(diag,i) - O2SCL_IX(offdiag,i-1) * gamma[i - 1];
267 gamma[i] = O2SCL_IX(offdiag,i) / alpha[i];
268 delta[i] = -delta[i - 1] * O2SCL_IX(offdiag,i-1) / alpha[i];
271 for (i = 0; i < N - 2; i++) {
272 sum += alpha[i] * delta[i] * delta[i];
275 alpha[N - 2] = diag[ (N - 2)] - O2SCL_IX(offdiag,N-3) * gamma[N - 3];
277 gamma[N - 2] = (offdiag[(N - 2)] - offdiag[(N - 3)] *
278 delta[N - 3]) / alpha[N - 2];
280 alpha[N - 1] = diag[(N - 1)] - sum - alpha[(N - 2)] *
281 gamma[N - 2] * gamma[N - 2];
286 for (i = 1; i < N - 1; i++) {
287 z[i] = O2SCL_IX(b,i) - z[i - 1] * gamma[i - 1];
290 for (i = 0; i < N - 2; i++) {
291 sum += delta[i] * z[i];
293 z[N - 1] = b[(N - 1)] - sum - gamma[N - 2] * z[N - 2];
294 for (i = 0; i < N; i++) {
295 c[i] = z[i] / alpha[i];
299 O2SCL_IX(x,N-1) = c[N - 1];
300 x[(N - 2)] = c[N - 2] - gamma[N - 2] * O2SCL_IX(x,N-1);
302 for (i = N - 3, j = 0; j <= N - 3; j++, i--) {
303 O2SCL_IX(x,i) = c[i] - gamma[i] * x[(i + 1)] -
304 delta[i] * O2SCL_IX(x,N-1);
336 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t,
337 class vec5_t,
class mem_t,
class mem_vec_t>
339 const vec3_t &belowdiag,
const vec4_t &rhs,
340 vec5_t &x,
size_t N, mem_t &m) {
343 mem_vec_t &alpha=m.v1;
354 zb[0] = O2SCL_IX(rhs,0);
355 if (O2SCL_IX(diag,0) != 0) {
356 beta = -O2SCL_IX(diag,0);
360 const double q = 1 - O2SCL_IX(abovediag,0)*O2SCL_IX(belowdiag,0)/
361 (O2SCL_IX(diag,0)*diag[1]);
362 if (fabs(q/beta) > 0.5 && fabs(q/beta) < 2) {
363 beta *= (fabs(q/beta) < 1) ? 0.5 : 2;
366 alpha[0] = O2SCL_IX(diag,0) - beta;
370 for (i = 1; i+1 < N; i++) {
371 const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1];
372 alpha[i] = O2SCL_IX(diag,i) - t*O2SCL_IX(abovediag,i-1);
373 zb[i] = O2SCL_IX(rhs,i) - t*zb[i-1];
384 const size_t i = N-1;
385 const double t = O2SCL_IX(belowdiag,i-1)/alpha[i-1];
386 alpha[i]=O2SCL_IX(diag,i)-O2SCL_IX(abovediag,i)*
387 O2SCL_IX(belowdiag,i)/beta-
388 t*O2SCL_IX(abovediag,i-1);
389 zb[i] = O2SCL_IX(rhs,i) - t*zb[i-1];
390 zu[i] = O2SCL_IX(abovediag,i) - t*zu[i-1];
402 w[N-1] = zu[N-1]/alpha[N-1];
403 O2SCL_IX(x,N-1) = zb[N-1]/alpha[N-1];
404 for (i = N - 2, j = 0; j <= N - 2; j++, i--) {
405 w[i] = (zu[i] - O2SCL_IX(abovediag,i) * w[i+1])/alpha[i];
406 O2SCL_IX(x,i) = (zb[i] - O2SCL_IX(abovediag,i) *
407 O2SCL_IX(x,i + 1))/alpha[i];
412 const double vw = w[0] + O2SCL_IX(belowdiag,N-1)/beta * w[N-1];
413 const double vx = O2SCL_IX(x,0) +
414 O2SCL_IX(belowdiag,N-1)/beta * O2SCL_IX(x,N-1);
424 for (i = 0; i < N; i++)
425 O2SCL_IX(x,i) -= vx/(1 + vw)*w[i];
433 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t>
435 const vec3_t &b, vec4_t &x,
size_t N) {
440 ubvector>(diag,offdiag,b,x,N,v4m);
446 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t,
449 const vec3_t &belowdiag,
const vec4_t &rhs,
450 vec5_t &x,
size_t N) {
455 ubvector>(diag,abovediag,belowdiag,rhs,x,N,v2m);
461 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t>
463 const vec3_t &b, vec4_t &x,
size_t N) {
468 ubvector>(diag,offdiag,b,x,N,v5m);
474 template<
class vec_t,
class vec2_t,
class vec3_t,
class vec4_t,
477 const vec3_t &belowdiag,
const vec4_t &rhs,
478 vec5_t &x,
size_t N) {
void solve_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
Solve a symmetric tridiagonal linear system with user-specified memory.
Allocation object for 5 arrays of equal size.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
The namespace for linear algebra classes and functions.
Allocation object for 4 arrays of equal size.
void solve_cyc_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag, const vec3_t &belowdiag, const vec4_t &rhs, vec5_t &x, size_t N, mem_t &m)
Solve an asymmetric cyclic tridiagonal linear system with user-specified memory.
void solve_tridiag_nonsym(const vec_t &diag, const vec2_t &abovediag, const vec3_t &belowdiag, const vec4_t &rhs, vec5_t &x, size_t N, mem_t &m)
Solve an asymmetric tridiagonal linear system with user-specified memory.
Allocation object for 2 arrays of equal size.
void solve_cyc_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
Solve a symmetric cyclic tridiagonal linear system with user specified memory.