ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m18.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar_m18.cpp 2703 2014-11-20 18:26:33Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00013 #ifdef __TURBOC__
00014   #pragma hdrstop
00015   #include <iostream.h>
00016 #endif
00017 
00018 #ifdef __ZTC__
00019   #include <iostream.hpp>
00020 #endif
00021 
00022 #ifndef OPT_LIB
00023   #include <cassert>
00024 #endif
00025 
00026 void cmdm_prod(void);
00027 
00032 dvar_matrix operator*(const dmatrix& cm1, const dvar_matrix& m2)
00033  {
00034    if (cm1.colmin() != m2.rowmin() || cm1.colmax() != m2.rowmax())
00035    {
00036      cerr << " Incompatible array bounds in "
00037      "dmatrix operator*(const dmatrix& x, const dvar_matrix& m)\n";
00038      ad_exit(21);
00039    }
00040    dmatrix cm2=value(m2);
00041    dmatrix tmp(cm1.rowmin(),cm1.rowmax(), m2.colmin(), m2.colmax());
00042 
00043 #ifdef OPT_LIB
00044    const size_t rowsize = (size_t)m2.rowsize();
00045 #else
00046    const int _rowsize = m2.rowsize();
00047    assert(_rowsize > 0);
00048    const size_t rowsize = (size_t)_rowsize;
00049 #endif
00050    try
00051    {
00052      double* temp_col = new double[rowsize];
00053 
00054      temp_col-=cm2.rowmin();
00055 
00056      for (int j=cm2.colmin(); j<=cm2.colmax(); j++)
00057      {
00058        for (int k=cm2.rowmin(); k<=cm2.rowmax(); k++)
00059        {
00060          temp_col[k] = cm2.elem(k,j);
00061        }
00062        for (int i=cm1.rowmin(); i<=cm1.rowmax(); i++)
00063        {
00064          double sum=0.0;
00065          const dvector& temp_row = cm1(i);
00066          for (int k=cm1.colmin(); k<=cm1.colmax(); k++)
00067          {
00068             sum+=temp_row(k) * (temp_col[k]);
00069            // sum+=temp_row(k) * cm2(k,j);
00070          }
00071          tmp(i,j)=sum;
00072        }
00073      }
00074      temp_col+=cm2.rowmin();
00075      delete [] temp_col;
00076      temp_col = 0;
00077    }
00078    catch (std::bad_alloc& e)
00079    {
00080      cerr << "Error[" << __FILE__ << ':' << __LINE__
00081           << "]: Unable to allocate array.\n";
00082      //ad_exit(21);
00083      throw e;
00084    }
00085 
00086    dvar_matrix vtmp=nograd_assign(tmp);
00087    save_identifier_string("TEST1");
00088    cm1.save_dmatrix_value();
00089    cm1.save_dmatrix_position();
00090    // m2.save_dvar_matrix_value();
00091    m2.save_dvar_matrix_position();
00092    vtmp.save_dvar_matrix_position();
00093    save_identifier_string("TEST6");
00094    gradient_structure::GRAD_STACK1->
00095             set_gradient_stack(cmdm_prod);
00096    return vtmp;
00097  }
00098 
00103 void cmdm_prod(void)
00104 {
00105   verify_identifier_string("TEST6");
00106   dvar_matrix_position vpos=restore_dvar_matrix_position();
00107   dmatrix dftmp=restore_dvar_matrix_derivatives(vpos);
00108   dvar_matrix_position m2pos=restore_dvar_matrix_position();
00109   //dmatrix cm2=restore_dvar_matrix_value(m2pos);
00110   dmatrix_position m1pos=restore_dmatrix_position();
00111   dmatrix cm1=restore_dmatrix_value(m1pos);
00112   verify_identifier_string("TEST1");
00113   //dmatrix dfm1(m1pos);
00114   dmatrix dfm2(m2pos);
00115   double dfsum;
00116   dfm2.initialize();
00117   for (int j=dfm2.colmin(); j<=dfm2.colmax(); j++)
00118   {
00119     for (int i=cm1.rowmin(); i<=cm1.rowmax(); i++)
00120     {
00121       //tmp.elem(i,j)=sum;
00122       dfsum=dftmp.elem(i,j);
00123       for (int k=cm1.colmin(); k<=cm1.colmax(); k++)
00124       {
00125         //sum+=cm1(i,k) * cm2(k,j);
00126        //dfm1.elem(i,k)+=dfsum * cm2.elem(k,j);
00127         dfm2.elem(k,j)+=dfsum * cm1.elem(i,k);
00128       }
00129     }
00130   }
00131   //dfm1.save_dmatrix_derivatives(m1pos);
00132   dfm2.save_dmatrix_derivatives(m2pos);
00133   // cout << "leaving dmdm_prod"<<endl;
00134 }