1
/*This module contributed by Doug Heisterkamp
2
Modified by Jim Hugunin
3
More modifications by Jeff Whitaker
4
*/
5
#define NPY_NO_DEPRECATED_API NPY_API_VERSION
6

7
#include "Python.h"
8
#include "numpy/arrayobject.h"
9
#include "npy_cblas.h"
10

11

12
#define FNAME(name) BLAS_FUNC(name)
13

14
typedef CBLAS_INT        fortran_int;
15

16
#ifdef HAVE_BLAS_ILP64
17

18
#if NPY_BITSOF_SHORT == 64
19
#define FINT_PYFMT       "h"
20
#elif NPY_BITSOF_INT == 64
21
#define FINT_PYFMT       "i"
22
#elif NPY_BITSOF_LONG == 64
23
#define FINT_PYFMT       "l"
24
#elif NPY_BITSOF_LONGLONG == 64
25
#define FINT_PYFMT       "L"
26
#else
27
#error No compatible 64-bit integer size. \
28
       Please contact NumPy maintainers and give detailed information about your \
29
       compiler and platform, or set NPY_USE_BLAS64_=0
30
#endif
31

32
#else
33
#define FINT_PYFMT       "i"
34
#endif
35

36
typedef struct { float r, i; } f2c_complex;
37
typedef struct { double r, i; } f2c_doublecomplex;
38
/* typedef long int (*L_fp)(); */
39

40
extern fortran_int FNAME(dgelsd)(fortran_int *m, fortran_int *n, fortran_int *nrhs,
41
                          double a[], fortran_int *lda, double b[], fortran_int *ldb,
42
                          double s[], double *rcond, fortran_int *rank,
43
                          double work[], fortran_int *lwork, fortran_int iwork[], fortran_int *info);
44

45
extern fortran_int FNAME(zgelsd)(fortran_int *m, fortran_int *n, fortran_int *nrhs,
46
                          f2c_doublecomplex a[], fortran_int *lda,
47
                          f2c_doublecomplex b[], fortran_int *ldb,
48
                          double s[], double *rcond, fortran_int *rank,
49
                          f2c_doublecomplex work[], fortran_int *lwork,
50
                          double rwork[], fortran_int iwork[], fortran_int *info);
51

52
extern fortran_int FNAME(dgeqrf)(fortran_int *m, fortran_int *n, double a[], fortran_int *lda,
53
                          double tau[], double work[],
54
                          fortran_int *lwork, fortran_int *info);
55

56
extern fortran_int FNAME(zgeqrf)(fortran_int *m, fortran_int *n, f2c_doublecomplex a[], fortran_int *lda,
57
                          f2c_doublecomplex tau[], f2c_doublecomplex work[],
58
                          fortran_int *lwork, fortran_int *info);
59

60
extern fortran_int FNAME(dorgqr)(fortran_int *m, fortran_int *n, fortran_int *k, double a[], fortran_int *lda,
61
                          double tau[], double work[],
62
                          fortran_int *lwork, fortran_int *info);
63

64
extern fortran_int FNAME(zungqr)(fortran_int *m, fortran_int *n, fortran_int *k, f2c_doublecomplex a[],
65
                          fortran_int *lda, f2c_doublecomplex tau[],
66
                          f2c_doublecomplex work[], fortran_int *lwork, fortran_int *info);
67

68
extern fortran_int FNAME(xerbla)(char *srname, fortran_int *info);
69

70
static PyObject *LapackError;
71

72
#define TRY(E) if (!(E)) return NULL
73

74
static int
75 1
check_object(PyObject *ob, int t, char *obname,
76
                        char *tname, char *funname)
77
{
78 1
    if (!PyArray_Check(ob)) {
79 0
        PyErr_Format(LapackError,
80
                     "Expected an array for parameter %s in lapack_lite.%s",
81
                     obname, funname);
82 0
        return 0;
83
    }
84 1
    else if (!PyArray_IS_C_CONTIGUOUS((PyArrayObject *)ob)) {
85 0
        PyErr_Format(LapackError,
86
                     "Parameter %s is not contiguous in lapack_lite.%s",
87
                     obname, funname);
88 0
        return 0;
89
    }
90 1
    else if (!(PyArray_TYPE((PyArrayObject *)ob) == t)) {
91 0
        PyErr_Format(LapackError,
92
                     "Parameter %s is not of type %s in lapack_lite.%s",
93
                     obname, tname, funname);
94 0
        return 0;
95
    }
96 1
    else if (PyArray_ISBYTESWAPPED((PyArrayObject *)ob)) {
97 0
        PyErr_Format(LapackError,
98
                     "Parameter %s has non-native byte order in lapack_lite.%s",
99
                     obname, funname);
100 0
        return 0;
101
    }
102
    else {
103
        return 1;
104
    }
105
}
106

107
#define CHDATA(p) ((char *) PyArray_DATA((PyArrayObject *)p))
108
#define SHDATA(p) ((short int *) PyArray_DATA((PyArrayObject *)p))
109
#define DDATA(p) ((double *) PyArray_DATA((PyArrayObject *)p))
110
#define FDATA(p) ((float *) PyArray_DATA((PyArrayObject *)p))
111
#define CDATA(p) ((f2c_complex *) PyArray_DATA((PyArrayObject *)p))
112
#define ZDATA(p) ((f2c_doublecomplex *) PyArray_DATA((PyArrayObject *)p))
113
#define IDATA(p) ((fortran_int *) PyArray_DATA((PyArrayObject *)p))
114

115
static PyObject *
116 0
lapack_lite_dgelsd(PyObject *NPY_UNUSED(self), PyObject *args)
117
{
118
    fortran_int lapack_lite_status;
119
    fortran_int m;
120
    fortran_int n;
121
    fortran_int nrhs;
122
    PyObject *a;
123
    fortran_int lda;
124
    PyObject *b;
125
    fortran_int ldb;
126
    PyObject *s;
127
    double rcond;
128
    fortran_int rank;
129
    PyObject *work;
130
    PyObject *iwork;
131
    fortran_int lwork;
132
    fortran_int info;
133

134 0
    TRY(PyArg_ParseTuple(args,
135
                         (FINT_PYFMT FINT_PYFMT FINT_PYFMT "O" FINT_PYFMT "O"
136
                          FINT_PYFMT "O" "d" FINT_PYFMT "O" FINT_PYFMT "O"
137
                          FINT_PYFMT ":dgelsd"),
138
                         &m,&n,&nrhs,&a,&lda,&b,&ldb,&s,&rcond,
139
                         &rank,&work,&lwork,&iwork,&info));
140

141 0
    TRY(check_object(a,NPY_DOUBLE,"a","NPY_DOUBLE","dgelsd"));
142 0
    TRY(check_object(b,NPY_DOUBLE,"b","NPY_DOUBLE","dgelsd"));
143 0
    TRY(check_object(s,NPY_DOUBLE,"s","NPY_DOUBLE","dgelsd"));
144 0
    TRY(check_object(work,NPY_DOUBLE,"work","NPY_DOUBLE","dgelsd"));
145
#ifndef NPY_UMATH_USE_BLAS64_
146 0
    TRY(check_object(iwork,NPY_INT,"iwork","NPY_INT","dgelsd"));
147
#else
148
    TRY(check_object(iwork,NPY_INT64,"iwork","NPY_INT64","dgelsd"));
149
#endif
150

151 0
    lapack_lite_status =
152 0
            FNAME(dgelsd)(&m,&n,&nrhs,DDATA(a),&lda,DDATA(b),&ldb,
153 0
                          DDATA(s),&rcond,&rank,DDATA(work),&lwork,
154 0
                          IDATA(iwork),&info);
155 0
    if (PyErr_Occurred()) {
156
        return NULL;
157
    }
158

159 0
    return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
160
                          ",s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
161
                          ",s:d,s:" FINT_PYFMT ",s:" FINT_PYFMT
162
                          ",s:" FINT_PYFMT "}"),
163
                         "dgelsd_",lapack_lite_status,"m",m,"n",n,"nrhs",nrhs,
164
                         "lda",lda,"ldb",ldb,"rcond",rcond,"rank",rank,
165
                         "lwork",lwork,"info",info);
166
}
167

168
static PyObject *
169 1
lapack_lite_dgeqrf(PyObject *NPY_UNUSED(self), PyObject *args)
170
{
171
        fortran_int lapack_lite_status;
172
        fortran_int m, n, lwork;
173
        PyObject *a, *tau, *work;
174
        fortran_int lda;
175
        fortran_int info;
176

177 1
        TRY(PyArg_ParseTuple(args,
178
                             (FINT_PYFMT FINT_PYFMT "O" FINT_PYFMT "OO"
179
                              FINT_PYFMT FINT_PYFMT ":dgeqrf"),
180
                             &m,&n,&a,&lda,&tau,&work,&lwork,&info));
181

182
        /* check objects and convert to right storage order */
183 1
        TRY(check_object(a,NPY_DOUBLE,"a","NPY_DOUBLE","dgeqrf"));
184 1
        TRY(check_object(tau,NPY_DOUBLE,"tau","NPY_DOUBLE","dgeqrf"));
185 1
        TRY(check_object(work,NPY_DOUBLE,"work","NPY_DOUBLE","dgeqrf"));
186

187 1
        lapack_lite_status =
188 1
                FNAME(dgeqrf)(&m, &n, DDATA(a), &lda, DDATA(tau),
189 1
                              DDATA(work), &lwork, &info);
190 1
        if (PyErr_Occurred()) {
191
            return NULL;
192
        }
193

194 1
        return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
195
                              ",s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT "}"),
196
                             "dgeqrf_",
197
                             lapack_lite_status,"m",m,"n",n,"lda",lda,
198
                             "lwork",lwork,"info",info);
199
}
200

201

202
static PyObject *
203 1
lapack_lite_dorgqr(PyObject *NPY_UNUSED(self), PyObject *args)
204
{
205
        fortran_int lapack_lite_status;
206
        fortran_int m, n, k, lwork;
207
        PyObject *a, *tau, *work;
208
        fortran_int lda;
209
        fortran_int info;
210

211 1
        TRY(PyArg_ParseTuple(args,
212
                             (FINT_PYFMT FINT_PYFMT FINT_PYFMT "O"
213
                              FINT_PYFMT "OO" FINT_PYFMT FINT_PYFMT
214
                              ":dorgqr"),
215
                             &m, &n, &k, &a, &lda, &tau, &work, &lwork, &info));
216 1
        TRY(check_object(a,NPY_DOUBLE,"a","NPY_DOUBLE","dorgqr"));
217 1
        TRY(check_object(tau,NPY_DOUBLE,"tau","NPY_DOUBLE","dorgqr"));
218 1
        TRY(check_object(work,NPY_DOUBLE,"work","NPY_DOUBLE","dorgqr"));
219 1
        lapack_lite_status =
220 1
            FNAME(dorgqr)(&m, &n, &k, DDATA(a), &lda, DDATA(tau), DDATA(work),
221
                          &lwork, &info);
222 1
        if (PyErr_Occurred()) {
223
            return NULL;
224
        }
225

226 1
        return Py_BuildValue("{s:i,s:i}","dorgqr_",lapack_lite_status,
227
                             "info",info);
228
}
229

230

231
static PyObject *
232 0
lapack_lite_zgelsd(PyObject *NPY_UNUSED(self), PyObject *args)
233
{
234
    fortran_int lapack_lite_status;
235
    fortran_int m;
236
    fortran_int n;
237
    fortran_int nrhs;
238
    PyObject *a;
239
    fortran_int lda;
240
    PyObject *b;
241
    fortran_int ldb;
242
    PyObject *s;
243
    double rcond;
244
    fortran_int rank;
245
    PyObject *work;
246
    fortran_int lwork;
247
    PyObject *rwork;
248
    PyObject *iwork;
249
    fortran_int info;
250 0
    TRY(PyArg_ParseTuple(args,
251
                         (FINT_PYFMT FINT_PYFMT FINT_PYFMT "O" FINT_PYFMT
252
                          "O" FINT_PYFMT "Od" FINT_PYFMT "O" FINT_PYFMT
253
                          "OO" FINT_PYFMT ":zgelsd"),
254
                         &m,&n,&nrhs,&a,&lda,&b,&ldb,&s,&rcond,
255
                         &rank,&work,&lwork,&rwork,&iwork,&info));
256

257 0
    TRY(check_object(a,NPY_CDOUBLE,"a","NPY_CDOUBLE","zgelsd"));
258 0
    TRY(check_object(b,NPY_CDOUBLE,"b","NPY_CDOUBLE","zgelsd"));
259 0
    TRY(check_object(s,NPY_DOUBLE,"s","NPY_DOUBLE","zgelsd"));
260 0
    TRY(check_object(work,NPY_CDOUBLE,"work","NPY_CDOUBLE","zgelsd"));
261 0
    TRY(check_object(rwork,NPY_DOUBLE,"rwork","NPY_DOUBLE","zgelsd"));
262
#ifndef NPY_UMATH_USE_BLAS64_
263 0
    TRY(check_object(iwork,NPY_INT,"iwork","NPY_INT","zgelsd"));
264
#else
265
    TRY(check_object(iwork,NPY_INT64,"iwork","NPY_INT64","zgelsd"));
266
#endif
267

268 0
    lapack_lite_status =
269 0
        FNAME(zgelsd)(&m,&n,&nrhs,ZDATA(a),&lda,ZDATA(b),&ldb,DDATA(s),&rcond,
270 0
                      &rank,ZDATA(work),&lwork,DDATA(rwork),IDATA(iwork),&info);
271 0
    if (PyErr_Occurred()) {
272
        return NULL;
273
    }
274

275 0
    return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
276
                          ",s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
277
                          ",s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
278
                          "}"),
279
                         "zgelsd_",
280
                         lapack_lite_status,"m",m,"n",n,"nrhs",nrhs,"lda",lda,
281
                         "ldb",ldb,"rank",rank,"lwork",lwork,"info",info);
282
}
283

284
static PyObject *
285 1
lapack_lite_zgeqrf(PyObject *NPY_UNUSED(self), PyObject *args)
286
{
287
        fortran_int lapack_lite_status;
288
        fortran_int m, n, lwork;
289
        PyObject *a, *tau, *work;
290
        fortran_int lda;
291
        fortran_int info;
292

293 1
        TRY(PyArg_ParseTuple(args,
294
                             (FINT_PYFMT FINT_PYFMT "O" FINT_PYFMT "OO"
295
                              FINT_PYFMT "" FINT_PYFMT ":zgeqrf"),
296
                             &m,&n,&a,&lda,&tau,&work,&lwork,&info));
297

298
/* check objects and convert to right storage order */
299 1
        TRY(check_object(a,NPY_CDOUBLE,"a","NPY_CDOUBLE","zgeqrf"));
300 1
        TRY(check_object(tau,NPY_CDOUBLE,"tau","NPY_CDOUBLE","zgeqrf"));
301 1
        TRY(check_object(work,NPY_CDOUBLE,"work","NPY_CDOUBLE","zgeqrf"));
302

303 1
        lapack_lite_status =
304 1
            FNAME(zgeqrf)(&m, &n, ZDATA(a), &lda, ZDATA(tau), ZDATA(work),
305
                          &lwork, &info);
306 1
        if (PyErr_Occurred()) {
307
            return NULL;
308
        }
309

310 1
        return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT
311
                              ",s:" FINT_PYFMT ",s:" FINT_PYFMT
312
                              ",s:" FINT_PYFMT ",s:" FINT_PYFMT "}"),
313
                             "zgeqrf_",lapack_lite_status,"m",m,"n",n,"lda",lda,"lwork",lwork,"info",info);
314
}
315

316

317
static PyObject *
318 1
lapack_lite_zungqr(PyObject *NPY_UNUSED(self), PyObject *args)
319
{
320
        fortran_int lapack_lite_status;
321
        fortran_int m, n, k, lwork;
322
        PyObject *a, *tau, *work;
323
        fortran_int lda;
324
        fortran_int info;
325

326 1
        TRY(PyArg_ParseTuple(args,
327
                             (FINT_PYFMT FINT_PYFMT FINT_PYFMT "O"
328
                              FINT_PYFMT "OO" FINT_PYFMT "" FINT_PYFMT
329
                              ":zungqr"),
330
                             &m, &n, &k, &a, &lda, &tau, &work, &lwork, &info));
331 1
        TRY(check_object(a,NPY_CDOUBLE,"a","NPY_CDOUBLE","zungqr"));
332 1
        TRY(check_object(tau,NPY_CDOUBLE,"tau","NPY_CDOUBLE","zungqr"));
333 1
        TRY(check_object(work,NPY_CDOUBLE,"work","NPY_CDOUBLE","zungqr"));
334

335

336 1
        lapack_lite_status =
337 1
            FNAME(zungqr)(&m, &n, &k, ZDATA(a), &lda, ZDATA(tau), ZDATA(work),
338
                          &lwork, &info);
339 1
        if (PyErr_Occurred()) {
340
            return NULL;
341
        }
342

343 1
        return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT "}"),
344
                             "zungqr_",lapack_lite_status,
345
                             "info",info);
346
}
347

348

349
static PyObject *
350 1
lapack_lite_xerbla(PyObject *NPY_UNUSED(self), PyObject *args)
351
{
352 1
    fortran_int info = -1;
353

354 1
    NPY_BEGIN_THREADS_DEF;
355 1
    NPY_BEGIN_THREADS;
356 1
    FNAME(xerbla)("test", &info);
357 1
    NPY_END_THREADS;
358

359 1
    if (PyErr_Occurred()) {
360
        return NULL;
361
    }
362 0
    Py_RETURN_NONE;
363
}
364

365

366
#define STR(x) #x
367
#define lameth(name) {STR(name), lapack_lite_##name, METH_VARARGS, NULL}
368
static struct PyMethodDef lapack_lite_module_methods[] = {
369
    lameth(dgelsd),
370
    lameth(dgeqrf),
371
    lameth(dorgqr),
372
    lameth(zgelsd),
373
    lameth(zgeqrf),
374
    lameth(zungqr),
375
    lameth(xerbla),
376
    { NULL,NULL,0, NULL}
377
};
378

379

380
static struct PyModuleDef moduledef = {
381
        PyModuleDef_HEAD_INIT,
382
        "lapack_lite",
383
        NULL,
384
        -1,
385
        lapack_lite_module_methods,
386
        NULL,
387
        NULL,
388
        NULL,
389
        NULL
390
};
391

392
/* Initialization function for the module */
393 1
PyMODINIT_FUNC PyInit_lapack_lite(void)
394
{
395
    PyObject *m,*d;
396 1
    m = PyModule_Create(&moduledef);
397 1
    if (m == NULL) {
398
        return NULL;
399
    }
400 1
    import_array();
401 1
    d = PyModule_GetDict(m);
402 1
    LapackError = PyErr_NewException("lapack_lite.LapackError", NULL, NULL);
403 1
    PyDict_SetItemString(d, "LapackError", LapackError);
404

405
#ifdef HAVE_BLAS_ILP64
406
    PyDict_SetItemString(d, "_ilp64", Py_True);
407
#else
408 1
    PyDict_SetItemString(d, "_ilp64", Py_False);
409
#endif
410

411 1
    return m;
412
}

Read our documentation on viewing source code .

Loading