-
Notifications
You must be signed in to change notification settings - Fork 0
/
gsvd.cpp
149 lines (112 loc) · 4.3 KB
/
gsvd.cpp
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
/**
* @date Tue Jan 10 21:14 2016 +0100
* @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
*
* @brief Binds the Generalized SVD
*/
#include "gsvd.h"
#include <bob.blitz/cppapi.h>
#include <bob.blitz/cleanup.h>
#include <bob.math/gsvd.h>
#include <bob.math/svd.h>
#include <bob.math/linear.h>
PyObject* py_gsvd (PyObject*, PyObject* args, PyObject* kwds) {
/* Parses input arguments in a single shot */
static const char* const_kwlist[] = { "A", "B", 0 /* Sentinel */ };
static char** kwlist = const_cast<char**>(const_kwlist);
PyBlitzArrayObject* A = 0;
PyBlitzArrayObject* B = 0;
if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&", kwlist,
&PyBlitzArray_Converter, &A,
&PyBlitzArray_Converter, &B
))
return 0;
auto A_ = make_safe(A);
auto B_ = make_safe(B);
if (A->ndim != 2 || A->type_num != NPY_FLOAT64) {
PyErr_Format(PyExc_TypeError, "`A` matrix only supports 2D 64-bit float array");
return 0;
}
if (B->ndim != 2 || B->type_num != NPY_FLOAT64) {
PyErr_Format(PyExc_TypeError, "`B` matrix only supports 2D 64-bit float array");
return 0;
}
auto A_bz = PyBlitzArrayCxx_AsBlitz<double,2>(A);
auto B_bz = PyBlitzArrayCxx_AsBlitz<double,2>(B);
const int M = A_bz->extent(0);
const int N = A_bz->extent(1);
const int P = B_bz->extent(0);
// Creating the output matrices
blitz::Array<double,2> U(M, M); U=0;
blitz::Array<double,2> V(P, P); V=0;
blitz::Array<double,2> Q(N, N); Q=0;
blitz::Array<double,2> zeroR(N, N); zeroR=0;
blitz::Array<double,2> X(N, N); X=0;
blitz::Array<double,2> C(N, N); C=0;
blitz::Array<double,2> S(N, N); S=0;
try {
bob::math::gsvd(*A_bz,*B_bz,U,V,zeroR,Q,X,C,S);
return Py_BuildValue("NNNNN",
PyBlitzArrayCxx_AsConstNumpy(U),
PyBlitzArrayCxx_AsConstNumpy(V),
PyBlitzArrayCxx_AsConstNumpy(X),
PyBlitzArrayCxx_AsConstNumpy(C),
PyBlitzArrayCxx_AsConstNumpy(S));
}
catch (std::exception& e) {
PyErr_SetString(PyExc_RuntimeError, e.what());
}
catch (...) {
PyErr_SetString(PyExc_RuntimeError, "gsvd failed: unknown exception caught");
}
return 0;
}
PyObject* py_svd (PyObject*, PyObject* args, PyObject* kwds) {
/* Parses input arguments in a single shot */
static const char* const_kwlist[] = { "A", 0 /* Sentinel */ };
static char** kwlist = const_cast<char**>(const_kwlist);
PyBlitzArrayObject* A = 0;
if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&", kwlist,
&PyBlitzArray_Converter, &A
))
return 0;
auto A_ = make_safe(A);
if (A->ndim != 2 || A->type_num != NPY_FLOAT64) {
PyErr_Format(PyExc_TypeError, "`A` matrix only supports 2D 64-bit float array");
return 0;
}
auto A_bz = PyBlitzArrayCxx_AsBlitz<double,2>(A);
int M = A_bz->extent(0);
int N = A_bz->extent(1);
// Creating the output matrices
// Prepares to call LAPACK function:
// We recall that FORTRAN/LAPACK is column-major order whereas blitz arrays
// are row-major order by default.
// If A = U.S.V^T, then A^T = V.S.U^T
blitz::Array<double,2> V(M, M); V=0;
blitz::Array<double,1> S(std::min(M,N)); S=0;
blitz::Array<double,2> U(N, N); U=0;
try {
bob::math::svd(*A_bz,V,S,U, true);
// S for the python output.
// LAPACK returns an 1d matrix of size n.
// In order to nicelly fit A = U S V,
// I need to create a matrix S_output of size MxN and fits the diagonal of S (NxN) on it.
blitz::Array<double,2> S_output(M, N); S_output=0;
blitz::Array<double,2> S_diag(N, N); S_diag=0;
bob::math::diag(S, S_diag);
S_output(blitz::Range(0,N-1), blitz::Range(0,N-1)) = S_diag;
//bob::math::diag(S, (reinterpret_cast<blitz::Array<double,2>>)S_output(blitz::Range(0,N-1), blitz::Range(0,N-1)));
return Py_BuildValue("NNN",
PyBlitzArrayCxx_AsConstNumpy(U),
PyBlitzArrayCxx_AsConstNumpy(S_output),
PyBlitzArrayCxx_AsConstNumpy(V));
}
catch (std::exception& e) {
PyErr_SetString(PyExc_RuntimeError, e.what());
}
catch (...) {
PyErr_SetString(PyExc_RuntimeError, "svd failed: unknown exception caught");
}
return 0;
}