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
|
#include "rb_lapack.h"
extern VOID chbtrd_(char* vect, char* uplo, integer* n, integer* kd, complex* ab, integer* ldab, real* d, real* e, complex* q, integer* ldq, complex* work, integer* info);
static VALUE
rblapack_chbtrd(int argc, VALUE *argv, VALUE self){
VALUE rblapack_vect;
char vect;
VALUE rblapack_uplo;
char uplo;
VALUE rblapack_kd;
integer kd;
VALUE rblapack_ab;
complex *ab;
VALUE rblapack_q;
complex *q;
VALUE rblapack_d;
real *d;
VALUE rblapack_e;
real *e;
VALUE rblapack_info;
integer info;
VALUE rblapack_ab_out__;
complex *ab_out__;
VALUE rblapack_q_out__;
complex *q_out__;
complex *work;
integer ldab;
integer n;
integer ldq;
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 d, e, info, ab, q = NumRu::Lapack.chbtrd( vect, uplo, kd, ab, q, [:usage => usage, :help => help])\n\n\nFORTRAN MANUAL\n SUBROUTINE CHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO )\n\n* Purpose\n* =======\n*\n* CHBTRD reduces a complex Hermitian band matrix A to real symmetric\n* tridiagonal form T by a unitary similarity transformation:\n* Q**H * A * Q = T.\n*\n\n* Arguments\n* =========\n*\n* VECT (input) CHARACTER*1\n* = 'N': do not form Q;\n* = 'V': form Q;\n* = 'U': update a matrix X, by forming X*Q.\n*\n* UPLO (input) CHARACTER*1\n* = 'U': Upper triangle of A is stored;\n* = 'L': Lower triangle of A is stored.\n*\n* N (input) INTEGER\n* The order of the matrix A. N >= 0.\n*\n* KD (input) INTEGER\n* The number of superdiagonals of the matrix A if UPLO = 'U',\n* or the number of subdiagonals if UPLO = 'L'. KD >= 0.\n*\n* AB (input/output) COMPLEX array, dimension (LDAB,N)\n* On entry, the upper or lower triangle of the Hermitian band\n* matrix A, stored in the first KD+1 rows of the array. The\n* j-th column of A is stored in the j-th column of the array AB\n* as follows:\n* if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;\n* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).\n* On exit, the diagonal elements of AB are overwritten by the\n* diagonal elements of the tridiagonal matrix T; if KD > 0, the\n* elements on the first superdiagonal (if UPLO = 'U') or the\n* first subdiagonal (if UPLO = 'L') are overwritten by the\n* off-diagonal elements of T; the rest of AB is overwritten by\n* values generated during the reduction.\n*\n* LDAB (input) INTEGER\n* The leading dimension of the array AB. LDAB >= KD+1.\n*\n* D (output) REAL array, dimension (N)\n* The diagonal elements of the tridiagonal matrix T.\n*\n* E (output) REAL array, dimension (N-1)\n* The off-diagonal elements of the tridiagonal matrix T:\n* E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.\n*\n* Q (input/output) COMPLEX array, dimension (LDQ,N)\n* On entry, if VECT = 'U', then Q must contain an N-by-N\n* matrix X; if VECT = 'N' or 'V', then Q need not be set.\n*\n* On exit:\n* if VECT = 'V', Q contains the N-by-N unitary matrix Q;\n* if VECT = 'U', Q contains the product X*Q;\n* if VECT = 'N', the array Q is not referenced.\n*\n* LDQ (input) INTEGER\n* The leading dimension of the array Q.\n* LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.\n*\n* WORK (workspace) COMPLEX array, dimension (N)\n*\n* INFO (output) INTEGER\n* = 0: successful exit\n* < 0: if INFO = -i, the i-th argument had an illegal value\n*\n\n* Further Details\n* ===============\n*\n* Modified by Linda Kaufman, Bell Labs.\n*\n* =====================================================================\n*\n\n");
return Qnil;
}
if (rb_hash_aref(rblapack_options, sUsage) == Qtrue) {
printf("%s\n", "USAGE:\n d, e, info, ab, q = NumRu::Lapack.chbtrd( vect, uplo, kd, ab, q, [:usage => usage, :help => help])\n");
return Qnil;
}
} else
rblapack_options = Qnil;
if (argc != 5 && argc != 5)
rb_raise(rb_eArgError,"wrong number of arguments (%d for 5)", argc);
rblapack_vect = argv[0];
rblapack_uplo = argv[1];
rblapack_kd = argv[2];
rblapack_ab = argv[3];
rblapack_q = argv[4];
if (argc == 5) {
} else if (rblapack_options != Qnil) {
} else {
}
vect = StringValueCStr(rblapack_vect)[0];
kd = NUM2INT(rblapack_kd);
if (!NA_IsNArray(rblapack_q))
rb_raise(rb_eArgError, "q (5th argument) must be NArray");
if (NA_RANK(rblapack_q) != 2)
rb_raise(rb_eArgError, "rank of q (5th argument) must be %d", 2);
ldq = NA_SHAPE0(rblapack_q);
n = NA_SHAPE1(rblapack_q);
if (NA_TYPE(rblapack_q) != NA_SCOMPLEX)
rblapack_q = na_change_type(rblapack_q, NA_SCOMPLEX);
q = NA_PTR_TYPE(rblapack_q, complex*);
uplo = StringValueCStr(rblapack_uplo)[0];
if (!NA_IsNArray(rblapack_ab))
rb_raise(rb_eArgError, "ab (4th argument) must be NArray");
if (NA_RANK(rblapack_ab) != 2)
rb_raise(rb_eArgError, "rank of ab (4th argument) must be %d", 2);
ldab = NA_SHAPE0(rblapack_ab);
if (NA_SHAPE1(rblapack_ab) != n)
rb_raise(rb_eRuntimeError, "shape 1 of ab must be the same as shape 1 of q");
if (NA_TYPE(rblapack_ab) != NA_SCOMPLEX)
rblapack_ab = na_change_type(rblapack_ab, NA_SCOMPLEX);
ab = NA_PTR_TYPE(rblapack_ab, complex*);
{
na_shape_t shape[1];
shape[0] = n;
rblapack_d = na_make_object(NA_SFLOAT, 1, shape, cNArray);
}
d = NA_PTR_TYPE(rblapack_d, real*);
{
na_shape_t shape[1];
shape[0] = n-1;
rblapack_e = na_make_object(NA_SFLOAT, 1, shape, cNArray);
}
e = NA_PTR_TYPE(rblapack_e, real*);
{
na_shape_t shape[2];
shape[0] = ldab;
shape[1] = n;
rblapack_ab_out__ = na_make_object(NA_SCOMPLEX, 2, shape, cNArray);
}
ab_out__ = NA_PTR_TYPE(rblapack_ab_out__, complex*);
MEMCPY(ab_out__, ab, complex, NA_TOTAL(rblapack_ab));
rblapack_ab = rblapack_ab_out__;
ab = ab_out__;
{
na_shape_t shape[2];
shape[0] = ldq;
shape[1] = n;
rblapack_q_out__ = na_make_object(NA_SCOMPLEX, 2, shape, cNArray);
}
q_out__ = NA_PTR_TYPE(rblapack_q_out__, complex*);
MEMCPY(q_out__, q, complex, NA_TOTAL(rblapack_q));
rblapack_q = rblapack_q_out__;
q = q_out__;
work = ALLOC_N(complex, (n));
chbtrd_(&vect, &uplo, &n, &kd, ab, &ldab, d, e, q, &ldq, work, &info);
free(work);
rblapack_info = INT2NUM(info);
return rb_ary_new3(5, rblapack_d, rblapack_e, rblapack_info, rblapack_ab, rblapack_q);
}
void
init_lapack_chbtrd(VALUE mLapack, VALUE sH, VALUE sU, VALUE zero){
sHelp = sH;
sUsage = sU;
rblapack_ZERO = zero;
rb_define_module_function(mLapack, "chbtrd", rblapack_chbtrd, -1);
}
|
