/* tensor/init_source.c * * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough * * 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 2 of the License, 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., 675 Mass Ave, Cambridge, MA 02139, USA. */ /* * This code follows as close as possible that of init_source.c in * the matrix directory */ /* * We should probably call a modified version of sys/pow_int.c, but * let's put a hack to show how it should work. * * This obviously returns x^n */ size_t quick_pow(size_t x, unsigned int n) { double value = 1.0; /* * repeated squaring method */ do { if (n & 1) value *= x; /* for n odd */ n >>= 1; /* n = n/2 */ x *= x; } while (n); return value; } /* * Allocate memory for a tensor and return a pointer to it. */ TYPE (gsl_tensor) * FUNCTION (gsl_tensor, alloc) (const unsigned int rank, const size_t dimension) { TYPE (gsl_tensor) * t; if (dimension == 0) { GSL_ERROR_VAL ("tensor dimension must be positive integer", GSL_EINVAL, 0); } t = gsl_tensor * malloc (sizeof (TYPE (gsl_tensor))); if (t == 0) { GSL_ERROR_VAL ("failed to allocate space for tensor struct", GSL_ENOMEM, 0); } t->data = (ATOMIC *) malloc (MULTIPLICITY * n * quick_pow(dimension, rank) * sizeof (ATOMIC)); if (t->data == 0) GSL_ERROR_VAL ("failed to allocate space for data", GSL_ENOMEM, 0); t->rank = rank; t->dimension = dimension; return t; } /* * Same as gsl_tensor_alloc, but put all elements to 0. */ TYPE (gsl_tensor) * FUNCTION (gsl_tensor, calloc) (const unsigned int rank, const size_t dimension) { size_t i; size_t n; TYPE (gsl_tensor) * t = FUNCTION (gsl_tensor, alloc) (rank, dimension); if (t == 0) return 0; /* initialize tensor to zero */ n = MULTIPLICITY * quick_pow(dimension, rank); for (i = 0; i < n; i++) t->data[i] = 0; return t; } /* * Copy constructor, oh, I mean, copy from an existing tensor. */ TYPE (gsl_tensor) * FUNCTION (gsl_tensor, copy) (TYPE(gsl_tensor) * tt) { size_t n; TYPE (gsl_tensor) * t; t = (TYPE (gsl_tensor) *) malloc (sizeof (TYPE (gsl_tensor))); if (t == 0) GSL_ERROR_VAL ("failed to allocate space for tensor struct", GSL_ENOMEM, 0); n = MULTIPLICITY * quick_pow(tt->dimension, tt->rank); memcpy(t->data, tt->data, n); t->rank = tt->rank; t->dimension = tt->dimension; return t; } /* * Free memory. */ void FUNCTION (gsl_tensor, free) (TYPE (gsl_tensor) * t) { free(t->data); free(t); } /* Operations */ /* * t = 0 (all elements = 0) */ void FUNCTION (gsl_tensor, set_zero) (TYPE (gsl_tensor) * t) { ATOMIC * const data = t->data; size_t n; n = MULTIPLICITY * quick_pow(t->dimension, t->rank); for (i = 0; i < n; i += MULTIPLICITY) *(BASE *) (data + MULTIPLICITY * i) = zero; } /* * t_ijk = x (for i,j,k = 0,1,...,dimension-1) */ void FUNCTION (gsl_tensor, set_all) (TYPE (gsl_tensor) * t, BASE x) { size_t n; ATOMIC * const data = m->data; n = MULTIPLICITY * quick_pow(t->dimension, t->rank); for (i = 0; i < n; i += MULTIPLICITY) *(BASE *) (data + MULTIPLICITY * i) = x; }