Skip to content

implement :det_exact and :inverse_exact for list/yale matrices #334

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

implement :det_exact and :inverse_exact for list/yale matrices #334

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
216 changes: 198 additions & 18 deletions ext/nmatrix/math.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,15 +207,14 @@ namespace nm {
* Calculate the determinant for a dense matrix (A [elements]) of size 2 or 3. Return the result.
*/
template <typename DType>
void det_exact(const int M, const void* A_elements, const int lda, void* result_arg) {
void det_exact_from_dense(const int M, const void* A_elements, const int lda, void* result_arg) {
DType* result = reinterpret_cast<DType*>(result_arg);
const DType* A = reinterpret_cast<const DType*>(A_elements);

typename LongDType<DType>::type x, y;

if (M == 2) {
*result = A[0] * A[lda+1] - A[1] * A[lda];

} else if (M == 3) {
x = A[lda+1] * A[2*lda+2] - A[lda+2] * A[2*lda+1]; // ei - fh
y = A[lda] * A[2*lda+2] - A[lda+2] * A[2*lda]; // fg - di
Expand All @@ -230,6 +229,53 @@ namespace nm {
}
}

/*
* Calculate the determinant for a yale matrix (storage) of size 2 or 3. Return the result.
*/
template <typename DType>
void det_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, void* result_arg) {
DType* result = reinterpret_cast<DType*>(result_arg);
IType* ija = reinterpret_cast<IType *>(storage->ija);
DType* a = reinterpret_cast<DType*>(storage->a);
IType col_pos = storage->shape[0] + 1;
if (M == 2) {
if (ija[2] - ija[0] == 2) {
*result = a[0] * a[1] - a[col_pos] * a[col_pos+1];
}
else { *result = a[0] * a[1]; }
} else if (M == 3) {
DType m[3][3];
for (int i = 0; i < 3; ++i) {
m[i][i] = a[i];
switch(ija[i+1] - ija[i]) {
case 2:
m[i][ija[col_pos]] = a[col_pos];
m[i][ija[col_pos+1]] = a[col_pos+1];
col_pos += 2;
break;
case 1:
m[i][(i+1)%3] = m[i][(i+2)%3] = 0;
m[i][ija[col_pos]] = a[col_pos];
++col_pos;
break;
case 0:
m[i][(i+1)%3] = m[i][(i+2)%3] = 0;
break;
default:
rb_raise(rb_eArgError, "some value in IJA is incorrect!");
}
}
*result =
m[0][0] * m[1][1] * m[2][2] + m[0][1] * m[1][2] * m[2][0] + m[0][2] * m[1][0] * m[2][1]
- m[0][0] * m[1][2] * m[2][1] - m[0][1] * m[1][0] * m[2][2] - m[0][2] * m[1][1] * m[2][0];

} else if (M < 2) {
rb_raise(rb_eArgError, "can only calculate exact determinant of a square matrix of size 2 or larger");
} else {
rb_raise(rb_eNotImpError, "exact determinant calculation needed for matrices larger than 3x3");
}
}

/*
* Solve a system of linear equations using forward-substution followed by
* back substution from the LU factorization of the matrix of co-efficients.
Expand Down Expand Up @@ -425,20 +471,24 @@ namespace nm {
delete[] u;
}

void raise_not_invertible_error() {
rb_raise(nm_eNotInvertibleError,
"matrix must have non-zero determinant to be invertible (not getting this error does not mean matrix is invertible if you're dealing with floating points)");
}

/*
* Calculate the exact inverse for a dense matrix (A [elements]) of size 2 or 3. Places the result in B_elements.
*/
template <typename DType>
void inverse_exact(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb) {
void inverse_exact_from_dense(const int M, const void* A_elements,
const int lda, void* B_elements, const int ldb) {

const DType* A = reinterpret_cast<const DType*>(A_elements);
DType* B = reinterpret_cast<DType*>(B_elements);

if (M == 2) {
DType det = A[0] * A[lda+1] - A[1] * A[lda];
if (det == 0) {
rb_raise(nm_eNotInvertibleError,
"matrix must have non-zero determinant to be invertible (not getting this error does not mean matrix is invertible if you're dealing with floating points)");
}
if (det == 0) { raise_not_invertible_error(); }
B[0] = A[lda+1] / det;
B[1] = -A[1] / det;
B[ldb] = -A[lda] / det;
Expand All @@ -447,11 +497,8 @@ namespace nm {
} else if (M == 3) {
// Calculate the exact determinant.
DType det;
det_exact<DType>(M, A_elements, lda, reinterpret_cast<void*>(&det));
if (det == 0) {
rb_raise(nm_eNotInvertibleError,
"matrix must have non-zero determinant to be invertible (not getting this error does not mean matrix is invertible if you're dealing with floating points)");
}
det_exact_from_dense<DType>(M, A_elements, lda, reinterpret_cast<void*>(&det));
if (det == 0) { raise_not_invertible_error(); }

B[0] = ( A[lda+1] * A[2*lda+2] - A[lda+2] * A[2*lda+1]) / det; // A = ei - fh
B[1] = (- A[1] * A[2*lda+2] + A[2] * A[2*lda+1]) / det; // D = -bi + ch
Expand All @@ -469,6 +516,111 @@ namespace nm {
}
}

template <typename DType>
void inverse_exact_from_yale(const int M, const YALE_STORAGE* storage,
const int lda, YALE_STORAGE* inverse, const int ldb) {

// inverse is a clone of storage
const DType* a = reinterpret_cast<const DType*>(storage->a);
const IType* ija = reinterpret_cast<const IType *>(storage->ija);
DType* b = reinterpret_cast<DType*>(inverse->a);
IType* ijb = reinterpret_cast<IType *>(inverse->ija);
IType col_pos = storage->shape[0] + 1;
// Calculate the exact determinant.
DType det;

if (M == 2) {
IType ndnz = ija[2] - ija[0];
if (ndnz == 2) {
det = a[0] * a[1] - a[col_pos] * a[col_pos+1];
}
else { det = a[0] * a[1]; }
if (det == 0) { raise_not_invertible_error(); }
b[0] = a[1] / det;
b[1] = -a[0] / det;
if (ndnz == 2) {
b[col_pos] = -a[col_pos] / det;
b[col_pos+1] = -a[col_pos+1] / det;
}
else if (ndnz == 1) {
b[col_pos] = -a[col_pos] / det;
}

} else if (M == 3) {
DType A[lda*3];
for (int i = 0; i < lda; ++i) {
A[i*3+i] = a[i];
switch (ija[i+1] - ija[i]) {
case 2:
A[i*3 + ija[col_pos]] = a[col_pos];
A[i*3 + ija[col_pos+1]] = a[col_pos+1];
col_pos += 2;
break;
case 1:
A[i*3 + (i+1)%3] = A[i*3 + (i+2)%3] = 0;
A[i*3 + ija[col_pos]] = a[col_pos];
col_pos += 1;
break;
case 0:
A[i*3 + (i+1)%3] = A[i*3 + (i+2)%3] = 0;
break;
default:
rb_raise(rb_eArgError, "some value in IJA is incorrect!");
}
}
det =
A[0] * A[lda+1] * A[2*lda+2] + A[1] * A[lda+2] * A[2*lda] + A[2] * A[lda] * A[2*lda+1]
- A[0] * A[lda+2] * A[2*lda+1] - A[1] * A[lda] * A[2*lda+2] - A[2] * A[lda+1] * A[2*lda];
if (det == 0) { raise_not_invertible_error(); }

DType B[3*ldb];
B[0] = ( A[lda+1] * A[2*lda+2] - A[lda+2] * A[2*lda+1]) / det; // A = ei - fh
B[1] = (- A[1] * A[2*lda+2] + A[2] * A[2*lda+1]) / det; // D = -bi + ch
B[2] = ( A[1] * A[lda+2] - A[2] * A[lda+1]) / det; // G = bf - ce
B[ldb] = (- A[lda] * A[2*lda+2] + A[lda+2] * A[2*lda]) / det; // B = -di + fg
B[ldb+1] = ( A[0] * A[2*lda+2] - A[2] * A[2*lda]) / det; // E = ai - cg
B[ldb+2] = (- A[0] * A[lda+2] + A[2] * A[lda]) / det; // H = -af + cd
B[2*ldb] = ( A[lda] * A[2*lda+1] - A[lda+1] * A[2*lda]) / det; // C = dh - eg
B[2*ldb+1]= ( -A[0] * A[2*lda+1] + A[1] * A[2*lda]) / det; // F = -ah + bg
B[2*ldb+2]= ( A[0] * A[lda+1] - A[1] * A[lda]) / det; // I = ae - bd

// Calculate the size of ijb and b, then reallocate them.
IType ndnz = 0;
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
if (j != i && B[i*ldb + j] != 0) { ++ndnz; }
}
}
inverse->ndnz = ndnz;
col_pos = 4; // shape[0] + 1
inverse->capacity = 4 + ndnz;
NM_REALLOC_N(inverse->a, DType, 4 + ndnz);
NM_REALLOC_N(inverse->ija, IType, 4 + ndnz);
b = reinterpret_cast<DType*>(inverse->a);
ijb = reinterpret_cast<IType *>(inverse->ija);

for (int i = 0; i < 3; ++i) {
ijb[i] = col_pos;
for (int j = 0; j < 3; ++j) {
if (j == i) {
b[i] = B[i*ldb + j];
}
else if (B[i*ldb + j] != 0) {
b[col_pos] = B[i*ldb + j];
ijb[col_pos] = j;
++col_pos;
}
}
}
b[3] = 0;
ijb[3] = col_pos;
} else if (M == 1) {
b[0] = 1 / a[0];
} else {
rb_raise(rb_eNotImpError, "exact inverse calculation needed for matrices larger than 3x3");
}
}

/*
* Function signature conversion for calling CBLAS' gesvd functions as directly as possible.
*/
Expand Down Expand Up @@ -1849,14 +2001,27 @@ static VALUE nm_clapack_laswp(VALUE self, VALUE n, VALUE a, VALUE lda, VALUE k1,


/*
* C accessor for calculating an exact determinant.
* C accessor for calculating an exact determinant. Dense matrix version.
*/
void nm_math_det_exact(const int M, const void* elements, const int lda, nm::dtype_t dtype, void* result) {
NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact, void, const int M, const void* A_elements, const int lda, void* result_arg);
void nm_math_det_exact_from_dense(const int M, const void* elements, const int lda,
nm::dtype_t dtype, void* result) {
NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_dense, void, const int M,
const void* A_elements, const int lda, void* result_arg);

ttable[dtype](M, elements, lda, result);
}

/*
* C accessor for calculating an exact determinant. Yale matrix version.
*/
void nm_math_det_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda,
nm::dtype_t dtype, void* result) {
NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_yale, void, const int M,
const YALE_STORAGE* storage, const int lda, void* result_arg);

ttable[dtype](M, storage, lda, result);
}

/*
* C accessor for solving a system of linear equations.
*/
Expand Down Expand Up @@ -1898,14 +2063,29 @@ void nm_math_inverse(const int M, void* a_elements, nm::dtype_t dtype) {
}

/*
* C accessor for calculating an exact inverse.
* C accessor for calculating an exact inverse. Dense matrix version.
*/
void nm_math_inverse_exact(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb, nm::dtype_t dtype) {
NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact, void, const int, const void*, const int, void*, const int);
void nm_math_inverse_exact_from_dense(const int M, const void* A_elements,
const int lda, void* B_elements, const int ldb, nm::dtype_t dtype) {

NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_dense, void,
const int, const void*, const int, void*, const int);

ttable[dtype](M, A_elements, lda, B_elements, ldb);
}

/*
* C accessor for calculating an exact inverse. Yale matrix version.
*/
void nm_math_inverse_exact_from_yale(const int M, const YALE_STORAGE* storage,
const int lda, YALE_STORAGE* inverse, const int ldb, nm::dtype_t dtype) {

NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_yale, void,
const int, const YALE_STORAGE*, const int, YALE_STORAGE*, const int);

ttable[dtype](M, storage, lda, inverse, ldb);
}

/*
* Transpose an array of elements that represent a row-major dense matrix. Does not allocate anything, only does an memcpy.
*/
Expand Down
10 changes: 8 additions & 2 deletions ext/nmatrix/math/math.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,8 +113,14 @@ extern "C" {
void nm_math_solve(VALUE lu, VALUE b, VALUE x, VALUE ipiv);
void nm_math_inverse(const int M, void* A_elements, nm::dtype_t dtype);
void nm_math_hessenberg(VALUE a);
void nm_math_det_exact(const int M, const void* elements, const int lda, nm::dtype_t dtype, void* result);
void nm_math_inverse_exact(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb, nm::dtype_t dtype);
void nm_math_det_exact_from_dense(const int M, const void* elements,
const int lda, nm::dtype_t dtype, void* result);
void nm_math_det_exact_from_yale(const int M, const YALE_STORAGE* storage,
const int lda, nm::dtype_t dtype, void* result);
void nm_math_inverse_exact_from_dense(const int M, const void* A_elements,
const int lda, void* B_elements, const int ldb, nm::dtype_t dtype);
void nm_math_inverse_exact_from_yale(const int M, const YALE_STORAGE* storage,
const int lda, YALE_STORAGE* inverse, const int ldb, nm::dtype_t dtype);
}


Expand Down
Loading