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
|
--- dbleGSVD.cc.orig 2015-01-24 19:54:11 UTC
+++ dbleGSVD.cc
@@ -14,14 +14,14 @@
// You should have received a copy of the GNU General Public License along with
// this program; if not, see <http://www.gnu.org/licenses/>.
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "dbleGSVD.h"
#include <iostream>
-#include "dbleGSVD.h"
+#include <octave/oct.h>
+
#include "f77-fcn.h"
+#include "lo-error.h"
/*
uncomment those lines to monitor k and l
@@ -29,6 +29,13 @@
#include "pager.h"
*/
+#if defined (OCTAVE_HAVE_F77_INT_TYPE)
+# define TO_F77_INT(x) octave::to_f77_int (x)
+#else
+# define TO_F77_INT(x) (x)
+typedef octave_idx_type F77_INT;
+#endif
+
extern "C"
{
F77_RET_T
@@ -37,26 +44,26 @@ extern "C"
F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1
F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1
F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1
- const octave_idx_type&, // M (input) INTEGER
- const octave_idx_type&, // N (input) INTEGER
- const octave_idx_type&, // P (input) INTEGER
- octave_idx_type &, // K (output) INTEGER
- octave_idx_type &, // L (output) INTEGER
+ const F77_INT&, // M (input) INTEGER
+ const F77_INT&, // N (input) INTEGER
+ const F77_INT&, // P (input) INTEGER
+ F77_INT &, // K (output) INTEGER
+ F77_INT &, // L (output) INTEGER
double*, // A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
- const octave_idx_type&, // LDA (input) INTEGER
+ const F77_INT&, // LDA (input) INTEGER
double*, // B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
- const octave_idx_type&, // LDB (input) INTEGER
+ const F77_INT&, // LDB (input) INTEGER
double*, // ALPHA (output) DOUBLE PRECISION array, dimension (N)
double*, // BETA (output) DOUBLE PRECISION array, dimension (N)
double*, // U (output) DOUBLE PRECISION array, dimension (LDU,M)
- const octave_idx_type&, // LDU (input) INTEGER
+ const F77_INT&, // LDU (input) INTEGER
double*, // V (output) DOUBLE PRECISION array, dimension (LDV,P)
- const octave_idx_type&, // LDV (input) INTEGER
+ const F77_INT&, // LDV (input) INTEGER
double*, // Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
- const octave_idx_type&, // LDQ (input) INTEGER
+ const F77_INT&, // LDQ (input) INTEGER
double*, // WORK (workspace) DOUBLE PRECISION array
- int*, // IWORK (workspace/output) INTEGER array, dimension (N)
- octave_idx_type& // INFO (output)INTEGER
+ F77_INT*, // IWORK (workspace/output) INTEGER array, dimension (N)
+ F77_INT& // INFO (output)INTEGER
F77_CHAR_ARG_LEN_DECL
F77_CHAR_ARG_LEN_DECL
F77_CHAR_ARG_LEN_DECL
@@ -117,11 +124,11 @@ GSVD::R_matrix (void) const
octave_idx_type
GSVD::init (const Matrix& a, const Matrix& b, GSVD::type gsvd_type)
{
- octave_idx_type info;
+ F77_INT info;
- octave_idx_type m = a.rows ();
- octave_idx_type n = a.cols ();
- octave_idx_type p = b.rows ();
+ F77_INT m = TO_F77_INT (a.rows ());
+ F77_INT n = TO_F77_INT (a.cols ());
+ F77_INT p = TO_F77_INT (b.rows ());
Matrix atmp = a;
double *tmp_dataA = atmp.fortran_vec ();
@@ -129,17 +136,17 @@ GSVD::init (const Matrix& a, const Matri
Matrix btmp = b;
double *tmp_dataB = btmp.fortran_vec ();
- // octave_idx_type min_mn = m < n ? m : n;
+ // F77_INT min_mn = m < n ? m : n;
char jobu = 'U';
char jobv = 'V';
char jobq = 'Q';
- octave_idx_type nrow_u = m;
- octave_idx_type nrow_v = p;
- octave_idx_type nrow_q = n;
+ F77_INT nrow_u = m;
+ F77_INT nrow_v = p;
+ F77_INT nrow_q = n;
- octave_idx_type k, l;
+ F77_INT k, l;
switch (gsvd_type)
{
@@ -183,14 +190,14 @@ GSVD::init (const Matrix& a, const Matri
}
double *q = right_sm.fortran_vec ();
- octave_idx_type lwork = 3*n;
+ F77_INT lwork = 3*n;
lwork = lwork > m ? lwork : m;
lwork = (lwork > p ? lwork : p) + n;
Array<double> work (dim_vector (lwork, 1));
Array<double> alpha (dim_vector (n, 1));
Array<double> beta (dim_vector (n, 1));
- Array<int> iwork (dim_vector (n, 1));
+ Array<F77_INT> iwork (dim_vector (n, 1));
F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
F77_CONST_CHAR_ARG2 (&jobv, 1),
@@ -213,13 +220,13 @@ GSVD::init (const Matrix& a, const Matri
if (info > 0) {
(*current_liboctave_error_handler) ("dggsvd.f: Jacobi-type procedure failed to converge.");
} else {
- octave_idx_type i, j;
+ F77_INT i, j;
if (GSVD::std == gsvd_type) {
R.resize(k+l, k+l);
- int astart = n-k-l;
+ F77_INT astart = n-k-l;
if (m - k - l >= 0) {
- int astart = n-k-l;
+ F77_INT astart = n-k-l;
/*
* R is stored in A(1:K+L,N-K-L+1:N)
*/
|
