Skip to content

Commit

Permalink
Fixes #3
Browse files Browse the repository at this point in the history 
Adding CHOL(finish) when "common" struct is not needed anymore
  • Loading branch information
@sanurielf
sanurielf committed Nov 10, 2021
1 parent cf2de59 commit fc0c30a
Show file tree
Hide file tree
Showing 4 changed files with 81 additions and 61 deletions.
4 changes: 2 additions & 2 deletions setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
48 changes: 34 additions & 14 deletions src/C/cholmod.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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();

Expand All @@ -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));
Expand Down Expand Up @@ -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,
Expand All @@ -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->minor<n) PY_ERR(PyExc_ArithmeticError, "singular matrix");
if ((int_t) L->minor<n) PY_ERR(PyExc_ArithmeticError, "singular matrix");

if (sys < 0 || sys > 8)
PY_ERR(PyExc_ValueError, "invalid value for sys");
Expand Down Expand Up @@ -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("");
}
Expand Down Expand Up @@ -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,
Expand All @@ -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->minor<n) PY_ERR(PyExc_ArithmeticError, "singular matrix");
if ((int_t) L->minor<n) PY_ERR(PyExc_ArithmeticError, "singular matrix");

if (sys < 0 || sys > 8)
PY_ERR(PyExc_ValueError, "invalid value for sys");
Expand All @@ -566,13 +574,15 @@ 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));
memcpy(SP_ROW(X), Xc->i, ((int_t *)Xc->p)[Xc->ncol]*sizeof(int_t));
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;
}

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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 {
Expand All @@ -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;
}
}
Expand All @@ -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)
Expand All @@ -697,10 +712,12 @@ static PyObject* linsolve(PyObject *self, PyObject *args,
break;

default:
CHOL(finish)(&Common);
PyErr_Warn(PyExc_UserWarning, "");
}

if (L->minor<n) {
if ((int_t) L->minor<n) {
CHOL(finish)(&Common);
CHOL(free_factor)(&L, &Common);
PY_ERR(PyExc_ArithmeticError, "singular matrix");
}
Expand All @@ -709,6 +726,7 @@ static PyObject* linsolve(PyObject *self, PyObject *args,
if (Common.status == CHOLMOD_OUT_OF_MEMORY) {
CHOL(free_factor)(&L, &Common);
CHOL(free_dense)(&b, &Common);
CHOL(finish)(&Common);
return PyErr_NoMemory();
}
b_old = b->x;
Expand All @@ -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)]);
Expand All @@ -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("");
}

Expand Down Expand Up @@ -826,7 +846,7 @@ static PyObject* splinsolve(PyObject *self, PyObject *args,
PyErr_Warn(PyExc_UserWarning, "");
}

if (L->minor<n) {
if ((int_t) L->minor<n) {
CHOL(free_factor)(&L, &Common);
PY_ERR(PyExc_ArithmeticError, "singular matrix");
}
Expand Down Expand Up @@ -883,8 +903,8 @@ static PyObject* diag(PyObject *self, PyObject *args)
matrix *d=NULL;
cholmod_factor *L;
const char *descr;
int k, strt, incx=1, incy, nrows, ncols;

int strt, incx=1, incy, nrows, ncols;
int_t k;
if (!set_options()) return NULL;
if (!PyArg_ParseTuple(args, "O", &F)) return NULL;

Expand All @@ -904,7 +924,7 @@ static PyObject* diag(PyObject *self, PyObject *args)
COMPLEX))) return NULL;

strt = 0;
for (k=0; k<L->nsuper; 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] */
Expand Down
42 changes: 21 additions & 21 deletions src/C/klu.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand All @@ -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),
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -454,20 +454,20 @@ 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),
Pt, Qt, SP_VALD(R), rt, &Common);
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));
Expand All @@ -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,
Expand All @@ -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",
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit fc0c30a

Please sign in to comment.