ExampleNLPDirectRTOps.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 "Teuchos_ConfigDefs.hpp"
00032 #include "ExampleNLPDirectRTOps.h"
00033 #include "RTOp_obj_null_vtbl.h"
00034 #include "RTOp_obj_index_vtbl.h"
00035 
00036 /* Implementation for RTOp_TOp_explnlp2_c_eval */
00037 
00038 static int explnlp2_c_eval_apply_op(
00039   const struct RTOp_RTOp_vtbl_t* vtbl, const void* obj_data
00040   , const int num_vecs, const struct RTOp_SubVector vecs[]
00041   , const int num_targ_vecs, const struct RTOp_MutableSubVector targ_vecs[]
00042   , RTOp_ReductTarget targ_obj )
00043 {
00044   /* c */
00045   size_t                 sub_dim;
00046   RTOp_value_type        *c_val;
00047   ptrdiff_t              c_val_s;
00048   /* xD */
00049   const RTOp_value_type  *xD_val;
00050   ptrdiff_t              xD_val_s;
00051   /* xI */
00052   const RTOp_value_type  *xI_val;
00053   ptrdiff_t              xI_val_s;
00054 
00055   register RTOp_index_type  k;
00056 
00057   /*
00058    *Validate the input
00059    */
00060   if( num_vecs != 2 || vecs == NULL )
00061     return RTOp_ERR_INVALID_NUM_VECS;
00062   if( num_targ_vecs != 1 || targ_vecs == NULL )
00063     return RTOp_ERR_INVALID_NUM_TARG_VECS;
00064   if( targ_vecs[0].sub_dim != vecs[0].sub_dim
00065     || targ_vecs[0].sub_dim != vecs[1].sub_dim
00066     || targ_vecs[0].global_offset != vecs[0].global_offset
00067     || targ_vecs[0].global_offset != vecs[1].global_offset )
00068     return RTOp_ERR_INCOMPATIBLE_VECS;
00069 
00070   /*
00071    * Get pointers to data
00072    */
00073 
00074   /* c */
00075   sub_dim       = targ_vecs[0].sub_dim;
00076   c_val         = targ_vecs[0].values;
00077   c_val_s       = targ_vecs[0].values_stride;
00078   /* xD */
00079   xD_val         = vecs[0].values;
00080   xD_val_s       = vecs[0].values_stride;
00081   /* xI */
00082   xI_val         = vecs[1].values;
00083   xI_val_s       = vecs[1].values_stride;
00084 
00085   /*
00086    * Compute c(j) = xI(i) * ( xD(i) - 1 ) - 10 * xD(i)
00087    */
00088 
00089   if( c_val_s == 1 && xD_val_s == 1 && xI_val_s == 1 ) {
00090     /* Slightly faster loop for unit stride vectors */
00091     for( k = 0; k < sub_dim; ++k, ++xI_val )
00092       *c_val++ = (*xD_val++) * (*xI_val - 1.0) - 10.0 * (*xI_val);
00093   }
00094   else {
00095     /* More general implementation for non-unit strides */
00096     for( k = 0; k < sub_dim; ++k, c_val+=c_val_s, xD_val+=xD_val_s, xI_val+=xI_val_s )
00097       *c_val = (*xD_val) * (*xI_val - 1.0) - 10.0 * (*xI_val);
00098   }
00099 
00100   return 0; /* success? */
00101 }
00102 
00103 const struct RTOp_RTOp_vtbl_t RTOp_TOp_explnlp2_c_eval_vtbl =
00104 {
00105   &RTOp_obj_null_vtbl
00106   ,&RTOp_obj_null_vtbl
00107   ,"TOp_explnlp2_c_eval"
00108   ,NULL
00109   ,explnlp2_c_eval_apply_op
00110   ,NULL
00111   ,NULL
00112 };
00113 
00114 int RTOp_TOp_explnlp2_c_eval_construct( struct RTOp_RTOp* op )
00115 {
00116   op->obj_data  = NULL;
00117   op->vtbl      = &RTOp_TOp_explnlp2_c_eval_vtbl;
00118   return 0;
00119 }
00120 
00121 int RTOp_TOp_explnlp2_c_eval_destroy( struct RTOp_RTOp* op )
00122 {
00123   op->obj_data  = NULL;
00124   op->vtbl      = NULL;
00125   return 0;
00126 }
00127 
00128 /* Implementation for RTOp_TOp_explnlp2_calc_py_D */
00129 
00130 static int explnlp2_calc_py_D_apply_op(
00131   const struct RTOp_RTOp_vtbl_t* vtbl, const void* obj_data
00132   , const int num_vecs, const struct RTOp_SubVector vecs[]
00133   , const int num_targ_vecs, const struct RTOp_MutableSubVector targ_vecs[]
00134   , RTOp_ReductTarget targ_obj )
00135 {
00136   size_t                 sub_dim;
00137   /* xD */
00138   const RTOp_value_type  *xD_val;
00139   ptrdiff_t              xD_val_s;
00140   /* xI */
00141   const RTOp_value_type  *xI_val;
00142   ptrdiff_t              xI_val_s;
00143   /* c */
00144   const RTOp_value_type  *c_val;
00145   ptrdiff_t              c_val_s;
00146   /* d */
00147   RTOp_value_type        *d_val;
00148   ptrdiff_t              d_val_s;
00149   /* py */
00150   RTOp_value_type        *py_val;
00151   ptrdiff_t              py_val_s;
00152 
00153   register RTOp_index_type  k;
00154   int                       all_unit_stride = 0;
00155   RTOp_value_type           denom;
00156 
00157   /* task */
00158   int task;
00159   assert(obj_data);
00160   task = *(int*)obj_data;
00161   assert(0 <= task && task <= 2);
00162 
00163   /*
00164    * Validate the input
00165    */
00166   if( ( (task == 0 || task == 1) && num_vecs != 2 )
00167     || ( (task == 2) && num_vecs != 3 )
00168     || vecs == NULL )
00169     return RTOp_ERR_INVALID_NUM_VECS;
00170   if( ( (task == 0 || task == 1) && num_targ_vecs != 1 )
00171     || ( (task == 2) && num_targ_vecs != 2 )
00172     || targ_vecs == NULL )
00173     return RTOp_ERR_INVALID_NUM_TARG_VECS;
00174   if( targ_vecs[0].sub_dim != vecs[0].sub_dim
00175     || targ_vecs[0].sub_dim != vecs[1].sub_dim
00176     || ( task == 2 && ( targ_vecs[0].sub_dim != vecs[2].sub_dim ) )
00177     || ( task == 2 && ( targ_vecs[0].sub_dim != targ_vecs[1].sub_dim ) )
00178     || targ_vecs[0].global_offset != vecs[0].global_offset
00179     || targ_vecs[0].global_offset != vecs[1].global_offset
00180     || ( task == 2 && (targ_vecs[0].global_offset != vecs[2].global_offset ) )
00181     || ( task == 2 && ( targ_vecs[0].global_offset != targ_vecs[1].global_offset ) ) )
00182     return RTOp_ERR_INCOMPATIBLE_VECS;
00183 
00184   /*
00185    * Get pointers to data
00186    */
00187 
00188   sub_dim = vecs[0].sub_dim;
00189 
00190   k = 0;
00191   /* xD */
00192   xD_val         = vecs[k].values;
00193   xD_val_s       = vecs[k].values_stride;
00194   ++k;
00195   if( task == 1 || task == 2 ) {
00196     /* xI */
00197     xI_val         = vecs[k].values;
00198     xI_val_s       = vecs[k].values_stride;
00199     ++k;
00200   }
00201   if( task == 0 || task == 2 ) {
00202     /* c */
00203     c_val         = vecs[k].values;
00204     c_val_s       = vecs[k].values_stride;
00205     ++k;
00206   }
00207   k = 0;
00208   if( task == 1 || task == 2 ) {
00209     /* d */
00210     d_val         = targ_vecs[k].values;
00211     d_val_s       = targ_vecs[k].values_stride;
00212     ++k;
00213   }
00214   if( task == 0 || task == 2 ) {
00215     /* py */
00216     py_val         = targ_vecs[k].values;
00217     py_val_s       = targ_vecs[k].values_stride;
00218     ++k;
00219   }
00220 
00221   /* Determine if all the vectors have unit stride! */
00222   all_unit_stride = 1;
00223   for( k = 0; k < num_vecs && !all_unit_stride; ++k )
00224     if( vecs[k].values_stride != 1 )
00225       all_unit_stride = 0;
00226   for( k = 0; k < num_targ_vecs && !all_unit_stride; ++k )
00227     if( targ_vecs[k].values_stride != 1 )
00228       all_unit_stride = 0;
00229 
00230   /*
00231    * Compute py and/or D
00232    */
00233 
00234   if( all_unit_stride) {
00235     if(task == 0) {
00236       /* Compute py only */
00237       for( k = 0; k < sub_dim; ++k )
00238         *py_val++ = *c_val++ / ( 1.0 - *xI_val++ );
00239     }
00240     if(task == 1) {
00241       /* Compute D only */
00242       for( k = 0; k < sub_dim; ++k )
00243         *d_val++ = ( *xD_val++ - 10.0 ) / ( 1.0 - *xI_val++ );
00244     }
00245     if(task == 2) {
00246       /* Compute py and D */
00247       for( k = 0; k < sub_dim; ++k ) {
00248         denom = ( 1.0 - *xI_val++ );
00249         *d_val++ = ( *xD_val++ - 10.0 ) / denom;
00250         *py_val++ = *c_val++ / denom;
00251       }
00252     }
00253   }
00254   else {
00255     assert(0); /* ToDo: Implement if needed! */
00256   }
00257 
00258   return 0;
00259 }
00260 
00261 const struct RTOp_RTOp_vtbl_t RTOp_TOp_explnlp2_calc_py_D_vtbl =
00262 {
00263   &RTOp_obj_index_vtbl
00264   ,&RTOp_obj_null_vtbl
00265   ,"TOp_explnlp2_calc_py_D"
00266   ,NULL
00267   ,explnlp2_calc_py_D_apply_op
00268   ,NULL
00269   ,NULL
00270 };
00271 
00272 int RTOp_TOp_explnlp2_calc_py_D_construct( int task, struct RTOp_RTOp* op )
00273 {
00274   int result;
00275 #ifdef RTOp_DEBUG
00276   assert( 0 <= task && task <= 2 );
00277 #endif  
00278   op->obj_data  = NULL;
00279   op->vtbl      = &RTOp_TOp_explnlp2_calc_py_D_vtbl;
00280   result = op->vtbl->obj_data_vtbl->obj_create( NULL, NULL, &op->obj_data );
00281   if(result != 0) return result;
00282 #ifdef RTOp_DEBUG
00283   assert(op->obj_data);
00284 #endif
00285   *((int*)op->obj_data) = task;
00286   return 0;
00287 }
00288 
00289 int RTOp_TOp_explnlp2_calc_py_D_set_task( int task, struct RTOp_RTOp* op )
00290 {
00291 #ifdef RTOp_DEBUG
00292   assert( 0 <= task && task <= 2 );
00293   assert(op->obj_data);
00294 #endif
00295   *((int*)op->obj_data) = task;
00296   return 0;
00297 }
00298 
00299 int RTOp_TOp_explnlp2_calc_py_D_destroy( struct RTOp_RTOp* op )
00300 {
00301   int result;
00302   result = op->vtbl->reduct_vtbl->obj_free( NULL, NULL, &op->obj_data );
00303   if(result != 0) return result;
00304   op->obj_data  = NULL;
00305   op->vtbl      = NULL;
00306   return 0;
00307 }

Generated on Tue Jul 13 09:36:22 2010 for MOOCHO by  doxygen 1.4.7