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
|
#include "rb_lapack.h"
extern VOID claed0_(integer* qsiz, integer* n, real* d, real* e, complex* q, integer* ldq, complex* qstore, integer* ldqs, real* rwork, integer* iwork, integer* info);
static VALUE
rblapack_claed0(int argc, VALUE *argv, VALUE self){
VALUE rblapack_qsiz;
integer qsiz;
VALUE rblapack_d;
real *d;
VALUE rblapack_e;
real *e;
VALUE rblapack_q;
complex *q;
VALUE rblapack_info;
integer info;
VALUE rblapack_d_out__;
real *d_out__;
VALUE rblapack_e_out__;
real *e_out__;
VALUE rblapack_q_out__;
complex *q_out__;
complex *qstore;
real *rwork;
integer *iwork;
integer n;
integer ldq;
integer ldqs;
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, d, e, q = NumRu::Lapack.claed0( qsiz, d, e, q, [:usage => usage, :help => help])\n\n\nFORTRAN MANUAL\n SUBROUTINE CLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, IWORK, INFO )\n\n* Purpose\n* =======\n*\n* Using the divide and conquer method, CLAED0 computes all eigenvalues\n* of a symmetric tridiagonal matrix which is one diagonal block of\n* those from reducing a dense or band Hermitian matrix and\n* corresponding eigenvectors of the dense or band matrix.\n*\n\n* Arguments\n* =========\n*\n* QSIZ (input) INTEGER\n* The dimension of the unitary matrix used to reduce\n* the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.\n*\n* N (input) INTEGER\n* The dimension of the symmetric tridiagonal matrix. N >= 0.\n*\n* D (input/output) REAL array, dimension (N)\n* On entry, the diagonal elements of the tridiagonal matrix.\n* On exit, the eigenvalues in ascending order.\n*\n* E (input/output) REAL array, dimension (N-1)\n* On entry, the off-diagonal elements of the tridiagonal matrix.\n* On exit, E has been destroyed.\n*\n* Q (input/output) COMPLEX array, dimension (LDQ,N)\n* On entry, Q must contain an QSIZ x N matrix whose columns\n* unitarily orthonormal. It is a part of the unitary matrix\n* that reduces the full dense Hermitian matrix to a\n* (reducible) symmetric tridiagonal matrix.\n*\n* LDQ (input) INTEGER\n* The leading dimension of the array Q. LDQ >= max(1,N).\n*\n* IWORK (workspace) INTEGER array,\n* the dimension of IWORK must be at least\n* 6 + 6*N + 5*N*lg N\n* ( lg( N ) = smallest integer k\n* such that 2^k >= N )\n*\n* RWORK (workspace) REAL array,\n* dimension (1 + 3*N + 2*N*lg N + 3*N**2)\n* ( lg( N ) = smallest integer k\n* such that 2^k >= N )\n*\n* QSTORE (workspace) COMPLEX array, dimension (LDQS, N)\n* Used to store parts of\n* the eigenvector matrix when the updating matrix multiplies\n* take place.\n*\n* LDQS (input) INTEGER\n* The leading dimension of the array QSTORE.\n* LDQS >= max(1,N).\n*\n* INFO (output) INTEGER\n* = 0: successful exit.\n* < 0: if INFO = -i, the i-th argument had an illegal value.\n* > 0: The algorithm failed to compute an eigenvalue while\n* working on the submatrix lying in rows and columns\n* INFO/(N+1) through mod(INFO,N+1).\n*\n\n* =====================================================================\n*\n* Warning: N could be as big as QSIZ!\n*\n\n");
return Qnil;
}
if (rb_hash_aref(rblapack_options, sUsage) == Qtrue) {
printf("%s\n", "USAGE:\n info, d, e, q = NumRu::Lapack.claed0( qsiz, d, e, q, [:usage => usage, :help => help])\n");
return Qnil;
}
} else
rblapack_options = Qnil;
if (argc != 4 && argc != 4)
rb_raise(rb_eArgError,"wrong number of arguments (%d for 4)", argc);
rblapack_qsiz = argv[0];
rblapack_d = argv[1];
rblapack_e = argv[2];
rblapack_q = argv[3];
if (argc == 4) {
} else if (rblapack_options != Qnil) {
} else {
}
qsiz = NUM2INT(rblapack_qsiz);
if (!NA_IsNArray(rblapack_q))
rb_raise(rb_eArgError, "q (4th argument) must be NArray");
if (NA_RANK(rblapack_q) != 2)
rb_raise(rb_eArgError, "rank of q (4th 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*);
if (!NA_IsNArray(rblapack_d))
rb_raise(rb_eArgError, "d (2th argument) must be NArray");
if (NA_RANK(rblapack_d) != 1)
rb_raise(rb_eArgError, "rank of d (2th argument) must be %d", 1);
if (NA_SHAPE0(rblapack_d) != n)
rb_raise(rb_eRuntimeError, "shape 0 of d must be the same as shape 1 of q");
if (NA_TYPE(rblapack_d) != NA_SFLOAT)
rblapack_d = na_change_type(rblapack_d, NA_SFLOAT);
d = NA_PTR_TYPE(rblapack_d, real*);
ldqs = MAX(1,n);
if (!NA_IsNArray(rblapack_e))
rb_raise(rb_eArgError, "e (3th argument) must be NArray");
if (NA_RANK(rblapack_e) != 1)
rb_raise(rb_eArgError, "rank of e (3th argument) must be %d", 1);
if (NA_SHAPE0(rblapack_e) != (n-1))
rb_raise(rb_eRuntimeError, "shape 0 of e must be %d", n-1);
if (NA_TYPE(rblapack_e) != NA_SFLOAT)
rblapack_e = na_change_type(rblapack_e, NA_SFLOAT);
e = NA_PTR_TYPE(rblapack_e, real*);
{
na_shape_t shape[1];
shape[0] = n;
rblapack_d_out__ = na_make_object(NA_SFLOAT, 1, shape, cNArray);
}
d_out__ = NA_PTR_TYPE(rblapack_d_out__, real*);
MEMCPY(d_out__, d, real, NA_TOTAL(rblapack_d));
rblapack_d = rblapack_d_out__;
d = d_out__;
{
na_shape_t shape[1];
shape[0] = n-1;
rblapack_e_out__ = na_make_object(NA_SFLOAT, 1, shape, cNArray);
}
e_out__ = NA_PTR_TYPE(rblapack_e_out__, real*);
MEMCPY(e_out__, e, real, NA_TOTAL(rblapack_e));
rblapack_e = rblapack_e_out__;
e = e_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__;
qstore = ALLOC_N(complex, (ldqs)*(n));
rwork = ALLOC_N(real, (1 + 3*n + 2*n*LG(n) + 3*pow(n,2)));
iwork = ALLOC_N(integer, (6 + 6*n + 5*n*LG(n)));
claed0_(&qsiz, &n, d, e, q, &ldq, qstore, &ldqs, rwork, iwork, &info);
free(qstore);
free(rwork);
free(iwork);
rblapack_info = INT2NUM(info);
return rb_ary_new3(4, rblapack_info, rblapack_d, rblapack_e, rblapack_q);
}
void
init_lapack_claed0(VALUE mLapack, VALUE sH, VALUE sU, VALUE zero){
sHelp = sH;
sUsage = sU;
rblapack_ZERO = zero;
rb_define_module_function(mLapack, "claed0", rblapack_claed0, -1);
}
|
