|
Packit |
67cb25 |
/* spgetset.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 void *tree_find(const gsl_spmatrix *m, const size_t i, const size_t j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_spmatrix_get(const gsl_spmatrix *m, const size_t i, const size_t j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (i >= m->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (j >= m->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (m->nz == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* no non-zero elements added to matrix */
|
|
Packit |
67cb25 |
return (0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (GSL_SPMATRIX_ISTRIPLET(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* traverse binary tree to search for (i,j) element */
|
|
Packit |
67cb25 |
void *ptr = tree_find(m, i, j);
|
|
Packit |
67cb25 |
double x = ptr ? *(double *) ptr : 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (GSL_SPMATRIX_ISCCS(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t *mi = m->i;
|
|
Packit |
67cb25 |
const size_t *mp = m->p;
|
|
Packit |
67cb25 |
size_t p;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* loop over column j and search for row index i */
|
|
Packit |
67cb25 |
for (p = mp[j]; p < mp[j + 1]; ++p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (mi[p] == i)
|
|
Packit |
67cb25 |
return m->data[p];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (GSL_SPMATRIX_ISCRS(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t *mj = m->i;
|
|
Packit |
67cb25 |
const size_t *mp = m->p;
|
|
Packit |
67cb25 |
size_t p;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* loop over row i and search for column index j */
|
|
Packit |
67cb25 |
for (p = mp[i]; p < mp[i + 1]; ++p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (mj[p] == j)
|
|
Packit |
67cb25 |
return m->data[p];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL("unknown sparse matrix type", GSL_EINVAL, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* element not found; return 0 */
|
|
Packit |
67cb25 |
return 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_get() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_spmatrix_set()
|
|
Packit |
67cb25 |
Add an element to a matrix in triplet form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: m - spmatrix
|
|
Packit |
67cb25 |
i - row index
|
|
Packit |
67cb25 |
j - column index
|
|
Packit |
67cb25 |
x - matrix value
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_spmatrix_set(gsl_spmatrix *m, const size_t i, const size_t j,
|
|
Packit |
67cb25 |
const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (!GSL_SPMATRIX_ISTRIPLET(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("matrix not in triplet representation", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (x == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* traverse binary tree to search for (i,j) element */
|
|
Packit |
67cb25 |
void *ptr = tree_find(m, i, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* just set the data element to 0; it would be easy to
|
|
Packit |
67cb25 |
* delete the node from the tree with avl_delete(), but
|
|
Packit |
67cb25 |
* we'd also have to delete it from the data arrays which
|
|
Packit |
67cb25 |
* is less simple
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if (ptr != NULL)
|
|
Packit |
67cb25 |
*(double *) ptr = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = GSL_SUCCESS;
|
|
Packit |
67cb25 |
void *ptr;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check if matrix needs to be realloced */
|
|
Packit |
67cb25 |
if (m->nz >= m->nzmax)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
s = gsl_spmatrix_realloc(2 * m->nzmax, m);
|
|
Packit |
67cb25 |
if (s)
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* store the triplet (i, j, x) */
|
|
Packit |
67cb25 |
m->i[m->nz] = i;
|
|
Packit |
67cb25 |
m->p[m->nz] = j;
|
|
Packit |
67cb25 |
m->data[m->nz] = x;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ptr = avl_insert(m->tree_data->tree, &m->data[m->nz]);
|
|
Packit |
67cb25 |
if (ptr != NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* found duplicate entry (i,j), replace with new x */
|
|
Packit |
67cb25 |
*((double *) ptr) = x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* no duplicate (i,j) found, update indices as needed */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* increase matrix dimensions if needed */
|
|
Packit |
67cb25 |
m->size1 = GSL_MAX(m->size1, i + 1);
|
|
Packit |
67cb25 |
m->size2 = GSL_MAX(m->size2, j + 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
++(m->nz);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_set() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double *
|
|
Packit |
67cb25 |
gsl_spmatrix_ptr(gsl_spmatrix *m, const size_t i, const size_t j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (i >= m->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("first index out of range", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (j >= m->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("second index out of range", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (GSL_SPMATRIX_ISTRIPLET(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* traverse binary tree to search for (i,j) element */
|
|
Packit |
67cb25 |
void *ptr = tree_find(m, i, j);
|
|
Packit |
67cb25 |
return (double *) ptr;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (GSL_SPMATRIX_ISCCS(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t *mi = m->i;
|
|
Packit |
67cb25 |
const size_t *mp = m->p;
|
|
Packit |
67cb25 |
size_t p;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* loop over column j and search for row index i */
|
|
Packit |
67cb25 |
for (p = mp[j]; p < mp[j + 1]; ++p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (mi[p] == i)
|
|
Packit |
67cb25 |
return &(m->data[p]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (GSL_SPMATRIX_ISCRS(m))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t *mj = m->i;
|
|
Packit |
67cb25 |
const size_t *mp = m->p;
|
|
Packit |
67cb25 |
size_t p;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* loop over row i and search for column index j */
|
|
Packit |
67cb25 |
for (p = mp[i]; p < mp[i + 1]; ++p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (mj[p] == j)
|
|
Packit |
67cb25 |
return &(m->data[p]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_NULL("unknown sparse matrix type", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* element not found; return 0 */
|
|
Packit |
67cb25 |
return NULL;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_spmatrix_ptr() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
tree_find()
|
|
Packit |
67cb25 |
Find node in tree corresponding to matrix entry (i,j). Adapted
|
|
Packit |
67cb25 |
from avl_find()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: m - spmatrix
|
|
Packit |
67cb25 |
i - row index
|
|
Packit |
67cb25 |
j - column index
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: pointer to tree node data if found, NULL if not found
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void *
|
|
Packit |
67cb25 |
tree_find(const gsl_spmatrix *m, const size_t i, const size_t j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const struct avl_table *tree = (struct avl_table *) m->tree_data->tree;
|
|
Packit |
67cb25 |
const struct avl_node *p;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (p = tree->avl_root; p != NULL; )
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t n = (double *) p->avl_data - m->data;
|
|
Packit |
67cb25 |
size_t pi = m->i[n];
|
|
Packit |
67cb25 |
size_t pj = m->p[n];
|
|
Packit |
67cb25 |
int cmp = gsl_spmatrix_compare_idx(i, j, pi, pj);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (cmp < 0)
|
|
Packit |
67cb25 |
p = p->avl_link[0];
|
|
Packit |
67cb25 |
else if (cmp > 0)
|
|
Packit |
67cb25 |
p = p->avl_link[1];
|
|
Packit |
67cb25 |
else /* |cmp == 0| */
|
|
Packit |
67cb25 |
return p->avl_data;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return NULL;
|
|
Packit |
67cb25 |
} /* tree_find() */
|