luisgulo/gambas-source-code
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
gambas-source-code/gb.gsl/src/c_matrix.c
Laurent Carlier f1deb8e330 Remove register keyword, deprecated, now reserved keyword with C++17 
Fixes "warning: 'register' storage class specifier is deprecated and incompatible with C++17 [-Wdeprecated-register]"
with clang

https://en.cppreference.com/w/cpp/keyword/register
2019-10-30 06:53:09 +00:00

1300 lines
27 KiB
C
Raw Blame History

/***************************************************************************
c_vector.c
gb.gsl component
(c) 2000-2017 Benoît Minisini <g4mba5@gmail.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 1, or (at your option)
any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
MA 02110-1301, USA.
***************************************************************************/
#define __C_MATRIX_C
#include "c_complex.h"
#include "c_vector.h"
#include "c_matrix.h"
#define THIS ((CMATRIX *)_object)
#define MAT(_m) ((gsl_matrix *)(_m)->matrix)
#define CMAT(_m) ((gsl_matrix_complex *)(_m)->matrix)
#define HEIGHT(_m) ((int)(MAT(_m)->size1))
#define WIDTH(_m) ((int)(MAT(_m)->size2))
#define COMPLEX(_m) ((_m)->complex)
//---- Matrix creation ------------------------------------------------------
static CMATRIX *MATRIX_create(int width, int height, bool complex, bool init)
{
CMATRIX *m = GB.Create(CLASS_Matrix, NULL,NULL);
if (complex)
m->matrix = init ? gsl_matrix_complex_calloc(height, width) : gsl_matrix_complex_alloc(height, width);
else
m->matrix = init ? gsl_matrix_calloc(height, width) : gsl_matrix_alloc(height, width);
m->complex = complex;
return m;
}
static CMATRIX *MATRIX_identity(int width, int height, bool complex)
{
CMATRIX *m = MATRIX_create(width, height, complex, FALSE);
if (complex)
gsl_matrix_complex_set_identity(CMAT(m));
else
gsl_matrix_set_identity(MAT(m));
return m;
}
static CMATRIX *MATRIX_create_from(void *matrix, bool complex)
{
CMATRIX *m = GB.Create(CLASS_Matrix, NULL,NULL);
m->matrix = matrix;
m->complex = complex;
return m;
}
static CMATRIX *MATRIX_copy(CMATRIX *_object)
{
CMATRIX *copy = MATRIX_create(WIDTH(THIS), HEIGHT(THIS), COMPLEX(THIS), FALSE);
if (COMPLEX(THIS))
gsl_matrix_complex_memcpy(CMAT(copy), CMAT(THIS));
else
gsl_matrix_memcpy(MAT(copy), MAT(THIS));
return copy;
}
#define MATRIX_make(_ma) (((_ma)->ob.ref <= 1) ? (_ma) : MATRIX_copy(_ma))
/*static CMATRIX *MATRIX_convert_to_complex(CMATRIX *_object)
{
CMATRIX *m = MATRIX_create(WIDTH(THIS), HEIGHT(THIS), TRUE, FALSE);
int i, j;
for (i = 0; i < HEIGHT(THIS); i++)
for (j = 0; j < WIDTH(THIS); j++)
gsl_matrix_complex_set(CMAT(m), i, j, gsl_complex_rect(gsl_matrix_get(MAT(THIS), i, j), 0));
return m;
}*/
static void MATRIX_ensure_complex(CMATRIX *_object)
{
gsl_matrix_complex *v;
int w = WIDTH(THIS);
int h = HEIGHT(THIS);
int i, j;
if (COMPLEX(THIS))
return;
v = gsl_matrix_complex_alloc(h, w);
for (i = 0; i < h; i++)
for (j = 0; j < w; j++)
gsl_matrix_complex_set(v, i, j, gsl_complex_rect(gsl_matrix_get(MAT(THIS), i, j), 0));
gsl_matrix_free(MAT(THIS));
THIS->matrix = v;
THIS->complex = TRUE;
}
/*static bool MATRIX_ensure_not_complex(CMATRIX *_object)
{
gsl_matrix *m;
int w = WIDTH(THIS);
int h = HEIGHT(THIS);
int i, j;
gsl_complex c;
if (!COMPLEX(THIS))
return FALSE;
for (i = 0; i < h; i++)
for (j = 0; j < w; j++)
{
c = gsl_matrix_complex_get(CMAT(THIS), i, j);
if (GSL_IMAG(c) != 0.0)
return TRUE;
}
m = gsl_matrix_alloc(h, w);
for (i = 0; i < h; i++)
for (j = 0; j < w; j++)
gsl_matrix_set(m, i, j, GSL_REAL(gsl_matrix_complex_get(CMAT(THIS), i, j)));
gsl_matrix_complex_free(CMAT(THIS));
THIS->matrix = m;
THIS->complex = FALSE;
return FALSE;
}*/
static void matrix_negative(void *m, bool complex)
{
uint i;
gsl_matrix *mm = (gsl_matrix *)m;
double *d = mm->data;
uint n = (uint)(mm->size1 * mm->size2);
if (complex)
n *= 2;
for (i = 0; i < n; i++)
d[i] = -d[i];
}
static bool matrix_determinant(CMATRIX *m, COMPLEX_VALUE *det)
{
int sign = 0;
int size = WIDTH(m);
if (size != HEIGHT(m))
return TRUE;
gsl_permutation *p = gsl_permutation_calloc(size);
if (COMPLEX(m))
{
gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(size, size);
gsl_matrix_complex_memcpy(tmp, CMAT(m));
gsl_linalg_complex_LU_decomp(tmp, p, &sign);
det->z = gsl_linalg_complex_LU_det(tmp, sign);
gsl_matrix_complex_free(tmp);
}
else
{
gsl_matrix *tmp = gsl_matrix_alloc(size, size);
gsl_matrix_memcpy(tmp, MAT(m));
gsl_linalg_LU_decomp(tmp, p, &sign);
det->x = gsl_linalg_LU_det(tmp, sign);
det->z.dat[1] = 0;
gsl_matrix_free(tmp);
}
gsl_permutation_free(p);
return FALSE;
}
static void *matrix_invert(void *m, bool complex)
{
int sign = 0;
int size = ((gsl_matrix *)m)->size1;
void *result;
if (size != ((gsl_matrix *)m)->size2)
return NULL;
gsl_permutation *p = gsl_permutation_calloc(size);
if (!complex)
{
gsl_matrix *tmp = gsl_matrix_alloc(size, size);
result = gsl_matrix_alloc(size, size);
gsl_matrix_memcpy(tmp, (gsl_matrix *)m);
gsl_linalg_LU_decomp(tmp, p, &sign);
if (gsl_linalg_LU_invert(tmp, p, (gsl_matrix *)result) != GSL_SUCCESS)
{
gsl_matrix_free(result);
return NULL;
}
gsl_matrix_free(tmp);
}
else
{
gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(size, size);
result = gsl_matrix_complex_alloc(size, size);
gsl_matrix_complex_memcpy(tmp, (gsl_matrix_complex *)m);
gsl_linalg_complex_LU_decomp(tmp, p, &sign);
if (gsl_linalg_complex_LU_invert(tmp, p, (gsl_matrix_complex *)result) != GSL_SUCCESS)
{
gsl_matrix_complex_free(result);
return NULL;
}
gsl_matrix_complex_free(tmp);
}
gsl_permutation_free(p);
return result;
}
static void matrix_add_identity(gsl_matrix *m, double f)
{
gsl_matrix *id = gsl_matrix_alloc(m->size1, m->size2);
gsl_matrix_set_identity(id);
gsl_matrix_scale(id, f);
gsl_matrix_add(m, id);
gsl_matrix_free(id);
}
static void matrix_complex_add_identity(gsl_matrix_complex *m, gsl_complex c)
{
gsl_matrix_complex *id = gsl_matrix_complex_alloc(m->size1, m->size2);
gsl_matrix_complex_set_identity(id);
gsl_matrix_complex_scale(id, c);
gsl_matrix_complex_add(m, id);
gsl_matrix_complex_free(id);
}
//---- Arithmetic operators -------------------------------------------------
#define IMPLEMENT_OP(_name) \
static CMATRIX *_name(CMATRIX *a, CMATRIX *b, bool invert) \
{ \
CMATRIX *m; \
\
if (COMPLEX(a) || COMPLEX(b)) \
{ \
MATRIX_ensure_complex(a); \
MATRIX_ensure_complex(b); \
m = MAKE_MATRIX(a); \
CFUNC(CMAT(m), CMAT(a), CMAT(b)); \
} \
else \
{ \
m = MAKE_MATRIX(a); \
FUNC(MAT(m), MAT(a), MAT(b)); \
} \
\
return m; \
}
#define IMPLEMENT_OP_FLOAT(_name) \
static CMATRIX *_name(CMATRIX *a, double f, bool invert) \
{ \
CMATRIX *m = MAKE_MATRIX(a); \
\
if (COMPLEX(a)) \
{ \
CFUNC(CMAT(m), CMAT(a), f); \
} \
else \
{ \
FUNC(MAT(m), MAT(a), f); \
} \
\
return m; \
}
#define IMPLEMENT_OP_OTHER(_name) \
static CMATRIX *_name(CMATRIX *a, void *b, bool invert) \
{ \
CMATRIX *m = MAKE_MATRIX(a); \
\
if (GB.Is(b, CLASS_Complex)) \
{ \
MATRIX_ensure_complex(m); \
CFUNC(CMAT(m), CMAT(a), ((CCOMPLEX *)b)->number); \
return m; \
} \
else \
return NULL; \
}
#define MAKE_MATRIX(_a) MATRIX_make(_a)
#define FUNC(_m, _a, _b) gsl_matrix_add(_m, _b)
#define CFUNC(_m, _a, _b) gsl_matrix_complex_add(_m, _b)
IMPLEMENT_OP(_add)
#undef FUNC
#undef CFUNC
#define FUNC(_m, _a, _f) matrix_add_identity(_m, _f)
#define CFUNC(_m, _a, _f) matrix_complex_add_identity(_m, gsl_complex_rect(_f, 0))
IMPLEMENT_OP_FLOAT(_addf)
#undef FUNC
#undef CFUNC
#define CFUNC(_m, _a, _c) matrix_complex_add_identity(_m, _c)
IMPLEMENT_OP_OTHER(_addo)
#undef CFUNC
#define FUNC(_m, _a, _b) gsl_matrix_sub(_m, _b)
#define CFUNC(_m, _a, _b) gsl_matrix_complex_sub(_m, _b)
IMPLEMENT_OP(_sub)
#undef FUNC
#undef CFUNC
#define FUNC(_m, _a, _f) \
if (invert) \
{ \
matrix_negative(_m, FALSE); \
matrix_add_identity(_m, _f); \
} \
else \
matrix_add_identity(_m, -(_f));
#define CFUNC(_m, _a, _f) \
if (invert) \
{ \
matrix_negative(_m, TRUE); \
matrix_complex_add_identity(_m, gsl_complex_rect(_f, 0)); \
} \
else \
matrix_complex_add_identity(_m, gsl_complex_rect(-(_f), 0));
IMPLEMENT_OP_FLOAT(_subf)
#undef FUNC
#undef CFUNC
#define CFUNC(_m, _a, _c) \
if (invert) \
matrix_negative(_m, TRUE); \
else \
gsl_complex_negative(_c); \
matrix_complex_add_identity(_m, _c);
IMPLEMENT_OP_OTHER(_subo)
#undef CFUNC
#define FUNC(_m, _a, _f) gsl_matrix_scale(_m, _f)
#define CFUNC(_m, _a, _f) gsl_matrix_complex_scale(_m, gsl_complex_rect(_f, 0))
IMPLEMENT_OP_FLOAT(_mulf)
#undef FUNC
#undef CFUNC
#define CFUNC(_m, _a, _c) gsl_matrix_complex_scale(_m, _c)
IMPLEMENT_OP_OTHER(_mulo)
#undef CFUNC
#undef MAKE_MATRIX
#define MAKE_MATRIX(_a) MATRIX_copy(_a)
#define FUNC(_m, _a, _b) gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, _a, _b, 0.0, _m);
#define CFUNC(_m, _a, _b) gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, COMPLEX_one, _a, _b, COMPLEX_zero, _m);
IMPLEMENT_OP(_mul)
#undef FUNC
#undef CFUNC
#define FUNC(_m, _a, _b) \
{ \
gsl_matrix *inv = matrix_invert(_b, FALSE); \
if (!inv) \
return NULL; \
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, _a, inv, 0.0, _m); \
gsl_matrix_free(inv); \
}
#define CFUNC(_m, _a, _b) \
{ \
gsl_matrix_complex *inv = matrix_invert(_b, TRUE); \
if (!inv) \
{ \
GB.Error(GB_ERR_ZERO); \
return NULL; \
} \
gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, COMPLEX_one, _a, inv, COMPLEX_zero, _m); \
gsl_matrix_complex_free(inv); \
}
IMPLEMENT_OP(_div)
#undef FUNC
#undef CFUNC
static CMATRIX *_divf(CMATRIX *a, double f, bool invert)
{
bool complex = COMPLEX(a);
CMATRIX *m;
if (invert)
{
void *inv = matrix_invert(MAT(a), complex);
if (!inv)
{
GB.Error(GB_ERR_ZERO);
return NULL;
}
m = MATRIX_create_from(inv, complex);
}
else
{
if (f == 0.0)
{
GB.Error(GB_ERR_ZERO);
return NULL;
}
f = 1 / f;
m = MATRIX_make(a);
}
if (complex)
gsl_matrix_complex_scale(CMAT(m), gsl_complex_rect(f, 0));
else
gsl_matrix_scale(MAT(m), f);
return m;
}
static CMATRIX *_divo(CMATRIX *a, void *b, bool invert)
{
bool complex = COMPLEX(a);
CMATRIX *m;
gsl_complex c;
if (!GB.Is(b, CLASS_Complex))
return NULL;
c = ((CCOMPLEX *)b)->number;
if (invert)
{
void *inv = matrix_invert(MAT(a), complex);
if (!inv)
{
GB.Error(GB_ERR_ZERO);
return NULL;
}
m = MATRIX_create_from(inv, complex);
}
else
{
if (GSL_REAL(c) == 0 && GSL_IMAG(c) == 0)
{
GB.Error(GB_ERR_ZERO);
return NULL;
}
c = gsl_complex_inverse(c);
m = MATRIX_make(a);
}
MATRIX_ensure_complex(m);
gsl_matrix_complex_scale(CMAT(m), c);
return m;
}
static int _equal(CMATRIX *a, CMATRIX *b)
{
if (WIDTH(a) != WIDTH(b) || HEIGHT(a) != HEIGHT(b))
return FALSE;
if (COMPLEX(a) || COMPLEX(b))
{
MATRIX_ensure_complex(a);
MATRIX_ensure_complex(b);
return gsl_matrix_complex_equal(CMAT(a), CMAT(b));
}
else
return gsl_matrix_equal(MAT(a), MAT(b));
}
static int _equalf(CMATRIX *a, double f)
{
bool result;
if (COMPLEX(a))
{
if (f == 0.0)
return gsl_matrix_complex_isnull(CMAT(a));
gsl_matrix_complex *m = gsl_matrix_complex_alloc(WIDTH(a), HEIGHT(a));
gsl_matrix_complex_set_identity(m);
gsl_matrix_complex_scale(m, gsl_complex_rect(f, 0));
result = gsl_matrix_complex_equal(CMAT(a), m);
gsl_matrix_complex_free(m);
}
else
{
if (f == 0.0)
return gsl_matrix_isnull(MAT(a));
gsl_matrix *m = gsl_matrix_alloc(WIDTH(a), HEIGHT(a));
gsl_matrix_set_identity(m);
gsl_matrix_scale(m, f);
result = gsl_matrix_equal(MAT(a), m);
gsl_matrix_free(m);
}
return result;
}
static int _equalo(CMATRIX *a, void *b)
{
bool result;
CCOMPLEX *c;
if (!GB.Is(b, CLASS_Complex))
return -1;
c = (CCOMPLEX *)b;
if (GSL_IMAG(c->number) == 0.0)
return _equalf(a, GSL_REAL(c->number));
if (!COMPLEX(a))
return FALSE;
gsl_matrix_complex *m = gsl_matrix_complex_alloc(WIDTH(a), HEIGHT(a));
gsl_matrix_complex_set_identity(m);
gsl_matrix_complex_scale(m, c->number);
result = gsl_matrix_complex_equal(CMAT(a), m);
gsl_matrix_complex_free(m);
return result;
}
static CMATRIX *_neg(CMATRIX *a)
{
CMATRIX *m = MATRIX_make(a);
matrix_negative(m->matrix, m->complex);
return m;
}
static CMATRIX *_powi(CMATRIX *m, int n)
{
if (n == 1)
return m;
CMATRIX *m2 = _mul(m, m, FALSE);
CMATRIX *r;
if ((n & 1) == 0)
{
n /= 2;
if (n > 1)
r = _powi(m2, n);
else
r = m2;
}
else
{
n /= 2;
if (n > 1)
r = _powi(m2, n);
else
r = m2;
m2 = _mul(r, m, FALSE);
GB.Unref(POINTER(&r));
r = m2;
}
GB.Unref(POINTER(&m));
return r;
}
static CMATRIX *_powf(CMATRIX *a, double f, bool invert)
{
if (invert || f != (double)(int)f)
return NULL;
CMATRIX *m = NULL;
int n = (int)f;
if (n == 0)
{
m = MATRIX_make(a);
if (COMPLEX(m))
gsl_matrix_complex_set_identity(CMAT(m));
else
gsl_matrix_set_identity(MAT(m));
}
else if (n == 1)
{
m = a;
}
else if (n > 1)
{
m = _powi(MATRIX_copy(a), n);
}
else if (n < 0)
{
void *inv = matrix_invert(a->matrix, COMPLEX(a));
if (inv == NULL)
{
GB.Error(GB_ERR_ZERO);
return NULL;
}
m = _powi(MATRIX_create_from(inv, COMPLEX(a)), (-n));
}
return m;
}
static GB_OPERATOR_DESC _operator =
{
.equal = (void *)_equal,
.equalf = (void *)_equalf,
.equalo = (void *)_equalo,
.add = (void *)_add,
.addf = (void *)_addf,
.addo = (void *)_addo,
.sub = (void *)_sub,
.subf = (void *)_subf,
.subo = (void *)_subo,
.mul = (void *)_mul,
.mulf = (void *)_mulf,
.mulo = (void *)_mulo,
.div = (void *)_div,
.divf = (void *)_divf,
.divo = (void *)_divo,
.powf = (void *)_powf,
.neg = (void *)_neg
};
//---- Conversions ----------------------------------------------------------
static char *_to_string(CMATRIX *_object, bool local)
{
char *result = NULL;
int i, j;
int w = WIDTH(THIS);
int h = HEIGHT(THIS);
char *str;
int len;
result = GB.AddChar(result, '[');
for (i = 0; i < h; i++)
{
if (i)
{
if (!local)
result = GB.AddChar(result, ',');
}
result = GB.AddChar(result, '[');
for (j = 0; j < w; j++)
{
if (j)
result = GB.AddChar(result, local ? ' ' : ',');
if (!COMPLEX(THIS))
{
GB.NumberToString(local, gsl_matrix_get(MAT(THIS), i, j), NULL, &str, &len);
result = GB.AddString(result, str, len);
}
else
{
str = COMPLEX_to_string(gsl_matrix_complex_get(CMAT(THIS), i, j), local);
result = GB.AddString(result, str, GB.StringLength(str));
GB.FreeString(&str);
}
}
result = GB.AddChar(result, ']');
}
result = GB.AddChar(result, ']');
return result;
}
static bool _convert(CMATRIX *_object, GB_TYPE type, GB_VALUE *conv)
{
if (THIS)
{
if (!COMPLEX(THIS))
{
switch (type)
{
/*case GB_T_FLOAT:
conv->_float.value = gsl_blas_dnrm2(MAT(THIS));
return FALSE;
case GB_T_SINGLE:
conv->_single.value = gsl_blas_dnrm2(MAT(THIS));
return FALSE;
case GB_T_INTEGER:
case GB_T_SHORT:
case GB_T_BYTE:
conv->_integer.value = gsl_blas_dnrm2(MAT(THIS));
return FALSE;
case GB_T_LONG:
conv->_long.value = gsl_blas_dnrm2(MAT(THIS));
return FALSE;*/
case GB_T_STRING:
case GB_T_CSTRING:
conv->_string.value.addr = _to_string(THIS, type == GB_T_CSTRING);
conv->_string.value.start = 0;
conv->_string.value.len = GB.StringLength(conv->_string.value.addr);
return FALSE;
default:
return TRUE;
}
}
else
{
switch (type)
{
/*case GB_T_FLOAT:
conv->_float.value = gsl_blas_dznrm2(CMAT(THIS));
return FALSE;
case GB_T_SINGLE:
conv->_single.value = gsl_blas_dznrm2(CMAT(THIS));
return FALSE;
case GB_T_INTEGER:
case GB_T_SHORT:
case GB_T_BYTE:
conv->_integer.value = gsl_blas_dznrm2(CMAT(THIS));
return FALSE;
case GB_T_LONG:
conv->_long.value = gsl_blas_dznrm2(CMAT(THIS));
return FALSE;*/
case GB_T_STRING:
case GB_T_CSTRING:
conv->_string.value.addr = _to_string(THIS, type == GB_T_CSTRING);
conv->_string.value.start = 0;
conv->_string.value.len = GB.StringLength(conv->_string.value.addr);
return FALSE;
default:
return TRUE;
}
}
}
else if (type >= GB_T_OBJECT)
{
/*if (type == CLASS_Complex)
{
CCOMPLEX *c = (CCOMPLEX *)conv->_object.value;
CMATRIX *m = MATRIX_create(2, 2, FALSE, FALSE);
gsl_matrix_set(MAT(m), 0, 0, GSL_REAL(c->number));
gsl_matrix_set(MAT(m), 1, 1, GSL_REAL(c->number));
gsl_matrix_set(MAT(m), 0, 1, -GSL_IMAG(c->number));
gsl_matrix_set(MAT(m), 1, 0, GSL_IMAG(c->number));
conv->_object.value = m;
return FALSE;
}
else*/ if (GB.Is(conv->_object.value, CLASS_Array))
{
GB_ARRAY array = (GB_ARRAY)conv->_object.value;
GB_ARRAY array2;
int height = GB.Array.Count(array);
int width = 0;
GB_TYPE atype = GB.Array.Type(array);
GB_TYPE atype2;
int i, j, w;
GB_VALUE temp;
void *data;
CMATRIX *m;
CCOMPLEX *c;
bool complex;
if (atype >= GB_T_OBJECT)
{
complex = FALSE;
for (i = 0; i < height; i++)
{
data = GB.Array.Get(array, i);
array2 = *((void **)data);
if (!array2 || !GB.Is(array2, CLASS_Array))
return TRUE;
w = GB.Array.Count(array2);
if (w > width)
width = w;
atype2 = GB.Array.Type(array2);
if (atype2 == GB_T_VARIANT || atype2 == CLASS_Complex)
complex = TRUE;
else if (!(atype2 > GB_T_BOOLEAN && atype2 <= GB_T_FLOAT))
return TRUE;
}
//fprintf(stderr, "create: %d %d %d\n", width, height, complex);
m = MATRIX_create(width, height, complex, TRUE);
for (i = 0; i < height; i++)
{
array2 = *((void **)GB.Array.Get(array, i));
atype2 = GB.Array.Type(array2);
w = GB.Array.Count(array2);
if (atype2 > GB_T_BOOLEAN && atype2 <= GB_T_FLOAT)
{
for (j = 0; j < w; j++)
{
data = GB.Array.Get(array2, j);
GB.ReadValue(&temp, data, atype2);
GB.Conv(&temp, GB_T_FLOAT);
if (complex)
gsl_matrix_complex_set(CMAT(m), i, j, gsl_complex_rect(temp._float.value, 0));
else
gsl_matrix_set(MAT(m), i, j, temp._float.value);
}
}
else if (atype2 == GB_T_VARIANT)
{
for (j = 0; j < w; j++)
{
GB.ReadValue(&temp, GB.Array.Get(array2, j), atype2);
GB.BorrowValue(&temp);
GB.Conv(&temp, CLASS_Complex);
c = temp._object.value;
if (c)
gsl_matrix_complex_set(CMAT(m), i, j, c->number);
else
gsl_matrix_complex_set(CMAT(m), i, j, COMPLEX_zero);
GB.ReleaseValue(&temp);
}
}
else if (atype2 == CLASS_Complex)
{
for (j = 0; j < w; j++)
{
c = *((CCOMPLEX **)GB.Array.Get(array2, j));
if (c)
gsl_matrix_complex_set(CMAT(m), i, j, c->number);
else
gsl_matrix_complex_set(CMAT(m), i, j, COMPLEX_zero);
}
}
}
conv->_object.value = m;
return FALSE;
}
}
}
/*else if (type > GB_T_BOOLEAN && type <= GB_T_FLOAT)
{
CMATRIX *m = MATRIX_create(2, 2, FALSE, TRUE);
double value;
if (type == GB_T_FLOAT)
value = conv->_float.value;
else if (type == GB_T_SINGLE)
value = conv->_single.value;
else
value = conv->_integer.value;
gsl_matrix_set(MAT(m), 0, 0, value);
gsl_matrix_set(MAT(m), 1, 1, value);
conv->_object.value = m;
return FALSE;
}*/
return TRUE;
}
//---------------------------------------------------------------------------
BEGIN_METHOD(Matrix_new, GB_INTEGER height; GB_INTEGER width; GB_BOOLEAN complex)
bool complex = VARGOPT(complex, FALSE);
int h = VARGOPT(height, 2);
int w = VARGOPT(width, 2);
if (h < 1) h = 1;
if (w < 1) w = 1;
THIS->complex = complex;
if (!complex)
THIS->matrix = gsl_matrix_calloc(h, w);
else
THIS->matrix = gsl_matrix_complex_calloc(h, w);
END_METHOD
BEGIN_METHOD_VOID(Matrix_free)
if (!COMPLEX(THIS))
gsl_matrix_free(MAT(THIS));
else
gsl_matrix_complex_free(CMAT(THIS));
END_METHOD
BEGIN_PROPERTY(Matrix_Width)
GB.ReturnInteger(WIDTH(THIS));
END_PROPERTY
BEGIN_PROPERTY(Matrix_Height)
GB.ReturnInteger(HEIGHT(THIS));
END_PROPERTY
BEGIN_METHOD_VOID(Matrix_Copy)
GB.ReturnObject(MATRIX_copy(THIS));
END_METHOD
BEGIN_METHOD(Matrix_get, GB_INTEGER i; GB_INTEGER j)
int w = WIDTH(THIS), h = HEIGHT(THIS);
int i = VARG(i), j = VARG(j);
if (i < 0 || i >= h || j < 0 || j >= w)
{
GB.Error(GB_ERR_BOUND);
return;
}
if (!COMPLEX(THIS))
GB.ReturnFloat(gsl_matrix_get(MAT(THIS), i, j));
else
GB.ReturnObject(COMPLEX_create(gsl_matrix_complex_get(CMAT(THIS), i, j)));
GB.ReturnConvVariant();
END_METHOD
BEGIN_METHOD(Matrix_put, GB_VARIANT value; GB_INTEGER i; GB_INTEGER j)
int w = WIDTH(THIS), h = HEIGHT(THIS);
int i = VARG(i), j = VARG(j);
GB_VALUE *value = (GB_VALUE *)ARG(value);
int type;
COMPLEX_VALUE cv;
if (i < 0 || i >= h || j < 0 || j >= w)
{
GB.Error(GB_ERR_BOUND);
return;
}
type = COMPLEX_get_value(value, &cv);
if (type == CGV_ERR)
return;
if (type == CGV_COMPLEX)
{
MATRIX_ensure_complex(THIS);
gsl_matrix_complex_set(CMAT(THIS), i, j, cv.z);
}
else
{
if (COMPLEX(THIS))
gsl_matrix_complex_set(CMAT(THIS), i, j, cv.z);
else
gsl_matrix_set(MAT(THIS), i, j, cv.x);
}
END_METHOD
BEGIN_PROPERTY(Matrix_Handle)
GB.ReturnPointer(THIS->matrix);
END_PROPERTY
BEGIN_METHOD(Matrix_Identity, GB_INTEGER width; GB_INTEGER height; GB_BOOLEAN complex)
GB.ReturnObject(MATRIX_identity(VARGOPT(width, 2), VARGOPT(height, 2), VARGOPT(complex, FALSE)));
END_METHOD
BEGIN_METHOD(Matrix_Row, GB_INTEGER row)
int row = VARG(row);
if (row < 0 || row >= HEIGHT(THIS))
{
GB.Error(GB_ERR_BOUND);
return;
}
bool complex = COMPLEX(THIS);
CVECTOR *v = VECTOR_create(WIDTH(THIS), complex, FALSE);
if (complex)
gsl_matrix_complex_get_row(CVEC(v), CMAT(THIS), row);
else
gsl_matrix_get_row(VEC(v), MAT(THIS), row);
GB.ReturnObject(v);
END_METHOD
BEGIN_METHOD(Matrix_Column, GB_INTEGER column)
int column = VARG(column);
if (column < 0 || column >= WIDTH(THIS))
{
GB.Error(GB_ERR_BOUND);
return;
}
bool complex = COMPLEX(THIS);
CVECTOR *v = VECTOR_create(HEIGHT(THIS), complex, FALSE);
if (complex)
gsl_matrix_complex_get_col(CVEC(v), CMAT(THIS), column);
else
gsl_matrix_get_col(VEC(v), MAT(THIS), column);
GB.ReturnObject(v);
END_METHOD
BEGIN_METHOD(Matrix_SetRow, GB_INTEGER row; GB_OBJECT vector)
int row = VARG(row);
if (row < 0 || row >= HEIGHT(THIS))
{
GB.Error(GB_ERR_BOUND);
return;
}
CVECTOR *v = (CVECTOR *)VARG(vector);
if (GB.CheckObject(v))
return;
if (SIZE(v) != WIDTH(THIS))
{
GB.Error("Vector size does not match matrix width");
return;
}
bool complex = COMPLEX(THIS);
if (complex)
{
VECTOR_ensure_complex(v);
gsl_matrix_complex_set_row(CMAT(THIS), row, CVEC(v));
}
else
{
if (VECTOR_ensure_not_complex(v))
{
GB.Error(GB_ERR_TYPE, "Float", "Complex");
return;
}
gsl_matrix_set_row(MAT(THIS), row, VEC(v));
}
END_METHOD
BEGIN_METHOD(Matrix_SetColumn, GB_INTEGER column; GB_OBJECT vector)
int column = VARG(column);
if (column < 0 || column >= WIDTH(THIS))
{
GB.Error(GB_ERR_BOUND);
return;
}
CVECTOR *v = (CVECTOR *)VARG(vector);
if (GB.CheckObject(v))
return;
if (SIZE(v) != HEIGHT(THIS))
{
GB.Error("Vector size does not match matrix height");
return;
}
bool complex = COMPLEX(THIS);
if (complex)
{
VECTOR_ensure_complex(v);
gsl_matrix_complex_set_col(CMAT(THIS), column, CVEC(v));
}
else
{
if (VECTOR_ensure_not_complex(v))
{
GB.Error(GB_ERR_TYPE, "Float", "Complex");
return;
}
gsl_matrix_set_col(MAT(THIS), column, VEC(v));
}
END_METHOD
BEGIN_METHOD_VOID(Matrix_Determinant)
COMPLEX_VALUE cv;
if (matrix_determinant(THIS, &cv))
{
GB.Error("Matrix is not square");
return;
}
if (COMPLEX(THIS))
GB.ReturnObject(COMPLEX_create(cv.z));
else
GB.ReturnFloat(cv.x);
GB.ReturnConvVariant();
END_METHOD
BEGIN_METHOD(Matrix_call, GB_OBJECT vector)
CVECTOR *v = VARG(vector);
CVECTOR *result;
if (GB.CheckObject(v))
return;
if (COMPLEX(THIS) || v->complex)
{
MATRIX_ensure_complex(THIS);
VECTOR_ensure_complex(v);
result = VECTOR_create(SIZE(v), TRUE, FALSE);
gsl_blas_zgemv(CblasNoTrans, COMPLEX_one, CMAT(THIS), CVEC(v), COMPLEX_zero, CVEC(result));
GB.ReturnObject(result);
}
else
{
result = VECTOR_create(SIZE(v), FALSE, FALSE);
gsl_blas_dgemv(CblasNoTrans, 1.0, MAT(THIS), VEC(v), 0.0, VEC(result));
GB.ReturnObject(result);
}
END_METHOD
BEGIN_METHOD_VOID(Matrix_Transpose)
if (!COMPLEX(THIS))
{
gsl_matrix *m = gsl_matrix_alloc(WIDTH(THIS), HEIGHT(THIS));
gsl_matrix_transpose_memcpy(m, MAT(THIS));
GB.ReturnObject(MATRIX_create_from(m, FALSE));
}
else
{
gsl_matrix_complex *m = gsl_matrix_complex_alloc(WIDTH(THIS), HEIGHT(THIS));
gsl_matrix_complex_transpose_memcpy(m, CMAT(THIS));
GB.ReturnObject(MATRIX_create_from(m, TRUE));
}
END_METHOD
BEGIN_METHOD_VOID(Matrix_Conjugate)
CMATRIX *m = MATRIX_copy(THIS);
if (COMPLEX(THIS))
{
int i, j;
for (i = 0; i < HEIGHT(m); i++)
for (j = 0; j < WIDTH(m); j++)
gsl_matrix_complex_set(CMAT(m), i, j, gsl_complex_conjugate(gsl_matrix_complex_get(CMAT(m), i, j)));
}
GB.ReturnObject(m);
END_METHOD
BEGIN_METHOD_VOID(Matrix_Invert)
void *m = matrix_invert(THIS->matrix, COMPLEX(THIS));
if (!m)
GB.ReturnNull();
else
GB.ReturnObject(MATRIX_create_from(m, COMPLEX(THIS)));
END_METHOD
BEGIN_METHOD(Matrix_ToString, GB_BOOLEAN local)
GB.ReturnString(GB.FreeStringLater(_to_string(THIS, VARGOPT(local, FALSE))));
END_METHOD
//---------------------------------------------------------------------------
GB_DESC MatrixDesc[] =
{
GB_DECLARE("Matrix", sizeof(CMATRIX)),
GB_METHOD("_new", NULL, Matrix_new, "[(Width)i(Height)i(Complex)b]"),
GB_METHOD("_free", NULL, Matrix_free, NULL),
//GB_STATIC_METHOD("_call", "Vector", Matrix_call, "(Value)f."),
GB_METHOD("Copy", "Matrix", Matrix_Copy, NULL),
GB_METHOD("ToString", "s", Matrix_ToString, "[(Local)b]"),
GB_STATIC_METHOD("Identity", "Matrix", Matrix_Identity, "[(Width)i(Height)i(Complex)b]"),
GB_PROPERTY_READ("Width", "i", Matrix_Width),
GB_PROPERTY_READ("Height", "i", Matrix_Height),
GB_PROPERTY_READ("Handle", "p", Matrix_Handle),
GB_METHOD("Det", "v", Matrix_Determinant, NULL),
GB_METHOD("_get", "v", Matrix_get, "(I)i(J)i"),
GB_METHOD("_put", NULL, Matrix_put, "(Value)v(I)i(J)i"),
GB_METHOD("_call", "Vector", Matrix_call, "(Vector)Vector"),
//GB_METHOD("Equal", "b", Matrix_Equal, "(Matrix)Matrix;"),
GB_METHOD("Row", "Vector", Matrix_Row, "(Row)i"),
GB_METHOD("Column", "Vector", Matrix_Column, "(Column)i"),
GB_METHOD("SetRow", NULL, Matrix_SetRow, "(Row)i(Vector)Vector;"),
GB_METHOD("SetColumn", NULL, Matrix_SetColumn, "(Column)i(Vector)Vector;"),
//GB_METHOD("Scale", "Matrix", Matrix_Scale, "(Value)v"),
GB_METHOD("Trans", "Matrix", Matrix_Transpose, NULL),
GB_METHOD("Conj", "Matrix", Matrix_Conjugate, NULL),
GB_METHOD("Inv", "Matrix", Matrix_Invert, NULL),
GB_INTERFACE("_convert", &_convert),
GB_INTERFACE("_operator", &_operator),
GB_END_DECLARE
};
Reference in a new issue View git blame Copy permalink