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