ADMB Documentation  11.1.2432
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m18.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar_m18.cpp 1919 2014-04-22 22:02:01Z 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 void cmdm_prod(void);
00023 
00028 dvar_matrix operator*(const dmatrix& cm1, const dvar_matrix& m2)
00029  {
00030    if (cm1.colmin() != m2.rowmin() || cm1.colmax() != m2.rowmax())
00031    {
00032      cerr << " Incompatible array bounds in "
00033      "dmatrix operator*(const dmatrix& x, const dvar_matrix& m)\n";
00034      ad_exit(21);
00035    }
00036    dmatrix cm2=value(m2);
00037    dmatrix tmp(cm1.rowmin(),cm1.rowmax(), m2.colmin(), m2.colmax());
00038    double sum;
00039    double * temp_col=(double *) malloc(m2.rowsize()*sizeof(double));
00040    temp_col-=cm2.rowmin();
00041 
00042    for (int j=cm2.colmin(); j<=cm2.colmax(); j++)
00043    {
00044      for (int k=cm2.rowmin(); k<=cm2.rowmax(); k++)
00045      {
00046        temp_col[k] = cm2.elem(k,j);
00047      }
00048 
00049      for (int i=cm1.rowmin(); i<=cm1.rowmax(); i++)
00050      {
00051        sum=0.0;
00052        const dvector& temp_row = cm1(i);
00053        for (int k=cm1.colmin(); k<=cm1.colmax(); k++)
00054        {
00055           sum+=temp_row(k) * (temp_col[k]);
00056          // sum+=temp_row(k) * cm2(k,j);
00057        }
00058        tmp(i,j)=sum;
00059      }
00060    }
00061 
00062    temp_col+=cm2.rowmin();
00063    free ((char*)temp_col);
00064    dvar_matrix vtmp=nograd_assign(tmp);
00065    save_identifier_string("TEST1");
00066    cm1.save_dmatrix_value();
00067    cm1.save_dmatrix_position();
00068    // m2.save_dvar_matrix_value();
00069    m2.save_dvar_matrix_position();
00070    vtmp.save_dvar_matrix_position();
00071    save_identifier_string("TEST6");
00072    gradient_structure::GRAD_STACK1->
00073             set_gradient_stack(cmdm_prod);
00074    return vtmp;
00075  }
00076 
00081 void cmdm_prod(void)
00082 {
00083   verify_identifier_string("TEST6");
00084   dvar_matrix_position vpos=restore_dvar_matrix_position();
00085   dmatrix dftmp=restore_dvar_matrix_derivatives(vpos);
00086   dvar_matrix_position m2pos=restore_dvar_matrix_position();
00087   //dmatrix cm2=restore_dvar_matrix_value(m2pos);
00088   dmatrix_position m1pos=restore_dmatrix_position();
00089   dmatrix cm1=restore_dmatrix_value(m1pos);
00090   verify_identifier_string("TEST1");
00091   //dmatrix dfm1(m1pos);
00092   dmatrix dfm2(m2pos);
00093   double dfsum;
00094   dfm2.initialize();
00095   for (int j=dfm2.colmin(); j<=dfm2.colmax(); j++)
00096   {
00097     for (int i=cm1.rowmin(); i<=cm1.rowmax(); i++)
00098     {
00099       //tmp.elem(i,j)=sum;
00100       dfsum=dftmp.elem(i,j);
00101       for (int k=cm1.colmin(); k<=cm1.colmax(); k++)
00102       {
00103         //sum+=cm1(i,k) * cm2(k,j);
00104        //dfm1.elem(i,k)+=dfsum * cm2.elem(k,j);
00105         dfm2.elem(k,j)+=dfsum * cm1.elem(i,k);
00106       }
00107     }
00108   }
00109   //dfm1.save_dmatrix_derivatives(m1pos);
00110   dfm2.save_dmatrix_derivatives(m2pos);
00111   // cout << "leaving dmdm_prod"<<endl;
00112 }