diff --git a/setup.py b/setup.py index a57a75aa..5e5cfe2f 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,7 @@ # @Author: Uriel Sandoval # @Date: 2021-09-20 08:53:47 # @Last Modified by: Uriel Sandoval -# @Last Modified time: 2021-09-20 10:24:22 +# @Last Modified time: 2021-11-03 22:44:06 from setuptools import setup, Extension from glob import glob import os, sys @@ -14,7 +14,7 @@ # Default names of BLAS and LAPACK libraries BLAS_LIB = ['blas'] LAPACK_LIB = ['lapack'] -FFTW_LIB = ['fftw3'] +FFTW_LIB = ['fftw3'] BLAS_EXTRA_LINK_ARGS = [] # Set environment variable BLAS_NOUNDERSCORES=1 if your BLAS/LAPACK do diff --git a/src/C/cholmod.c b/src/C/cholmod.c index 78d10f83..58ffac07 100644 --- a/src/C/cholmod.c +++ b/src/C/cholmod.c @@ -283,6 +283,7 @@ static PyObject* symbolic(PyObject *self, PyObject *args, return NULL; } } + CHOL(finish)(&Common); return (PyObject *) PyCapsule_New((void *) L, SP_ID(A)==DOUBLE ? (uplo == 'L' ? descrdFs_L : descrdFs_U) : (uplo == 'L' ? descrzFs_L : descrzFs_U), @@ -324,6 +325,7 @@ static PyObject* numeric(PyObject *self, PyObject *args) PyObject *F; cholmod_factor *Lc; cholmod_sparse *Ac = NULL; + int status; char uplo; const char *descr; @@ -356,10 +358,12 @@ static PyObject* numeric(PyObject *self, PyObject *args) Lc = (cholmod_factor *) PyCapsule_GetPointer(F, descr); if (!(Ac = pack(A, uplo))) return PyErr_NoMemory(); + status = Common.status; CHOL(factorize) (Ac, Lc, &Common); CHOL(free_sparse)(&Ac, &Common); + CHOL(finish)(&Common); - if (Common.status < 0) switch (Common.status) { + if (status < 0) switch (status) { case CHOLMOD_OUT_OF_MEMORY: return PyErr_NoMemory(); @@ -368,7 +372,7 @@ static PyObject* numeric(PyObject *self, PyObject *args) return NULL; } - if (Common.status > 0) switch (Common.status) { + if (status > 0) switch (status) { case CHOLMOD_NOT_POSDEF: PyErr_SetObject(PyExc_ArithmeticError, Py_BuildValue("i", Lc->minor)); @@ -426,7 +430,8 @@ static PyObject* solve(PyObject *self, PyObject *args, PyObject *kwrds) { matrix *B; PyObject *F; - int i, n, oB=0, ldB=0, nrhs=-1, sys=0; + int i, oB=0, ldB=0, nrhs=-1, sys=0; + int_t n; const char *descr; char *kwlist[] = {"F", "B", "sys", "nrhs", "ldB", "offsetB", NULL}; int sysvalues[] = { CHOLMOD_A, CHOLMOD_LDLt, CHOLMOD_LD, @@ -448,7 +453,7 @@ static PyObject* solve(PyObject *self, PyObject *args, PyObject *kwrds) PY_ERR(PyExc_ValueError, "called with symbolic factor"); n = L->n; - if (L->minorminor 8) PY_ERR(PyExc_ValueError, "invalid value for sys"); @@ -480,13 +485,15 @@ static PyObject* solve(PyObject *self, PyObject *args, PyObject *kwrds) PyErr_SetString(PyExc_ValueError, "solve step failed"); CHOL(free_dense)(&x, &Common); CHOL(free_dense)(&b, &Common); - return NULL; + CHOL(finish)(&Common); + return NULL; } memcpy(b->x, x->x, n*E_SIZE[MAT_ID(B)]); CHOL(free_dense)(&x, &Common); } b->x = b_old; CHOL(free_dense)(&b, &Common); + CHOL(finish)(&Common); return Py_BuildValue(""); } @@ -521,7 +528,8 @@ static PyObject* spsolve(PyObject *self, PyObject *args, cholmod_sparse *Bc=NULL, *Xc=NULL; PyObject *F; cholmod_factor *L; - int n, sys=0; + int sys=0; + int_t n; const char *descr; char *kwlist[] = {"F", "B", "sys", NULL}; int sysvalues[] = {CHOLMOD_A, CHOLMOD_LDLt, CHOLMOD_LD, @@ -542,7 +550,7 @@ static PyObject* spsolve(PyObject *self, PyObject *args, if (L->xtype == CHOLMOD_PATTERN) PY_ERR(PyExc_ValueError, "called with symbolic factor"); n = L->n; - if (L->minorminor 8) PY_ERR(PyExc_ValueError, "invalid value for sys"); @@ -566,6 +574,7 @@ static PyObject* spsolve(PyObject *self, PyObject *args, ((int_t*)Xc->p)[Xc->ncol], (L->xtype == CHOLMOD_REAL ? DOUBLE : COMPLEX)))) { CHOL(free_sparse)(&Xc, &Common); + CHOL(finish)(&Common); return NULL; } memcpy(SP_COL(X), Xc->p, (Xc->ncol+1)*sizeof(int_t)); @@ -573,6 +582,7 @@ static PyObject* spsolve(PyObject *self, PyObject *args, memcpy(SP_VAL(X), Xc->x, ((int_t *) Xc->p)[Xc->ncol]*E_SIZE[SP_ID(X)]); CHOL(free_sparse)(&Xc, &Common); + CHOL(finish)(&Common); return (PyObject *) X; } @@ -610,7 +620,8 @@ static PyObject* linsolve(PyObject *self, PyObject *args, { spmatrix *A; matrix *B, *P=NULL; - int i, n, oB=0, ldB=0, nrhs=-1; + int i, oB=0, ldB=0, nrhs=-1; + int_t n; cholmod_sparse *Ac=NULL; cholmod_factor *L=NULL; cholmod_dense *x=NULL, *b=NULL; @@ -654,6 +665,7 @@ static PyObject* linsolve(PyObject *self, PyObject *args, free_matrix(Ac); CHOL(free_sparse)(&Ac, &Common); CHOL(free_factor)(&L, &Common); + CHOL(finish)(&Common); if (Common.status == CHOLMOD_OUT_OF_MEMORY) return PyErr_NoMemory(); else { @@ -662,18 +674,19 @@ static PyObject* linsolve(PyObject *self, PyObject *args, return NULL; } } - CHOL(factorize) (Ac, L, &Common); CHOL(free_sparse)(&Ac, &Common); if (Common.status < 0) { CHOL(free_factor)(&L, &Common); switch (Common.status) { case CHOLMOD_OUT_OF_MEMORY: + CHOL(finish)(&Common); return PyErr_NoMemory(); default: PyErr_SetString(PyExc_ValueError, "factorization " "failed"); + CHOL(finish)(&Common); return NULL; } } @@ -682,10 +695,12 @@ static PyObject* linsolve(PyObject *self, PyObject *args, PyErr_SetObject(PyExc_ArithmeticError, Py_BuildValue("i", L->minor)); CHOL(free_factor)(&L, &Common); + CHOL(finish)(&Common); return NULL; break; case CHOLMOD_DSMALL: + CHOL(finish)(&Common); /* This never happens unless we change the default value * of Common.dbound (0.0). */ if (L->is_ll) @@ -697,10 +712,12 @@ static PyObject* linsolve(PyObject *self, PyObject *args, break; default: + CHOL(finish)(&Common); PyErr_Warn(PyExc_UserWarning, ""); } - if (L->minorminorx; @@ -721,6 +739,7 @@ static PyObject* linsolve(PyObject *self, PyObject *args, b->x = b_old; CHOL(free_dense)(&b, &Common); CHOL(free_dense)(&x, &Common); + CHOL(finish)(&Common); return NULL; } memcpy(b->x, x->x, SP_NROWS(A)*E_SIZE[MAT_ID(B)]); @@ -729,6 +748,7 @@ static PyObject* linsolve(PyObject *self, PyObject *args, b->x = b_old; CHOL(free_dense)(&b, &Common); CHOL(free_factor)(&L, &Common); + CHOL(finish)(&Common); return Py_BuildValue(""); } @@ -826,7 +846,7 @@ static PyObject* splinsolve(PyObject *self, PyObject *args, PyErr_Warn(PyExc_UserWarning, ""); } - if (L->minorminornsuper; k++){ + for (k=0; k<(int_t) L->nsuper; k++){ /* x[L->px[k], .... ,L->px[k+1]-1] is a dense lower-triangular * nrowx times ncols matrix. We copy its diagonal to * d[strt, ..., strt+ncols-1] */ diff --git a/src/C/klu.c b/src/C/klu.c index b4b7fb5c..d67c4074 100644 --- a/src/C/klu.c +++ b/src/C/klu.c @@ -176,7 +176,7 @@ static PyObject* linsolve(PyObject *self, PyObject *args, if (Common.status == KLU_SINGULAR) PyErr_SetString(PyExc_ArithmeticError, "singular matrix"); else { - snprintf(klu_error, 20, "%s %i", "KLU ERROR", + snprintf(klu_error, 20, "%s %i", "KLU ERROR", (int) Common.status); PyErr_SetString(PyExc_ValueError, klu_error); } @@ -193,7 +193,7 @@ static PyObject* linsolve(PyObject *self, PyObject *args, } else { if (trans == 'N') - KLUZ(solve)(Symbolic, Numeric, n, nrhs, (double *) MAT_BUFZ(B), + KLUZ(solve)(Symbolic, Numeric, n, nrhs, (double *) MAT_BUFZ(B), &Common); else KLUZ(tsolve)(Symbolic, Numeric, n, nrhs, (double *) MAT_BUFZ(B), @@ -265,7 +265,7 @@ static PyObject* symbolic(PyObject *self, PyObject *args) Symbolic = KLUD(analyze)(n, SP_COL(A), SP_ROW(A), &Common); if (Common.status == KLU_OK) { /* Symbolic is the same for DOUBLE and COMPLEX cases. Only make - * difference in the Capsule descriptor to avoid user errors + * difference in the Capsule descriptor to avoid user errors */ if (SP_ID(A) == DOUBLE) return (PyObject *) PyCapsule_New( @@ -401,12 +401,12 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) KLU(common) Common; double *Lx, *Lz, *Ux, *Uz, *Fx, *Fz; - if (!PyArg_ParseTuple(args, "OOO", &A, &Fs, &Fn)) + if (!PyArg_ParseTuple(args, "OOO", &A, &Fs, &Fn)) return NULL; if (!SpMatrix_Check(A)) PY_ERR_TYPE("A must be a sparse matrix"); - + nn = SP_NROWS(A); @@ -454,12 +454,12 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) /* We store r blocks as c array, then we can transform it into Python list */ rt = malloc((symbolic->nblocks+1) * sizeof(int_t)); - /* KLU_extract does not handle packed complex arrays, thus we create + /* KLU_extract does not handle packed complex arrays, thus we create * two auxiliary arrays for L, U and F */ switch (SP_ID(A)) { case DOUBLE: - status = KLUD(extract)(numeric, symbolic, + status = KLUD(extract)(numeric, symbolic, SP_COL(L), SP_ROW(L), SP_VALD(L), SP_COL(U), SP_ROW(U), SP_VALD(U), SP_COL(F), SP_ROW(F), SP_VALD(F), @@ -467,7 +467,7 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) break; case COMPLEX: - /* KLU_extract does not handle packed complex arrays, thus we create + /* KLU_extract does not handle packed complex arrays, thus we create * two auxiliary arrays for L, U and F */ Lx = malloc(lnz * sizeof(double)); @@ -477,7 +477,7 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) Fx = malloc(fnz * sizeof(double)); Fz = malloc(fnz * sizeof(double)); - status = KLUZ(extract)(numeric, symbolic, + status = KLUZ(extract)(numeric, symbolic, SP_COL(L), SP_ROW(L), Lx, Lz, SP_COL(U), SP_ROW(U), Ux, Uz, SP_COL(F), SP_ROW(F), Fx, Fz, @@ -504,7 +504,7 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) break; - } + } if (status != 1){ snprintf(klu_error, 20, "%s %i", "KLU ERROR", @@ -554,13 +554,13 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) free(Qt); if (!(res = PyTuple_New(7))) return PyErr_NoMemory(); - PyTuple_SET_ITEM(res, 0, L); - PyTuple_SET_ITEM(res, 1, U); - PyTuple_SET_ITEM(res, 2, P); - PyTuple_SET_ITEM(res, 3, Q); - PyTuple_SET_ITEM(res, 4, R); - PyTuple_SET_ITEM(res, 5, F); - PyTuple_SET_ITEM(res, 6, r); + PyTuple_SET_ITEM(res, 0, (PyObject *) L); + PyTuple_SET_ITEM(res, 1, (PyObject *) U); + PyTuple_SET_ITEM(res, 2, (PyObject *) P); + PyTuple_SET_ITEM(res, 3, (PyObject *) Q); + PyTuple_SET_ITEM(res, 4, (PyObject *) R); + PyTuple_SET_ITEM(res, 5, (PyObject *) F); + PyTuple_SET_ITEM(res, 6, (PyObject *) r); return res; @@ -691,7 +691,7 @@ static PyObject* solve(PyObject *self, PyObject *args, PyObject *kwrds) } -static char doc_get_det[] = +static char doc_get_det[] = "Returns determinant of a KLU symbolic/numeric object\n" "d = get_det(A, Fs, Fn)\n\n" "PURPOSE\n" @@ -719,10 +719,10 @@ static PyObject* get_det(PyObject *self, PyObject *args, PyObject *kwrds) { double det = 1, d_sign; #ifndef _MSC_VER double complex det_c = 1.0; - double complex *Uzdiag; + double complex *Uzdiag = NULL; #else _Dcomplex det_c = _Cbuild(1.0, 0.0); - _Dcomplex *Uzdiag; + _Dcomplex *Uzdiag = NULL; #endif int_t i, k, n, npiv, itmp, *P, *Q, *Wi; @@ -767,7 +767,7 @@ static PyObject* get_det(PyObject *self, PyObject *args, PyObject *kwrds) { Q = Fsptr->Q; Rs = Fptr->Rs; - if (SP_ID(A) == DOUBLE) + if (SP_ID(A) == DOUBLE) for (k = 0; k < n; k++) det *= (Udiag[k]*Rs[k]); else diff --git a/src/C/umfpack.c b/src/C/umfpack.c index 8df3073e..826de507 100644 --- a/src/C/umfpack.c +++ b/src/C/umfpack.c @@ -42,8 +42,8 @@ static char umfpack_error[40]; const char *descrdFs = "UMFPACK SYM D FACTOR"; const char *descrzFs = "UMFPACK SYM Z FACTOR"; -const char *descrdFn = "UMFPACK NUM D FACTOR"; -const char *descrzFn = "UMFPACK NUM Z FACTOR"; +const char *descrdFn = "UMFPACK NUM D FACTOR"; +const char *descrzFn = "UMFPACK NUM Z FACTOR"; PyDoc_STRVAR(umfpack__doc__,"Interface to the UMFPACK library.\n\n" "Routines for symbolic and numeric LU factorization of sparse\n" @@ -257,9 +257,9 @@ static PyObject* symbolic(PyObject *self, PyObject *args) UMFD(symbolic)(SP_NROWS(A), SP_NCOLS(A), SP_COL(A), SP_ROW(A), SP_VAL(A), &symbolic, NULL, info); if (info[UMFPACK_STATUS] == UMFPACK_OK) - return (PyObject *) PyCapsule_New( (void *) symbolic, - descrdFs, - (PyCapsule_Destructor) &free_umfpack_d_symbolic); + return (PyObject *) PyCapsule_New( (void *) symbolic, + descrdFs, + (PyCapsule_Destructor) &free_umfpack_d_symbolic); else UMFD(free_symbolic)(&symbolic); @@ -339,7 +339,7 @@ static PyObject* numeric(PyObject *self, PyObject *args) "factor of a 'z' matrix"); if (!(Fsptr = (void *) PyCapsule_GetPointer(Fs, descrzFs))) err_CO("Fs"); - UMFZ(numeric)(SP_COL(A), SP_ROW(A), SP_VAL(A), NULL, Fsptr, + UMFZ(numeric)(SP_COL(A), SP_ROW(A), SP_VAL(A), NULL, Fsptr, &numeric, NULL, info); if (info[UMFPACK_STATUS] == UMFPACK_OK) @@ -382,7 +382,7 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) void *numeric; int trans_ = 'N'; char trans='N'; - int_t i, n_row, n_col, nn, n_inner, lnz, unz, status, ignore1, ignore2, ignore3, + int_t i, n_row, n_col, nn, n_inner, lnz, unz, status, ignore1, ignore2, ignore3, *Ltp, *Ltj, *Up, *Ui, *Pt, *Qt, *Rp, *Ri, do_recip; double *Ltx, *Ux, *Rs; char *kwlist[] = {"A", "F", "trans", NULL}; @@ -392,7 +392,7 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) trans = (char) trans_; if (!SpMatrix_Check(A)) PY_ERR_TYPE("A must be a sparse matrix"); - + n_row = SP_NROWS(A); n_col = SP_NCOLS(A); @@ -412,13 +412,13 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) switch (SP_ID(A)) { case DOUBLE: numeric = (void *) PyCapsule_GetPointer(F, descrdFn); - status = UMFD(get_lunz)(&lnz, &unz, &ignore1, &ignore2, &ignore3, + status = UMFD(get_lunz)(&lnz, &unz, &ignore1, &ignore2, &ignore3, numeric); break; case COMPLEX: numeric = (void *) PyCapsule_GetPointer(F, descrzFn); - status = UMFZ(get_lunz)(&lnz, &unz, &ignore1, &ignore2, &ignore3, + status = UMFZ(get_lunz)(&lnz, &unz, &ignore1, &ignore2, &ignore3, numeric); break; @@ -471,7 +471,7 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) case COMPLEX: - status = UMFZ(get_numeric)(Ltp, Ltj, Ltx, NULL, Up, Ui, Ux, NULL, Pt, + status = UMFZ(get_numeric)(Ltp, Ltj, Ltx, NULL, Up, Ui, Ux, NULL, Pt, Qt, NULL, NULL, &do_recip, Rs, numeric); break; } @@ -519,13 +519,13 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) /* convert L from row form to column form */ switch (SP_ID(A)) { case DOUBLE: - status = UMFD(transpose)(n_inner, n_row, Ltp, Ltj, Ltx, NULL, NULL, + status = UMFD(transpose)(n_inner, n_row, Ltp, Ltj, Ltx, NULL, NULL, SP_COL(L), SP_ROW(L), SP_VALD(L)); break; case COMPLEX: - status = UMFZ(transpose)(n_inner, n_row, Ltp, Ltj, Ltx, NULL, NULL, - NULL, SP_COL(L), SP_ROW(L), SP_VALD(L), NULL, 0); - break; + status = UMFZ(transpose)(n_inner, n_row, Ltp, Ltj, Ltx, NULL, NULL, + NULL, SP_COL(L), SP_ROW(L), SP_VALD(L), NULL, 0); + break; } @@ -545,15 +545,15 @@ static PyObject* get_numeric(PyObject *self, PyObject *args, PyObject *kwrds) } if (!(res = PyTuple_New(5))) return PyErr_NoMemory(); - PyTuple_SET_ITEM(res, 0, L); - PyTuple_SET_ITEM(res, 1, U); - PyTuple_SET_ITEM(res, 2, P); - PyTuple_SET_ITEM(res, 3, Q); - PyTuple_SET_ITEM(res, 4, R); + PyTuple_SET_ITEM(res, 0, (PyObject *) L); + PyTuple_SET_ITEM(res, 1, (PyObject *) U); + PyTuple_SET_ITEM(res, 2, (PyObject *) P); + PyTuple_SET_ITEM(res, 3, (PyObject *) Q); + PyTuple_SET_ITEM(res, 4, (PyObject *) R); return res; - + } static char doc_solve[] = @@ -668,7 +668,7 @@ static PyObject* solve(PyObject *self, PyObject *args, PyObject *kwrds) -static char doc_get_det[] = +static char doc_get_det[] = "Returns determinant of a UMFPACK symbolic/numeric object\n" "d = get_det(A, Fs, Fn)\n\n" "PURPOSE\n" @@ -698,7 +698,7 @@ static PyObject* get_det(PyObject *self, PyObject *args, PyObject *kwrds) { if (SP_ID(A) == DOUBLE){ TypeCheck_Capsule(Fn, descrdFn, "F is not the UMFPACK numeric factor " "of a 'd' matrix"); - status = UMFD(get_determinant)(&dx, NULL, + status = UMFD(get_determinant)(&dx, NULL, (void *) PyCapsule_GetPointer(Fn, descrdFn), info); } @@ -722,7 +722,7 @@ static PyObject* get_det(PyObject *self, PyObject *args, PyObject *kwrds) { return Py_BuildValue("d", dx); else return PyComplex_FromDoubles(dx, dz); - + } static PyMethodDef umfpack_functions[] = {