1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
|
#include "rb_lapack.h"
extern VOID dlagts_(integer* job, integer* n, doublereal* a, doublereal* b, doublereal* c, doublereal* d, integer* in, doublereal* y, doublereal* tol, integer* info);
static VALUE
rblapack_dlagts(int argc, VALUE *argv, VALUE self){
VALUE rblapack_job;
integer job;
VALUE rblapack_a;
doublereal *a;
VALUE rblapack_b;
doublereal *b;
VALUE rblapack_c;
doublereal *c;
VALUE rblapack_d;
doublereal *d;
VALUE rblapack_in;
integer *in;
VALUE rblapack_y;
doublereal *y;
VALUE rblapack_tol;
doublereal tol;
VALUE rblapack_info;
integer info;
VALUE rblapack_y_out__;
doublereal *y_out__;
integer n;
VALUE rblapack_options;
if (argc > 0 && TYPE(argv[argc-1]) == T_HASH) {
argc--;
rblapack_options = argv[argc];
if (rb_hash_aref(rblapack_options, sHelp) == Qtrue) {
printf("%s\n", "USAGE:\n info, y, tol = NumRu::Lapack.dlagts( job, a, b, c, d, in, y, tol, [:usage => usage, :help => help])\n\n\nFORTRAN MANUAL\n SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )\n\n* Purpose\n* =======\n*\n* DLAGTS may be used to solve one of the systems of equations\n*\n* (T - lambda*I)*x = y or (T - lambda*I)'*x = y,\n*\n* where T is an n by n tridiagonal matrix, for x, following the\n* factorization of (T - lambda*I) as\n*\n* (T - lambda*I) = P*L*U ,\n*\n* by routine DLAGTF. The choice of equation to be solved is\n* controlled by the argument JOB, and in each case there is an option\n* to perturb zero or very small diagonal elements of U, this option\n* being intended for use in applications such as inverse iteration.\n*\n\n* Arguments\n* =========\n*\n* JOB (input) INTEGER\n* Specifies the job to be performed by DLAGTS as follows:\n* = 1: The equations (T - lambda*I)x = y are to be solved,\n* but diagonal elements of U are not to be perturbed.\n* = -1: The equations (T - lambda*I)x = y are to be solved\n* and, if overflow would otherwise occur, the diagonal\n* elements of U are to be perturbed. See argument TOL\n* below.\n* = 2: The equations (T - lambda*I)'x = y are to be solved,\n* but diagonal elements of U are not to be perturbed.\n* = -2: The equations (T - lambda*I)'x = y are to be solved\n* and, if overflow would otherwise occur, the diagonal\n* elements of U are to be perturbed. See argument TOL\n* below.\n*\n* N (input) INTEGER\n* The order of the matrix T.\n*\n* A (input) DOUBLE PRECISION array, dimension (N)\n* On entry, A must contain the diagonal elements of U as\n* returned from DLAGTF.\n*\n* B (input) DOUBLE PRECISION array, dimension (N-1)\n* On entry, B must contain the first super-diagonal elements of\n* U as returned from DLAGTF.\n*\n* C (input) DOUBLE PRECISION array, dimension (N-1)\n* On entry, C must contain the sub-diagonal elements of L as\n* returned from DLAGTF.\n*\n* D (input) DOUBLE PRECISION array, dimension (N-2)\n* On entry, D must contain the second super-diagonal elements\n* of U as returned from DLAGTF.\n*\n* IN (input) INTEGER array, dimension (N)\n* On entry, IN must contain details of the matrix P as returned\n* from DLAGTF.\n*\n* Y (input/output) DOUBLE PRECISION array, dimension (N)\n* On entry, the right hand side vector y.\n* On exit, Y is overwritten by the solution vector x.\n*\n* TOL (input/output) DOUBLE PRECISION\n* On entry, with JOB .lt. 0, TOL should be the minimum\n* perturbation to be made to very small diagonal elements of U.\n* TOL should normally be chosen as about eps*norm(U), where eps\n* is the relative machine precision, but if TOL is supplied as\n* non-positive, then it is reset to eps*max( abs( u(i,j) ) ).\n* If JOB .gt. 0 then TOL is not referenced.\n*\n* On exit, TOL is changed as described above, only if TOL is\n* non-positive on entry. Otherwise TOL is unchanged.\n*\n* INFO (output) INTEGER\n* = 0 : successful exit\n* .lt. 0: if INFO = -i, the i-th argument had an illegal value\n* .gt. 0: overflow would occur when computing the INFO(th)\n* element of the solution vector x. This can only occur\n* when JOB is supplied as positive and either means\n* that a diagonal element of U is very small, or that\n* the elements of the right-hand side vector y are very\n* large.\n*\n\n* =====================================================================\n*\n\n");
return Qnil;
}
if (rb_hash_aref(rblapack_options, sUsage) == Qtrue) {
printf("%s\n", "USAGE:\n info, y, tol = NumRu::Lapack.dlagts( job, a, b, c, d, in, y, tol, [:usage => usage, :help => help])\n");
return Qnil;
}
} else
rblapack_options = Qnil;
if (argc != 8 && argc != 8)
rb_raise(rb_eArgError,"wrong number of arguments (%d for 8)", argc);
rblapack_job = argv[0];
rblapack_a = argv[1];
rblapack_b = argv[2];
rblapack_c = argv[3];
rblapack_d = argv[4];
rblapack_in = argv[5];
rblapack_y = argv[6];
rblapack_tol = argv[7];
if (argc == 8) {
} else if (rblapack_options != Qnil) {
} else {
}
job = NUM2INT(rblapack_job);
if (!NA_IsNArray(rblapack_in))
rb_raise(rb_eArgError, "in (6th argument) must be NArray");
if (NA_RANK(rblapack_in) != 1)
rb_raise(rb_eArgError, "rank of in (6th argument) must be %d", 1);
n = NA_SHAPE0(rblapack_in);
if (NA_TYPE(rblapack_in) != NA_LINT)
rblapack_in = na_change_type(rblapack_in, NA_LINT);
in = NA_PTR_TYPE(rblapack_in, integer*);
tol = NUM2DBL(rblapack_tol);
if (!NA_IsNArray(rblapack_a))
rb_raise(rb_eArgError, "a (2th argument) must be NArray");
if (NA_RANK(rblapack_a) != 1)
rb_raise(rb_eArgError, "rank of a (2th argument) must be %d", 1);
if (NA_SHAPE0(rblapack_a) != n)
rb_raise(rb_eRuntimeError, "shape 0 of a must be the same as shape 0 of in");
if (NA_TYPE(rblapack_a) != NA_DFLOAT)
rblapack_a = na_change_type(rblapack_a, NA_DFLOAT);
a = NA_PTR_TYPE(rblapack_a, doublereal*);
if (!NA_IsNArray(rblapack_y))
rb_raise(rb_eArgError, "y (7th argument) must be NArray");
if (NA_RANK(rblapack_y) != 1)
rb_raise(rb_eArgError, "rank of y (7th argument) must be %d", 1);
if (NA_SHAPE0(rblapack_y) != n)
rb_raise(rb_eRuntimeError, "shape 0 of y must be the same as shape 0 of in");
if (NA_TYPE(rblapack_y) != NA_DFLOAT)
rblapack_y = na_change_type(rblapack_y, NA_DFLOAT);
y = NA_PTR_TYPE(rblapack_y, doublereal*);
if (!NA_IsNArray(rblapack_b))
rb_raise(rb_eArgError, "b (3th argument) must be NArray");
if (NA_RANK(rblapack_b) != 1)
rb_raise(rb_eArgError, "rank of b (3th argument) must be %d", 1);
if (NA_SHAPE0(rblapack_b) != (n-1))
rb_raise(rb_eRuntimeError, "shape 0 of b must be %d", n-1);
if (NA_TYPE(rblapack_b) != NA_DFLOAT)
rblapack_b = na_change_type(rblapack_b, NA_DFLOAT);
b = NA_PTR_TYPE(rblapack_b, doublereal*);
if (!NA_IsNArray(rblapack_d))
rb_raise(rb_eArgError, "d (5th argument) must be NArray");
if (NA_RANK(rblapack_d) != 1)
rb_raise(rb_eArgError, "rank of d (5th argument) must be %d", 1);
if (NA_SHAPE0(rblapack_d) != (n-2))
rb_raise(rb_eRuntimeError, "shape 0 of d must be %d", n-2);
if (NA_TYPE(rblapack_d) != NA_DFLOAT)
rblapack_d = na_change_type(rblapack_d, NA_DFLOAT);
d = NA_PTR_TYPE(rblapack_d, doublereal*);
if (!NA_IsNArray(rblapack_c))
rb_raise(rb_eArgError, "c (4th argument) must be NArray");
if (NA_RANK(rblapack_c) != 1)
rb_raise(rb_eArgError, "rank of c (4th argument) must be %d", 1);
if (NA_SHAPE0(rblapack_c) != (n-1))
rb_raise(rb_eRuntimeError, "shape 0 of c must be %d", n-1);
if (NA_TYPE(rblapack_c) != NA_DFLOAT)
rblapack_c = na_change_type(rblapack_c, NA_DFLOAT);
c = NA_PTR_TYPE(rblapack_c, doublereal*);
{
na_shape_t shape[1];
shape[0] = n;
rblapack_y_out__ = na_make_object(NA_DFLOAT, 1, shape, cNArray);
}
y_out__ = NA_PTR_TYPE(rblapack_y_out__, doublereal*);
MEMCPY(y_out__, y, doublereal, NA_TOTAL(rblapack_y));
rblapack_y = rblapack_y_out__;
y = y_out__;
dlagts_(&job, &n, a, b, c, d, in, y, &tol, &info);
rblapack_info = INT2NUM(info);
rblapack_tol = rb_float_new((double)tol);
return rb_ary_new3(3, rblapack_info, rblapack_y, rblapack_tol);
}
void
init_lapack_dlagts(VALUE mLapack, VALUE sH, VALUE sU, VALUE zero){
sHelp = sH;
sUsage = sU;
rblapack_ZERO = zero;
rb_define_module_function(mLapack, "dlagts", rblapack_dlagts, -1);
}
|
