|
Packit |
67cb25 |
/* spmatrix.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2012 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_spmatrix.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "avl.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int compare_triplet(const void *pa, const void *pb, void *param);
|
|
Packit |
67cb25 |
static void *avl_spmalloc (size_t size, void *param);
|
|
Packit |
67cb25 |
static void avl_spfree (void *block, void *param);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static struct libavl_allocator avl_allocator_spmatrix =
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
avl_spmalloc,
|
|
Packit |
67cb25 |
avl_spfree
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_spmatrix_alloc()
|
|
Packit |
67cb25 |
Allocate a sparse matrix in triplet representation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: n1 - number of rows
|
|
Packit |
67cb25 |
n2 - number of columns
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: if (n1,n2) are not known at allocation time, they can each be
|
|
Packit |
67cb25 |
set to 1, and they will be expanded as elements are added to the matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_spmatrix *
|
|
Packit |
67cb25 |
gsl_spmatrix_alloc(const size_t n1, const size_t n2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double density = 0.1; /* estimate */
|
|
Packit |
67cb25 |
size_t nzmax = (size_t) floor(n1 * n2 * density);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (nzmax == 0)
|
|
Packit |
67cb25 |
nzmax = 10;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return gsl_spmatrix_alloc_nzmax(n1, n2, nzmax, GSL_SPMATRIX_TRIPLET);
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_alloc() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_spmatrix_alloc_nzmax()
|
|
Packit |
67cb25 |
Allocate a sparse matrix with given nzmax
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: n1 - number of rows
|
|
Packit |
67cb25 |
n2 - number of columns
|
|
Packit |
67cb25 |
nzmax - maximum number of matrix elements
|
|
Packit |
67cb25 |
sptype - type of matrix (triplet, CCS, CRS)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: if (n1,n2) are not known at allocation time, they can each be
|
|
Packit |
67cb25 |
set to 1, and they will be expanded as elements are added to the matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_spmatrix *
|
|
Packit |
67cb25 |
gsl_spmatrix_alloc_nzmax(const size_t n1, const size_t n2,
|
|
Packit |
67cb25 |
const size_t nzmax, const size_t sptype)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix *m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n1 == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("matrix dimension n1 must be positive integer",
|
|
Packit |
67cb25 |
GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (n2 == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL ("matrix dimension n2 must be positive integer",
|
|
Packit |
67cb25 |
GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = calloc(1, sizeof(gsl_spmatrix));
|
|
Packit |
67cb25 |
if (!m)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for spmatrix struct",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->size1 = n1;
|
|
Packit |
67cb25 |
m->size2 = n2;
|
|
Packit |
67cb25 |
m->nz = 0;
|
|
Packit |
67cb25 |
m->nzmax = GSL_MAX(nzmax, 1);
|
|
Packit |
67cb25 |
m->sptype = sptype;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->i = malloc(m->nzmax * sizeof(size_t));
|
|
Packit |
67cb25 |
if (!m->i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix_free(m);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for row indices",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (sptype == GSL_SPMATRIX_TRIPLET)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
m->tree_data = malloc(sizeof(gsl_spmatrix_tree));
|
|
Packit |
67cb25 |
if (!m->tree_data)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix_free(m);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for AVL tree struct",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->tree_data->n = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* allocate tree data structure */
|
|
Packit |
67cb25 |
m->tree_data->tree = avl_create(compare_triplet, (void *) m,
|
|
Packit |
67cb25 |
&avl_allocator_spmatrix);
|
|
Packit |
67cb25 |
if (!m->tree_data->tree)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix_free(m);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for AVL tree",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* preallocate nzmax tree nodes */
|
|
Packit |
67cb25 |
m->tree_data->node_array = malloc(m->nzmax * sizeof(struct avl_node));
|
|
Packit |
67cb25 |
if (!m->tree_data->node_array)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix_free(m);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for AVL tree nodes",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->p = malloc(m->nzmax * sizeof(size_t));
|
|
Packit |
67cb25 |
if (!m->p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix_free(m);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for column indices",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (sptype == GSL_SPMATRIX_CCS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
m->p = malloc((n2 + 1) * sizeof(size_t));
|
|
Packit |
67cb25 |
m->work = malloc(GSL_MAX(n1, n2) *
|
|
Packit |
67cb25 |
GSL_MAX(sizeof(size_t), sizeof(double)));
|
|
Packit |
67cb25 |
if (!m->p || !m->work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix_free(m);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for column pointers",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (sptype == GSL_SPMATRIX_CRS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
m->p = malloc((n1 + 1) * sizeof(size_t));
|
|
Packit |
67cb25 |
m->work = malloc(GSL_MAX(n1, n2) *
|
|
Packit |
67cb25 |
GSL_MAX(sizeof(size_t), sizeof(double)));
|
|
Packit |
67cb25 |
if (!m->p || !m->work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix_free(m);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for row pointers",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->data = malloc(m->nzmax * sizeof(double));
|
|
Packit |
67cb25 |
if (!m->data)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix_free(m);
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("failed to allocate space for data",
|
|
Packit |
67cb25 |
GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return m;
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_alloc_nzmax() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_spmatrix_free()
|
|
Packit |
67cb25 |
Free sparse matrix object
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_spmatrix_free(gsl_spmatrix *m)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (m->i)
|
|
Packit |
67cb25 |
free(m->i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (m->p)
|
|
Packit |
67cb25 |
free(m->p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (m->data)
|
|
Packit |
67cb25 |
free(m->data);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (m->work)
|
|
Packit |
67cb25 |
free(m->work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (m->tree_data)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (m->tree_data->tree)
|
|
Packit |
67cb25 |
avl_destroy(m->tree_data->tree, NULL);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (m->tree_data->node_array)
|
|
Packit |
67cb25 |
free(m->tree_data->node_array);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
free(m->tree_data);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
free(m);
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_free() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_spmatrix_realloc()
|
|
Packit |
67cb25 |
As elements are added to the sparse matrix, its possible that they
|
|
Packit |
67cb25 |
will exceed the previously specified nzmax - reallocate the matrix
|
|
Packit |
67cb25 |
with a new nzmax
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_spmatrix_realloc(const size_t nzmax, gsl_spmatrix *m)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = GSL_SUCCESS;
|
|
Packit |
67cb25 |
void *ptr;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (nzmax < m->nz)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("new nzmax is less than current nz", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ptr = realloc(m->i, nzmax * sizeof(size_t));
|
|
Packit |
67cb25 |
if (!ptr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("failed to allocate space for row indices", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->i = (size_t *) ptr;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (GSL_SPMATRIX_ISTRIPLET(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ptr = realloc(m->p, nzmax * sizeof(size_t));
|
|
Packit |
67cb25 |
if (!ptr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("failed to allocate space for column indices", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->p = (size_t *) ptr;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ptr = realloc(m->data, nzmax * sizeof(double));
|
|
Packit |
67cb25 |
if (!ptr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("failed to allocate space for data", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->data = (double *) ptr;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* rebuild binary tree */
|
|
Packit |
67cb25 |
if (GSL_SPMATRIX_ISTRIPLET(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* reset tree to empty state, but don't free root tree ptr */
|
|
Packit |
67cb25 |
avl_empty(m->tree_data->tree, NULL);
|
|
Packit |
67cb25 |
m->tree_data->n = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ptr = realloc(m->tree_data->node_array, nzmax * sizeof(struct avl_node));
|
|
Packit |
67cb25 |
if (!ptr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("failed to allocate space for AVL tree nodes", GSL_ENOMEM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m->tree_data->node_array = ptr;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* need to reinsert all tree elements since the m->data addresses
|
|
Packit |
67cb25 |
* have changed
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
for (n = 0; n < m->nz; ++n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ptr = avl_insert(m->tree_data->tree, &m->data[n]);
|
|
Packit |
67cb25 |
if (ptr != NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("detected duplicate entry", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* update to new nzmax */
|
|
Packit |
67cb25 |
m->nzmax = nzmax;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_realloc() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_spmatrix_set_zero(gsl_spmatrix *m)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
m->nz = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (GSL_SPMATRIX_ISTRIPLET(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* reset tree to empty state and node index pointer to 0 */
|
|
Packit |
67cb25 |
avl_empty(m->tree_data->tree, NULL);
|
|
Packit |
67cb25 |
m->tree_data->n = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_set_zero() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t
|
|
Packit |
67cb25 |
gsl_spmatrix_nnz(const gsl_spmatrix *m)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return m->nz;
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_nnz() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_spmatrix_compare_idx()
|
|
Packit |
67cb25 |
Comparison function for searching binary tree in triplet
|
|
Packit |
67cb25 |
representation.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
To detect duplicate elements in the tree, we want to determine
|
|
Packit |
67cb25 |
if there already exists an entry for (i,j) in the tree. Since
|
|
Packit |
67cb25 |
the actual tree node stores only the data elements data[n],
|
|
Packit |
67cb25 |
we will do pointer arithmetic to get from the given data[n]
|
|
Packit |
67cb25 |
to the row/column indices i[n] and j[n].
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This compare function will sort the tree first by row i,
|
|
Packit |
67cb25 |
and for equal rows, it will then sort by column j
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: ia - row index of element a
|
|
Packit |
67cb25 |
ja - column index of element a
|
|
Packit |
67cb25 |
ib - row index of element b
|
|
Packit |
67cb25 |
jb - column index of element b
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return:
|
|
Packit |
67cb25 |
-1 if pa < pb: (ia,ja) < (ib,jb)
|
|
Packit |
67cb25 |
+1 if pa > pb: (ia,ja) > (ib,jb)
|
|
Packit |
67cb25 |
0 if pa = pb: (ia,ja) == (ib,jb)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_spmatrix_compare_idx(const size_t ia, const size_t ja,
|
|
Packit |
67cb25 |
const size_t ib, const size_t jb)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (ia < ib)
|
|
Packit |
67cb25 |
return -1;
|
|
Packit |
67cb25 |
else if (ia > ib)
|
|
Packit |
67cb25 |
return 1;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* row indices are equal, sort by column index */
|
|
Packit |
67cb25 |
if (ja < jb)
|
|
Packit |
67cb25 |
return -1;
|
|
Packit |
67cb25 |
else if (ja > jb)
|
|
Packit |
67cb25 |
return 1;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return 0; /* row and column indices are equal */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_spmatrix_tree_rebuild()
|
|
Packit |
67cb25 |
When reading a triplet matrix from disk, or when
|
|
Packit |
67cb25 |
copying a triplet matrix, it is necessary to rebuild the
|
|
Packit |
67cb25 |
binary tree for element searches.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: m - triplet matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_spmatrix_tree_rebuild(gsl_spmatrix * m)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (!GSL_SPMATRIX_ISTRIPLET(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("m must be in triplet format", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* reset tree to empty state, but don't free root tree ptr */
|
|
Packit |
67cb25 |
avl_empty(m->tree_data->tree, NULL);
|
|
Packit |
67cb25 |
m->tree_data->n = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* insert all tree elements */
|
|
Packit |
67cb25 |
for (n = 0; n < m->nz; ++n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
void *ptr = avl_insert(m->tree_data->tree, &m->data[n]);
|
|
Packit |
67cb25 |
if (ptr != NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("detected duplicate entry", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
compare_triplet()
|
|
Packit |
67cb25 |
Comparison function for searching binary tree in triplet
|
|
Packit |
67cb25 |
representation.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
To detect duplicate elements in the tree, we want to determine
|
|
Packit |
67cb25 |
if there already exists an entry for (i,j) in the tree. Since
|
|
Packit |
67cb25 |
the actual tree node stores only the data elements data[n],
|
|
Packit |
67cb25 |
we will do pointer arithmetic to get from the given data[n]
|
|
Packit |
67cb25 |
to the row/column indices i[n] and j[n].
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This compare function will sort the tree first by row i,
|
|
Packit |
67cb25 |
and for equal rows, it will then sort by column j
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: pa - element 1 for comparison (double *)
|
|
Packit |
67cb25 |
pb - element 2 for comparison (double *)
|
|
Packit |
67cb25 |
param - parameter (gsl_spmatrix)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return:
|
|
Packit |
67cb25 |
-1 if pa < pb: (ia,ja) < (ib,jb)
|
|
Packit |
67cb25 |
+1 if pa > pb: (ia,ja) > (ib,jb)
|
|
Packit |
67cb25 |
0 if pa = pb: (ia,ja) == (ib,jb)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
compare_triplet(const void *pa, const void *pb, void *param)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix *m = (gsl_spmatrix *) param;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* pointer arithmetic to find indices in data array */
|
|
Packit |
67cb25 |
const size_t idxa = (const double *) pa - m->data;
|
|
Packit |
67cb25 |
const size_t idxb = (const double *) pb - m->data;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return gsl_spmatrix_compare_idx(m->i[idxa], m->p[idxa],
|
|
Packit |
67cb25 |
m->i[idxb], m->p[idxb]);
|
|
Packit |
67cb25 |
} /* compare_triplet() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void *
|
|
Packit |
67cb25 |
avl_spmalloc (size_t size, void *param)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_spmatrix *m = (gsl_spmatrix *) param;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (size != sizeof(struct avl_node))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("attemping to allocate incorrect node size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* return the next available avl_node slot; index
|
|
Packit |
67cb25 |
* m->tree_data->n keeps track of next open slot
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if (m->tree_data->n < m->nzmax)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* cast to char* for pointer arithmetic */
|
|
Packit |
67cb25 |
unsigned char *node_ptr = (unsigned char *) m->tree_data->node_array;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* offset in bytes for next node slot */
|
|
Packit |
67cb25 |
size_t offset = (m->tree_data->n)++ * sizeof(struct avl_node);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return node_ptr + offset;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* we should never get here - gsl_spmatrix_realloc() should
|
|
Packit |
67cb25 |
* be called before exceeding nzmax nodes
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("attemping to allocate tree node past nzmax", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
avl_spfree (void *block, void *param)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
(void)block;
|
|
Packit |
67cb25 |
(void)param;
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* do nothing - instead of allocating/freeing individual nodes,
|
|
Packit |
67cb25 |
* we malloc and free nzmax nodes at a time
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
}
|