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
140
141
142
143
144
145
146
147
148
149
150
151
|
#include "rb_lapack.h"
extern VOID zgelss_(integer* m, integer* n, integer* nrhs, doublecomplex* a, integer* lda, doublecomplex* b, integer* ldb, doublereal* s, doublereal* rcond, integer* rank, doublecomplex* work, integer* lwork, doublereal* rwork, integer* info);
static VALUE
rblapack_zgelss(int argc, VALUE *argv, VALUE self){
VALUE rblapack_a;
doublecomplex *a;
VALUE rblapack_b;
doublecomplex *b;
VALUE rblapack_rcond;
doublereal rcond;
VALUE rblapack_lwork;
integer lwork;
VALUE rblapack_s;
doublereal *s;
VALUE rblapack_rank;
integer rank;
VALUE rblapack_work;
doublecomplex *work;
VALUE rblapack_info;
integer info;
VALUE rblapack_a_out__;
doublecomplex *a_out__;
VALUE rblapack_b_out__;
doublecomplex *b_out__;
doublereal *rwork;
integer lda;
integer n;
integer m;
integer nrhs;
integer ldb;
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 s, rank, work, info, a, b = NumRu::Lapack.zgelss( a, b, rcond, [:lwork => lwork, :usage => usage, :help => help])\n\n\nFORTRAN MANUAL\n SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO )\n\n* Purpose\n* =======\n*\n* ZGELSS computes the minimum norm solution to a complex linear\n* least squares problem:\n*\n* Minimize 2-norm(| b - A*x |).\n*\n* using the singular value decomposition (SVD) of A. A is an M-by-N\n* matrix which may be rank-deficient.\n*\n* Several right hand side vectors b and solution vectors x can be\n* handled in a single call; they are stored as the columns of the\n* M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix\n* X.\n*\n* The effective rank of A is determined by treating as zero those\n* singular values which are less than RCOND times the largest singular\n* value.\n*\n\n* Arguments\n* =========\n*\n* M (input) INTEGER\n* The number of rows of the matrix A. M >= 0.\n*\n* N (input) INTEGER\n* The number of columns of the matrix A. N >= 0.\n*\n* NRHS (input) INTEGER\n* The number of right hand sides, i.e., the number of columns\n* of the matrices B and X. NRHS >= 0.\n*\n* A (input/output) COMPLEX*16 array, dimension (LDA,N)\n* On entry, the M-by-N matrix A.\n* On exit, the first min(m,n) rows of A are overwritten with\n* its right singular vectors, stored rowwise.\n*\n* LDA (input) INTEGER\n* The leading dimension of the array A. LDA >= max(1,M).\n*\n* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)\n* On entry, the M-by-NRHS right hand side matrix B.\n* On exit, B is overwritten by the N-by-NRHS solution matrix X.\n* If m >= n and RANK = n, the residual sum-of-squares for\n* the solution in the i-th column is given by the sum of\n* squares of the modulus of elements n+1:m in that column.\n*\n* LDB (input) INTEGER\n* The leading dimension of the array B. LDB >= max(1,M,N).\n*\n* S (output) DOUBLE PRECISION array, dimension (min(M,N))\n* The singular values of A in decreasing order.\n* The condition number of A in the 2-norm = S(1)/S(min(m,n)).\n*\n* RCOND (input) DOUBLE PRECISION\n* RCOND is used to determine the effective rank of A.\n* Singular values S(i) <= RCOND*S(1) are treated as zero.\n* If RCOND < 0, machine precision is used instead.\n*\n* RANK (output) INTEGER\n* The effective rank of A, i.e., the number of singular values\n* which are greater than RCOND*S(1).\n*\n* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))\n* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.\n*\n* LWORK (input) INTEGER\n* The dimension of the array WORK. LWORK >= 1, and also:\n* LWORK >= 2*min(M,N) + max(M,N,NRHS)\n* For good performance, LWORK should generally be larger.\n*\n* If LWORK = -1, then a workspace query is assumed; the routine\n* only calculates the optimal size of the WORK array, returns\n* this value as the first entry of the WORK array, and no error\n* message related to LWORK is issued by XERBLA.\n*\n* RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,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 for computing the SVD failed to converge;\n* if INFO = i, i off-diagonal elements of an intermediate\n* bidiagonal form did not converge to zero.\n*\n\n* =====================================================================\n*\n\n");
return Qnil;
}
if (rb_hash_aref(rblapack_options, sUsage) == Qtrue) {
printf("%s\n", "USAGE:\n s, rank, work, info, a, b = NumRu::Lapack.zgelss( a, b, rcond, [:lwork => lwork, :usage => usage, :help => help])\n");
return Qnil;
}
} else
rblapack_options = Qnil;
if (argc != 3 && argc != 4)
rb_raise(rb_eArgError,"wrong number of arguments (%d for 3)", argc);
rblapack_a = argv[0];
rblapack_b = argv[1];
rblapack_rcond = argv[2];
if (argc == 4) {
rblapack_lwork = argv[3];
} else if (rblapack_options != Qnil) {
rblapack_lwork = rb_hash_aref(rblapack_options, ID2SYM(rb_intern("lwork")));
} else {
rblapack_lwork = Qnil;
}
if (!NA_IsNArray(rblapack_a))
rb_raise(rb_eArgError, "a (1th argument) must be NArray");
if (NA_RANK(rblapack_a) != 2)
rb_raise(rb_eArgError, "rank of a (1th argument) must be %d", 2);
lda = NA_SHAPE0(rblapack_a);
n = NA_SHAPE1(rblapack_a);
if (NA_TYPE(rblapack_a) != NA_DCOMPLEX)
rblapack_a = na_change_type(rblapack_a, NA_DCOMPLEX);
a = NA_PTR_TYPE(rblapack_a, doublecomplex*);
rcond = NUM2DBL(rblapack_rcond);
m = lda;
if (!NA_IsNArray(rblapack_b))
rb_raise(rb_eArgError, "b (2th argument) must be NArray");
if (NA_RANK(rblapack_b) != 2)
rb_raise(rb_eArgError, "rank of b (2th argument) must be %d", 2);
if (NA_SHAPE0(rblapack_b) != m)
rb_raise(rb_eRuntimeError, "shape 0 of b must be lda");
nrhs = NA_SHAPE1(rblapack_b);
if (NA_TYPE(rblapack_b) != NA_DCOMPLEX)
rblapack_b = na_change_type(rblapack_b, NA_DCOMPLEX);
b = NA_PTR_TYPE(rblapack_b, doublecomplex*);
ldb = MAX(m, n);
if (rblapack_lwork == Qnil)
lwork = 3*MIN(m,n) + MAX(MAX(2*MIN(m,n),MAX(m,n)),nrhs);
else {
lwork = NUM2INT(rblapack_lwork);
}
{
na_shape_t shape[1];
shape[0] = MIN(m,n);
rblapack_s = na_make_object(NA_DFLOAT, 1, shape, cNArray);
}
s = NA_PTR_TYPE(rblapack_s, doublereal*);
{
na_shape_t shape[1];
shape[0] = MAX(1,lwork);
rblapack_work = na_make_object(NA_DCOMPLEX, 1, shape, cNArray);
}
work = NA_PTR_TYPE(rblapack_work, doublecomplex*);
{
na_shape_t shape[2];
shape[0] = lda;
shape[1] = n;
rblapack_a_out__ = na_make_object(NA_DCOMPLEX, 2, shape, cNArray);
}
a_out__ = NA_PTR_TYPE(rblapack_a_out__, doublecomplex*);
MEMCPY(a_out__, a, doublecomplex, NA_TOTAL(rblapack_a));
rblapack_a = rblapack_a_out__;
a = a_out__;
{
na_shape_t shape[2];
shape[0] = MAX(m, n);
shape[1] = nrhs;
rblapack_b_out__ = na_make_object(NA_DCOMPLEX, 2, shape, cNArray);
}
b_out__ = NA_PTR_TYPE(rblapack_b_out__, doublecomplex*);
{
VALUE __shape__[3];
__shape__[0] = m < n ? rb_range_new(rblapack_ZERO, INT2NUM(m), Qtrue) : Qtrue;
__shape__[1] = Qtrue;
__shape__[2] = rblapack_b;
na_aset(3, __shape__, rblapack_b_out__);
}
rblapack_b = rblapack_b_out__;
b = b_out__;
rwork = ALLOC_N(doublereal, (5*MIN(m,n)));
zgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, work, &lwork, rwork, &info);
free(rwork);
rblapack_rank = INT2NUM(rank);
rblapack_info = INT2NUM(info);
{
VALUE __shape__[2];
__shape__[0] = m < n ? Qtrue : rb_range_new(rblapack_ZERO, INT2NUM(n), Qtrue);
__shape__[1] = Qtrue;
rblapack_b = na_aref(2, __shape__, rblapack_b);
}
return rb_ary_new3(6, rblapack_s, rblapack_rank, rblapack_work, rblapack_info, rblapack_a, rblapack_b);
}
void
init_lapack_zgelss(VALUE mLapack, VALUE sH, VALUE sU, VALUE zero){
sHelp = sH;
sUsage = sU;
rblapack_ZERO = zero;
rb_define_module_function(mLapack, "zgelss", rblapack_zgelss, -1);
}
|
