Commit b20e9ec2 authored by Jiri (George) Lebl's avatar Jiri (George) Lebl Committed by George Lebl

fix possible memory corruption in the gauss routine


Wed Nov 14 20:13:29 2007  Jiri (George) Lebl <jirka@5z.com>

	* src/matop.c: fix possible memory corruption in the gauss
	  routine

	* src/matrixw.c: slight optimization to not copy zeros

	* src/funclib.c, lib/linear_algebra/*.gel: Implement NullSpace
	  and PivotColumns internally for speed.

	* src/gnome-genius.c: remove some unused vars

	* src/geniustests.txt, src/nullspacetests.txt: add tests for
	  PivotColumns and NullSpace


svn path=/trunk/; revision=612
parent e3ed0741
Wed Nov 14 20:13:29 2007 Jiri (George) Lebl <jirka@5z.com>
* src/matop.c: fix possible memory corruption in the gauss
routine
* src/matrixw.c: slight optimization to not copy zeros
* src/funclib.c, lib/linear_algebra/*.gel: Implement NullSpace
and PivotColumns internally for speed.
* src/gnome-genius.c: remove some unused vars
* src/geniustests.txt, src/nullspacetests.txt: add tests for
PivotColumns and NullSpace
Wed Nov 14 04:15:49 2007 Jiri (George) Lebl <jirka@5z.com>
* lib/linear_algebra/*.gel: OrthogonalComplement is with respect
......
......@@ -10,6 +10,7 @@ Changes to 1.0.1
these are floating point by appending .0 when they would look like integers.
* Add IsMatrixPositive, IsMatrixNonnegative, version, IsZero, IsIdentity,
DividePoly, IsSubset
* PivotColumns and NullSpace are built in for greater speed.
* CHANGE: Remove IsGaussianInteger, we already have IsGaussInteger alias
IsComplexInteger
* CHANGE: OrthogonalComplement is with respect to the Hermitian product
......
......@@ -129,10 +129,8 @@ char *fake = N_("Get the LU decomposition of A and store the result in the L and
char *fake = N_("Get the i-j minor of a matrix");
char *fake = N_("Return the columns that are not the pivot columns of a matrix");
char *fake = N_("Get the p Norm (or 2 Norm if no p is supplied) of a vector");
char *fake = N_("Get the nullspace of a matrix");
char *fake = N_("Get the nullity of a matrix");
char *fake = N_("Get the orthogonal complement of the columnspace");
char *fake = N_("Return pivot columns of a matrix, that is columns which have a leading 1 in rref form, also returns the row where they occur");
char *fake = N_("Projection of vector v onto subspace W given a sesquilinear form B (if not given use hermitian product)");
char *fake = N_("Get the QR decomposition of A, returns R and Q can be a reference");
char *fake = N_("Get the rank of a matrix");
......
......@@ -12,28 +12,7 @@ protect ("Conjugate Transpose")
#FIXME: We'll need more sophisticated tools for Tensors
# Pivot Columns
# Given a matrix in rref form, the columns which have a leading 1
# in some row are the pivot columns.
# (also returns in which row they occur)
SetHelp ("PivotColumns", "linear_algebra", "Return pivot columns of a matrix, that is columns which have a leading 1 in rref form, also returns the row where they occur")
function PivotColumns(M) =
(
N=rref(M);
col_dim_N=columns(N);
pivot_columns=.;
for loop=1 to min(rows(N),col_dim_N) do
(
for inner_loop=loop to col_dim_N do (
if N@(loop,inner_loop) != 0 then (
pivot_columns=[pivot_columns,[inner_loop;loop]];
break
)
)
);
pivot_columns
)
protect ("PivotColumns")
# PivotColumns is now built-in for speed
# Non-pivot columns
SetHelp ("NonPivotColumns", "linear_algebra", "Return the columns that are not the pivot columns of a matrix")
......
......@@ -62,60 +62,10 @@ Image = ColumnSpace
protect ("Image")
# Null space/kernel of a linear transform
# Okay, here's the idea:
# Row reduce a matrix. Then the non-pivot columns are basically
# the independent variables, and the pivot columns are the dependent ones.
# So if your row reduced matrix looks like this:
# [1 0 0 2 4]
# [0 0 1 -3 5]
# then to find a basis for the kernel, look at your non-pivot columns
# (4, 5)
# and for each non-pivot column, you get one vector.
# So take the fourth column, and start off with the vector [0,0,0,-1,0].'
# (so a -1 in the fourth place)
# Now in each pivot entry, you need to put a value to cancel what this
# -1 gives you -- so the pivot column entries are 2 and -3 (the entries
# of the fourth column that have a pivot to the left of them).
# So the first vector is [2,0,-3,-1,0], and the second is
# [4,0,5,0,-1]
# This is poorly explained (FIXME), but some examples should make it
# clear (find a good reference for this!)
SetHelp ("NullSpace", "linear_algebra", "Get the nullspace of a matrix")
function NullSpace (T)=
(
T = rref (T);
pivots = PivotColumns (T);
dim_image = columns (T);
number_of_pivots = columns (pivots);
# injective transforms, so no kernel
if dim_image == number_of_pivots or dim_image == 0 then return null;
non_pivots = SetMinus ([1:columns(T)], pivots@(1,));
null_space = zeros(dim_image,dim_image - number_of_pivots);
j = 1;
for loop in non_pivots do
(
new_vector=zeros(dim_image,1);
new_vector@(loop)=-1;
for inner_loop=1 to number_of_pivots do
if pivots@(1,inner_loop) < loop
then new_vector@(pivots@(1,inner_loop))
=T@(pivots@(2,inner_loop),loop)
else break;
null_space@(,j)=new_vector;
j = j+1
);
#ColumnSpace(null_space)
# Hmmm, no need to use ColumnSpace that does another bunch of operations,
# if the recipient wants to reduce this, they should do it themselves, what
# we give them are lin ind already
null_space
)
protect ("NullSpace")
# NullSpace is now built-in for speed
SetHelp ("Kernel", "linear_algebra", "Get the kernel (nullspace) of a linear transform")
Kernel = NullSpace
function Kernel(M) = NullSpace(M)
protect ("Kernel")
## Vector Subspace operations
......
......@@ -3244,6 +3244,210 @@ rref_op(GelCtx *ctx, GelETree * * a, gboolean *exception)
return n;
}
/* cols and rows should have enough space (min(cols,rows) of m)
* and m should be in at least ref (if not rref) form) and value only,
* returns the count. The values returned are zero based! */
static int
get_pivot_cols (GelMatrixW *m, int *cols, int *rows)
{
int i, j, w, h, mwh;
int cnt = 0;
w = gel_matrixw_width (m);
h = gel_matrixw_height (m);
mwh = MIN (w, h);
for (j = 0; j < mwh; j++) {
for (i = j; i < w; i++) {
GelETree *t = gel_matrixw_get_index (m, i, j);
if (t != NULL &&
! mpw_zero_p (t->val.value)) {
cols[cnt] = i;
rows[cnt] = j;
cnt++;
break;
}
}
}
return cnt;
}
/* PivotColumns
* Given a matrix in rref form, the columns which have a leading 1
* in some row are the pivot columns.
* (also returns in which row they occur) */
static GelETree *
PivotColumns_op(GelCtx *ctx, GelETree * * a, gboolean *exception)
{
GelETree *n;
GelMatrixW *m;
GelMatrix *nm;
gboolean copied_m = FALSE;
int *cols, *rows;
int cnt, mwh;
int i;
if G_UNLIKELY (a[0]->type == NULL_NODE)
return gel_makenum_null ();
if G_UNLIKELY ( ! check_argument_value_only_matrix (a, 0, "PivotColumns"))
return NULL;
m = a[0]->mat.matrix;
if ( ! m->rref) {
m = gel_matrixw_copy (m);
/* only do ref, not rref for speed */
gel_value_matrix_gauss (ctx, m, FALSE, FALSE, FALSE, NULL, NULL);
copied_m = TRUE;
}
mwh = MIN (gel_matrixw_width (m), gel_matrixw_height (m));
cols = g_new (int, mwh);
rows = g_new (int, mwh);
cnt = get_pivot_cols (m, cols, rows);
if (copied_m)
gel_matrixw_free (m);
if (cnt == 0) {
g_free (cols);
g_free (rows);
return gel_makenum_null ();
}
nm = gel_matrix_new ();
gel_matrix_set_size (nm, cnt, 2, FALSE /* padding */);
for (i = 0; i < cnt; i++) {
gel_matrix_index (nm, i, 0) = gel_makenum_ui (cols[i]+1);
gel_matrix_index (nm, i, 1) = gel_makenum_ui (rows[i]+1);
}
GET_NEW_NODE (n);
n->type = MATRIX_NODE;
n->mat.matrix = gel_matrixw_new_with_matrix (nm);
n->mat.quoted = FALSE;
g_free (cols);
g_free (rows);
return n;
}
/*
# Null space/kernel of a linear transform
# Okay, here's the idea:
# Row reduce a matrix. Then the non-pivot columns are basically
# the independent variables, and the pivot columns are the dependent ones.
# So if your row reduced matrix looks like this:
# [1 0 0 2 4]
# [0 0 1 -3 5]
# then to find a basis for the kernel, look at your non-pivot columns
# (4, 5)
# and for each non-pivot column, you get one vector.
# So take the fourth column, and start off with the vector [0,0,0,-1,0].'
# (so a -1 in the fourth place)
# Now in each pivot entry, you need to put a value to cancel what this
# -1 gives you -- so the pivot column entries are 2 and -3 (the entries
# of the fourth column that have a pivot to the left of them).
# So the first vector is [2,0,-3,-1,0], and the second is
# [4,0,5,0,-1]
# This is poorly explained (FIXME), but some examples should make it
# clear (find a good reference for this!)
*/
static GelETree *
NullSpace_op(GelCtx *ctx, GelETree * * a, gboolean *exception)
{
GelETree *n;
GelMatrixW *m;
GelMatrix *nm;
gboolean copied_m = FALSE;
int *pivot_cols, *pivot_rows;
int dim_image;
int number_of_pivots, mwh;
int i, ii, j, pi;
if G_UNLIKELY (a[0]->type == NULL_NODE)
return gel_makenum_null ();
if G_UNLIKELY ( ! check_argument_value_only_matrix (a, 0, "NullSpace"))
return NULL;
m = a[0]->mat.matrix;
if ( ! m->rref) {
m = gel_matrixw_copy (m);
gel_value_matrix_gauss (ctx, m, TRUE, FALSE, FALSE, NULL, NULL);
copied_m = TRUE;
}
dim_image = gel_matrixw_width (m);
mwh = MIN (dim_image, gel_matrixw_height (m));
pivot_cols = g_new (int, mwh);
pivot_rows = g_new (int, mwh);
number_of_pivots = get_pivot_cols (m, pivot_cols, pivot_rows);
if (dim_image == number_of_pivots) {
if (copied_m)
gel_matrixw_free (m);
g_free (pivot_cols);
g_free (pivot_rows);
return gel_makenum_null ();
}
nm = gel_matrix_new ();
gel_matrix_set_size (nm, dim_image - number_of_pivots, dim_image,
FALSE /* padding */);
j = 0;
/* Loop over nonpivots */
ii = 0;
for (i = 0; i < dim_image; i++) {
/* skip pivots */
if (ii < number_of_pivots &&
i == pivot_cols[ii]) {
ii++;
continue;
}
gel_matrix_index (nm, j, i) = gel_makenum_si (-1);
for (pi = 0; pi < number_of_pivots; pi++) {
if (pivot_cols[pi] < i) {
GelETree *t = gel_matrixw_get_index
(m, i, pivot_rows[pi]);
if (t != NULL)
gel_matrix_index (nm, j, pivot_cols[pi])
= copynode (t);
} else {
break;
}
}
j++;
}
if (copied_m)
gel_matrixw_free (m);
g_free (pivot_cols);
g_free (pivot_rows);
GET_NEW_NODE (n);
n->type = MATRIX_NODE;
n->mat.matrix = gel_matrixw_new_with_matrix (nm);
n->mat.quoted = FALSE;
return n;
}
static GelEFunc *
get_reference (GelETree *a, const char *argname, const char *func)
{
......@@ -3936,7 +4140,6 @@ DividePoly_op(GelCtx *ctx, GelETree * * a, gboolean *exception)
GelETree *n, *rn, *ql;
long size, sizeq;
int i, j;
mpw_t accu;
GelMatrixW *pm, *qm, *mn, *rem;
GelEFunc *retrem = NULL;
mpw_t tmp;
......@@ -5783,6 +5986,10 @@ gel_funclib_addall(void)
FUNC (det, 1, "M", "linear_algebra", N_("Get the determinant of a matrix"));
ALIAS (Determinant, 1, det);
FUNC (PivotColumns, 1, "M", "linear_algebra", N_("Return pivot columns of a matrix, that is columns which have a leading 1 in rref form, also returns the row where they occur"));
FUNC (NullSpace, 1, "M", "linear_algebra", N_("Get the nullspace of a matrix"))
FUNC (SetMatrixSize, 3, "M,rows,columns", "matrix", N_("Make new matrix of given size from old one"));
FUNC (IndexComplement, 2, "vec,msize", "matrix", N_("Return the index complement of a vector of indexes"));
FUNC (HermitianProduct, 2, "u,v", "matrix", N_("Get the hermitian product of two vectors"));
......
......@@ -871,4 +871,13 @@ DividePoly([1,2,3],[1,2,3],&r);r [0]
12345678911.4 12345678911.4
123456789112.4 1.23456789112e11
123456789112.6 1.23456789113e11
PivotColumns([1,1;0,1]) [1,2;1,2]
PivotColumns ([0,0,1,0;0,0,0,1]) [3,4;1,2]
PivotColumns ([0,0,1,0,0;0,0,0,1,0]) [3,4;1,2]
PivotColumns ([1,2,3;4,5,6;7,8,9]) [1,2;1,2]
PivotColumns ([1,2,3;1,2,3;1,2,3]) [1;1]
PivotColumns ([0,2,3;0,2,3;0,2,3]) [2;1]
NullSpace(null)+1 ((null)+1)
NullSpace(I(10))+1 ((null)+1)
load "nullspacetest.gel" true
load "longtest.gel" true
......@@ -2327,8 +2327,6 @@ build_program_menu (void)
int n = gtk_notebook_get_n_pages (GTK_NOTEBOOK (notebook));
int i;
GtkWidget *menu;
GtkWidget *item;
GtkWidget *w;
while (prog_menu_items != NULL) {
gtk_widget_destroy (prog_menu_items->data);
......
......@@ -377,6 +377,7 @@ gel_value_matrix_gauss (GelCtx *ctx, GelMatrixW *m, gboolean reduce, gboolean
GelETree *piv;
mpw_t tmp;
gboolean ret = TRUE;
gboolean made_private = FALSE;
if(detop)
mpw_set_ui(detop,1);
......@@ -392,6 +393,7 @@ gel_value_matrix_gauss (GelCtx *ctx, GelMatrixW *m, gboolean reduce, gboolean
! mpw_zero_p (t->val.value))
break;
}
if (j == h) {
ret = FALSE;
if(stopsing) {
......@@ -399,7 +401,14 @@ gel_value_matrix_gauss (GelCtx *ctx, GelMatrixW *m, gboolean reduce, gboolean
return FALSE;
}
continue;
} else if (j > d) {
}
if ( ! made_private) {
gel_matrixw_make_private(m);
made_private = TRUE;
}
if (j > d) {
swap_rows(m,j,d);
if(simul)
swap_rows(simul,j,d);
......@@ -407,7 +416,6 @@ gel_value_matrix_gauss (GelCtx *ctx, GelMatrixW *m, gboolean reduce, gboolean
mpw_neg(detop,detop);
}
gel_matrixw_make_private(m);
piv = gel_matrixw_index(m,i,d);
for (j = d+1; j < h; j++) {
......
......@@ -288,7 +288,10 @@ copy_internal_region (GelMatrixW *m, int w, int h)
int mi = m->regx ? m->regx[i] : i;
int mj = m->regy ? m->regy[j] : j;
GelETree *t = gel_matrix_index (old, mi, mj);
gel_matrix_index (m->m, i, j) = copynode (t);
if (t != NULL &&
(t->type != VALUE_NODE ||
! mpw_exact_zero_p (t->val.value)))
gel_matrix_index (m->m, i, j) = copynode (t);
}
}
......
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment