RTOpPack: Extra C/C++ Code for Vector Reduction/Transformation Operators Version of the Day
RTOp_TOp_set_sub_vector.c
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00006 //                  Copyright (2003) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 */
00030 
00031 #include "RTOp_TOp_set_sub_vector.h"
00032 #include "RTOp_obj_null_vtbl.h"
00033 
00034 #include <stdlib.h>
00035 
00036 /* Operator object data virtual function table */
00037 
00038 struct RTOp_TOp_set_sub_vector_state_t { /* operator object instance data */
00039   struct RTOp_SparseSubVector  sub_vec;   /* The sub vector to set! */
00040   int                    owns_mem;  /* if true then memory is sub_vec must be deleted! */
00041 };
00042 
00043 static int get_op_type_num_entries(
00044   const struct RTOp_obj_type_vtbl_t*   vtbl
00045   ,const void* obj_data
00046   ,int* num_values
00047   ,int* num_indexes
00048   ,int* num_chars
00049   )
00050 {
00051   const struct RTOp_TOp_set_sub_vector_state_t *state = NULL;
00052 #ifdef RTOp_DEBUG
00053   assert( num_values );
00054   assert( num_indexes );
00055   assert( num_chars );
00056   assert( obj_data );
00057 #endif
00058   state = (const struct RTOp_TOp_set_sub_vector_state_t*)obj_data;
00059   *num_values  = state->sub_vec.sub_nz;  /* values[] */
00060   *num_indexes =
00061     8  /* global_offset, sub_dim, sub_nz, values_stride, indices_stride, local_offset, is_sorted, owns_mem */
00062     + (state->sub_vec.indices ? state->sub_vec.sub_nz : 0 ); /* indices[] */
00063   *num_chars   = 0;
00064   return 0;
00065 }
00066 
00067 static int op_create(
00068   const struct RTOp_obj_type_vtbl_t* vtbl, const void* instance_data
00069   , RTOp_ReductTarget* obj )
00070 {
00071   struct RTOp_TOp_set_sub_vector_state_t *state = NULL;
00072   *obj = malloc(sizeof(struct RTOp_TOp_set_sub_vector_state_t));
00073   state = (struct RTOp_TOp_set_sub_vector_state_t*)*obj;
00074   RTOp_sparse_sub_vector_null( &state->sub_vec );
00075   state->owns_mem = 0;
00076   return 0;
00077 }
00078 
00079 static int op_free(
00080   const struct RTOp_obj_type_vtbl_t* vtbl, const void* dummy
00081   , void ** obj_data )
00082 {
00083   struct RTOp_TOp_set_sub_vector_state_t *state = NULL;
00084 #ifdef RTOp_DEBUG
00085   assert( *obj_data );
00086 #endif
00087   state = (struct RTOp_TOp_set_sub_vector_state_t*)*obj_data;
00088   if( state->owns_mem ) {
00089     free( (void*)state->sub_vec.values );
00090     free( (void*)state->sub_vec.indices );
00091   }
00092   free( *obj_data );
00093   *obj_data = NULL;
00094   return 0;
00095 }
00096 
00097 static int extract_op_state(
00098   const struct RTOp_obj_type_vtbl_t*   vtbl
00099   ,const void *       dummy
00100   ,void *             obj_data
00101   ,int                num_values
00102   ,RTOp_value_type    value_data[]
00103   ,int                num_indexes
00104   ,RTOp_index_type    index_data[]
00105   ,int                num_chars
00106   ,RTOp_char_type     char_data[]
00107   )
00108 {
00109   const struct RTOp_TOp_set_sub_vector_state_t *state = NULL;
00110   register RTOp_index_type k, j;
00111 #ifdef RTOp_DEBUG
00112   assert( obj_data );
00113 #endif
00114   state = (const struct RTOp_TOp_set_sub_vector_state_t*)obj_data;
00115 #ifdef RTOp_DEBUG
00116   assert( num_values == state->sub_vec.sub_nz );
00117   assert( num_indexes == 8 + (state->sub_vec.indices ? state->sub_vec.sub_nz : 0 ) );
00118   assert( num_chars == 0 );
00119   assert( value_data );
00120   assert( index_data );
00121 #endif
00122   for( k = 0; k < state->sub_vec.sub_nz; ++k )
00123     value_data[k] = state->sub_vec.values[k];
00124   index_data[k=0] = state->sub_vec.global_offset;
00125   index_data[++k] = state->sub_vec.sub_dim;
00126   index_data[++k] = state->sub_vec.sub_nz;
00127   index_data[++k] = state->sub_vec.values_stride;
00128   index_data[++k] = state->sub_vec.indices_stride;
00129   index_data[++k] = state->sub_vec.local_offset;
00130   index_data[++k] = state->sub_vec.is_sorted;
00131   index_data[++k] = state->owns_mem;
00132   if( state->sub_vec.indices ) {
00133     for( j = 0; j < state->sub_vec.sub_nz; ++j )
00134       index_data[++k] = state->sub_vec.indices[j];
00135   }
00136   return 0;
00137 }
00138 
00139 static int load_op_state(
00140   const struct RTOp_obj_type_vtbl_t*   vtbl
00141   ,const void *            dummy
00142   ,int                     num_values
00143   ,const RTOp_value_type   value_data[]
00144   ,int                     num_indexes
00145   ,const RTOp_index_type   index_data[]
00146   ,int                     num_chars
00147   ,const RTOp_char_type    char_data[]
00148   ,void **                 obj_data
00149   )
00150 {
00151   struct RTOp_TOp_set_sub_vector_state_t *state = NULL;
00152   register RTOp_index_type k, j;
00153   RTOp_value_type *values = NULL;
00154   /* Allocate the operator object's state data if it has not been already */
00155   if( *obj_data == NULL ) {
00156     *obj_data = (void*)malloc(sizeof(struct RTOp_TOp_set_sub_vector_state_t));
00157     state = (struct RTOp_TOp_set_sub_vector_state_t*)*obj_data;
00158     RTOp_sparse_sub_vector_null( &state->sub_vec );
00159     state->owns_mem = 0;
00160   }
00161   /* Get the operator object's state data */
00162 #ifdef RTOp_DEBUG
00163   assert( *obj_data );
00164 #endif
00165   state = (struct RTOp_TOp_set_sub_vector_state_t*)*obj_data;
00166 #ifdef RTOp_DEBUG
00167   if( num_indexes > 8 ) assert( num_values == num_indexes - 8 );
00168 #endif
00169   /* Reallocate storage if we have to */
00170   if( num_values != state->sub_vec.sub_nz || !state->owns_mem ) {
00171     /* Delete current storage if owned */
00172     if( state->owns_mem ) {
00173       free( (RTOp_value_type*)state->sub_vec.values );
00174       free( (RTOp_index_type*)state->sub_vec.indices );
00175     }
00176     /* We need to reallocate storage for values[] and perhaps */
00177     state->sub_vec.values = (RTOp_value_type*)malloc(num_values*sizeof(RTOp_value_type));
00178     if( num_indexes > 8 )
00179       state->sub_vec.indices = (RTOp_index_type*)malloc(num_values*sizeof(RTOp_index_type));
00180     state->owns_mem = 1;
00181   }
00182   /* Set the internal state! */
00183 #ifdef RTOp_DEBUG
00184   assert( num_chars == 0 );
00185   assert( value_data );
00186   assert( index_data );
00187 #endif
00188   for( values = (RTOp_value_type*)state->sub_vec.values, k = 0; k < num_values; ++k )
00189     *values++ = value_data[k];
00190   state->sub_vec.global_offset  = index_data[k=0];
00191   state->sub_vec.sub_dim        = index_data[++k];
00192   state->sub_vec.sub_nz         = index_data[++k];
00193   state->sub_vec.values_stride  = index_data[++k];
00194   state->sub_vec.indices_stride = index_data[++k];
00195   state->sub_vec.local_offset   = index_data[++k];
00196   state->sub_vec.is_sorted      = index_data[++k];
00197   state->owns_mem               = index_data[++k];
00198   if( num_indexes > 8 ) {
00199     for( j = 0; j < num_values; ++j )
00200       ((RTOp_index_type*)state->sub_vec.indices)[j] = index_data[++k];
00201   }
00202   return 0;
00203 }
00204 
00205 static struct RTOp_obj_type_vtbl_t  instance_obj_vtbl =
00206 {
00207   get_op_type_num_entries
00208   ,op_create
00209   ,NULL
00210   ,op_free
00211   ,extract_op_state
00212   ,load_op_state
00213 };
00214 
00215 /* Implementation functions */
00216 
00217 static int RTOp_TOp_set_sub_vector_apply_op(
00218   const struct RTOp_RTOp_vtbl_t* vtbl, const void* obj_data
00219   , const int num_vecs, const struct RTOp_SubVector vecs[]
00220   , const int num_targ_vecs, const struct RTOp_MutableSubVector targ_vecs[]
00221   , RTOp_ReductTarget targ_obj )
00222 {
00223   const struct RTOp_TOp_set_sub_vector_state_t *state = NULL;
00224   RTOp_index_type        v_global_offset;
00225   RTOp_index_type        v_sub_dim;
00226   RTOp_index_type        v_sub_nz;
00227   const RTOp_value_type  *v_val;
00228   ptrdiff_t              v_val_s;
00229   const RTOp_index_type  *v_ind;
00230   ptrdiff_t              v_ind_s;
00231   ptrdiff_t              v_l_off;
00232   int                    v_sorted;
00233   RTOp_index_type        z_global_offset;
00234   RTOp_index_type        z_sub_dim;
00235   RTOp_value_type        *z_val;
00236   ptrdiff_t              z_val_s;
00237   register RTOp_index_type  k, i;
00238   RTOp_index_type num_overlap;
00239 
00240   /* */
00241   /* Validate the input */
00242   /* */
00243   if( num_vecs != 0 )
00244     return RTOp_ERR_INVALID_NUM_VECS;
00245   if( num_targ_vecs != 1 )
00246     return RTOp_ERR_INVALID_NUM_TARG_VECS;
00247   assert(targ_vecs);
00248   assert(targ_obj == RTOp_REDUCT_OBJ_NULL);
00249 
00250   /* */
00251   /* Get pointers to data */
00252   /* */
00253 
00254   /* Get the sub-vector we are reading from */
00255   assert(obj_data);
00256   state = (const struct RTOp_TOp_set_sub_vector_state_t*)obj_data;
00257 
00258   /* v (the vector to read from) */
00259   v_global_offset  = state->sub_vec.global_offset;
00260   v_sub_dim        = state->sub_vec.sub_dim;
00261   v_sub_nz         = state->sub_vec.sub_nz;
00262   v_val            = state->sub_vec.values;
00263   v_val_s          = state->sub_vec.values_stride;
00264   v_ind            = state->sub_vec.indices;
00265   v_ind_s          = state->sub_vec.indices_stride;
00266   v_l_off          = state->sub_vec.local_offset;
00267   v_sorted         = state->sub_vec.is_sorted;
00268 
00269   /* z (the vector to set) */
00270   z_global_offset  = targ_vecs[0].global_offset;
00271   z_sub_dim        = targ_vecs[0].sub_dim;
00272   z_val            = targ_vecs[0].values;
00273   z_val_s          = targ_vecs[0].values_stride;
00274 
00275   /* */
00276   /* Set part of the sub-vector for this chunk. */
00277   /* */
00278 
00279   if( v_global_offset + v_sub_dim < z_global_offset + 1
00280     || z_global_offset + z_sub_dim < v_global_offset + 1 )
00281     return 0; /* The sub-vector that we are setting does not overlap with this vector chunk! */
00282 
00283   if( v_sub_nz == 0 )
00284     return 0; /* The sparse sub-vector we are reading from is empty? */
00285 
00286   /* Get the number of elements that overlap */
00287   if( v_global_offset <= z_global_offset ) {
00288     if( v_global_offset + v_sub_dim >= z_global_offset + z_sub_dim )
00289       num_overlap = z_sub_dim;
00290     else
00291       num_overlap = (v_global_offset + v_sub_dim) - z_global_offset;
00292   }
00293   else {
00294     if( z_global_offset + z_sub_dim >= v_global_offset + v_sub_dim )
00295       num_overlap = v_sub_dim;
00296     else
00297       num_overlap = (z_global_offset + z_sub_dim) - v_global_offset;
00298   }
00299 
00300   /* Set the part of the sub-vector that overlaps */
00301   if( v_ind != NULL ) {
00302     /* Sparse elements */
00303     /* Set the overlapping elements to zero first. */
00304     if( v_global_offset >= z_global_offset )
00305       z_val += (v_global_offset - z_global_offset) * z_val_s;
00306     for( k = 0; k < num_overlap; ++k, z_val += z_val_s )
00307       *z_val = 0.0;
00308     /* Now set the sparse entries */
00309     z_val = targ_vecs[0].values;
00310     for( k = 0; k < v_sub_nz; ++k, v_val += v_val_s, v_ind += v_ind_s ) {
00311       i = v_global_offset + v_l_off + (*v_ind);
00312       if( z_global_offset < i && i <= z_global_offset + z_sub_dim )
00313         z_val[ z_val_s * (i - z_global_offset - 1) ] = *v_val;
00314     }
00315     /* ToDo: Implement a faster version for v sorted and eliminate the */
00316     /* if statement in the loop. */
00317   }
00318   else {
00319     /* Dense elemements */
00320     if( v_global_offset <= z_global_offset )
00321       v_val += (z_global_offset - v_global_offset) * v_val_s;
00322     else
00323       z_val += (v_global_offset - z_global_offset) * z_val_s;
00324     for( k = 0; k < num_overlap; ++k, v_val += v_val_s, z_val += z_val_s )
00325       *z_val = *v_val;
00326   }
00327 
00328   return 0; /* success? */
00329 }
00330 
00331 /* Public interface */
00332 
00333 const struct RTOp_RTOp_vtbl_t RTOp_TOp_set_sub_vector_vtbl =
00334 {
00335   &instance_obj_vtbl
00336   ,&RTOp_obj_null_vtbl /* use null type for reduction target object */
00337   ,"RTOp_TOp_set_sub_vector"
00338   ,NULL /* use default from reduct_vtbl */
00339   ,RTOp_TOp_set_sub_vector_apply_op
00340   ,NULL /* reduce_reduct_obj */
00341   ,NULL /* get_reduct_op */
00342 };
00343 
00344 int RTOp_TOp_set_sub_vector_construct(
00345   const struct RTOp_SparseSubVector* sub_vec, struct RTOp_RTOp* op )
00346 {
00347 #ifdef RTOp_DEBUG
00348   assert(sub_vec);
00349   assert(op);
00350 #endif
00351   op->vtbl     = &RTOp_TOp_set_sub_vector_vtbl;
00352   op->vtbl->obj_data_vtbl->obj_create(NULL,NULL,&op->obj_data);
00353   return RTOp_TOp_set_sub_vector_set_sub_vec(sub_vec,op);
00354 }
00355 
00356 int RTOp_TOp_set_sub_vector_set_sub_vec(
00357   const struct RTOp_SparseSubVector* sub_vec, struct RTOp_RTOp* op )
00358 {
00359   struct RTOp_TOp_set_sub_vector_state_t *state = NULL;
00360 #ifdef RTOp_DEBUG
00361   assert( op->obj_data );
00362 #endif
00363   state = (struct RTOp_TOp_set_sub_vector_state_t*)op->obj_data;
00364   if( state->owns_mem ) {
00365     free( (void*)state->sub_vec.values );
00366     free( (void*)state->sub_vec.indices );
00367   }
00368   state->sub_vec = *sub_vec;  /* We do not own the arrays values[] and indices[] */
00369   state->owns_mem = 0;
00370   return 0;
00371 }
00372 
00373 int RTOp_TOp_set_sub_vector_destroy( struct RTOp_RTOp* op )
00374 {
00375   op_free(NULL,NULL,&op->obj_data);
00376   op->vtbl = NULL;
00377   return 0;
00378 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends