ADMB Documentation  11.1.2438
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m20.cpp
Go to the documentation of this file.
00001 
00007 #include <fvar.hpp>
00008 #ifdef __TURBOC__
00009   #pragma hdrstop
00010   #include <iomanip.h>
00011 #endif
00012 
00013 #ifdef __ZTC__
00014   #include <iomanip.hpp>
00015 #endif
00016 
00017 #define TINY 1.0e-20;
00018 void dfinvpret(void);
00019 
00020 // int min(int a,int b);
00021 void df_xdet(void);
00022 
00033 dvariable det(const dvar_matrix& aa)
00034 {
00035   int i,j,k,n;
00036   n=aa.colsize();
00037   int lb=aa.colmin();
00038   int ub=aa.colmax();
00039   ivector indx(lb,ub);
00040   if (lb!=aa.rowmin()||ub!=aa.colmax())
00041   {
00042     cerr << "Error matrix not square in det()"<<endl;
00043     ad_exit(1);
00044   }
00045   int One=1;
00046   indx.fill_seqadd(lb,One);
00047   double d;
00048   double big,dum,sum,temp;
00049   dvar_matrix_position dmp(aa,1);
00050   dmatrix bb=value(aa);
00051   dvector vv(lb,ub);
00052   dvector part_prod(lb,ub);
00053 
00054   d=1.0;
00055   for (i=lb;i<=ub;i++)
00056   {
00057     big=0.0;
00058     for (j=lb;j<=ub;j++)
00059     {
00060       temp=fabs(bb.elem(i,j));
00061       if (temp > big)
00062       {
00063         big=temp;
00064       }
00065     }
00066     if (big == 0.0)
00067     {
00068       cerr << "Error in matrix inverse -- matrix singular in inv(dmatrix)\n";
00069     }
00070     vv[i]=1.0/big;
00071   }
00072 
00073   for (j=lb;j<=ub;j++)
00074   {
00075     for (i=lb;i<j;i++)
00076     {
00077       sum=bb.elem(i,j);
00078       for (k=lb;k<i;k++)
00079       {
00080         sum = sum - bb.elem(i,k)*bb.elem(k,j);
00081       }
00082       //a[i][j]=sum;
00083       bb(i,j)=sum;
00084     }
00085     int imax = j;
00086     big=0.0;
00087     for (i=j;i<=ub;i++)
00088     {
00089       sum=bb.elem(i,j);
00090       for (k=lb;k<j;k++)
00091       {
00092         sum = sum - bb(i,k)*bb(k,j);
00093       }
00094       bb(i,j)=sum;
00095       dum=vv.elem(i)*fabs(sum);
00096       if ( dum >= big)
00097       {
00098         big=dum;
00099         imax=i;
00100       }
00101     }
00102     if (j != imax)
00103     {
00104       for (k=lb;k<=ub;k++)
00105       {
00106         dum=bb.elem(imax,k);
00107         bb.elem(imax,k)=bb.elem(j,k);
00108         bb.elem(j,k)=dum;
00109       }
00110       d = -1.*d;
00111       vv.elem(imax)=vv.elem(j);
00112 
00113       //if (j<ub)
00114       {
00115         int itemp=indx.elem(imax);
00116         indx.elem(imax)=indx.elem(j);
00117         indx.elem(j)=itemp;
00118       }
00119       //cout << "indx= " <<indx<<endl;
00120     }
00121 
00122     if (bb.elem(j,j) == 0.0)
00123     {
00124       bb(j,j)=TINY;
00125     }
00126 
00127     if (j != n)
00128     {
00129       dum=1.0/bb(j,j);
00130       for (i=j+1;i<=ub;i++)
00131       {
00132         bb.elem(i,j) *= dum;
00133       }
00134     }
00135   }
00136   double det=d;
00137 
00138   // SM Bug 129, issue appears to be at this line
00139   // part_prod is declared above as dvector(lb,ub)
00140   // cout<<"Bug 129 ="<<part_prod(lb)<<endl;
00141   // part_prod(1)=d*bb(1,1);  // replaced this line with:
00142   part_prod(lb) = d*bb(lb,lb);
00143   // cout<<"Ok got this far; det = "<<det<<endl;
00144   for (j=lb+1;j<=ub;j++)
00145   {
00146     part_prod(j)=part_prod(j-1)*bb(j,j);
00147   }
00148   det=part_prod(ub);
00149   dvariable rdet=nograd_assign(det);
00150   save_identifier_string("PLACE7");
00151   part_prod.save_dvector_value();
00152   part_prod.save_dvector_position();
00153   indx.save_ivector_value();
00154   indx.save_ivector_position();
00155   save_identifier_string("PLACE3");
00156   aa.save_dvar_matrix_position();
00157   save_identifier_string("PLACE2b");
00158   rdet.save_prevariable_position();
00159   save_identifier_string("PLACE2a");
00160   bb.save_dmatrix_value();
00161   save_identifier_string("PLACE2");
00162   bb.save_dmatrix_position();
00163   save_identifier_string("PLACE1");
00164   save_double_value(d);
00165   save_identifier_string("PLACE0");
00166   gradient_structure::GRAD_STACK1->
00167       set_gradient_stack(df_xdet);
00168   return rdet;
00169 }
00170 
00172 void df_xdet(void)
00173 {
00174   verify_identifier_string("PLACE0");
00175   double d=restore_double_value();
00176   verify_identifier_string("PLACE1");
00177   dmatrix_position bpos=restore_dmatrix_position();
00178   verify_identifier_string("PLACE2");
00179   dmatrix b=restore_dmatrix_value(bpos);
00180   verify_identifier_string("PLACE2a");
00181   //dvar_matrix_position rdet_pos=restore_prevariable_position();
00182   double dfdet=restore_prevariable_derivative();
00183   verify_identifier_string("PLACE2b");
00184   dvar_matrix_position a_pos=restore_dvar_matrix_position();
00185   verify_identifier_string("PLACE3");
00186   ivector_position indx_pos=restore_ivector_position();
00187   ivector indx=restore_ivector_value(indx_pos);
00188   dvector_position part_prod_pos=restore_dvector_position();
00189   dvector part_prod=restore_dvector_value(part_prod_pos);
00190   verify_identifier_string("PLACE7");
00191   int lb=b.colmin();
00192   int ub=b.colmax();
00193   dmatrix dfb(lb,ub,lb,ub);
00194 
00195   dvector dfpart_prod(lb,ub);
00196 
00197   #ifndef SAFE_INITIALIZE
00198     dfb.initialize();
00199     dfpart_prod.initialize();
00200   #endif
00201 
00202 
00203   dfpart_prod(ub)=dfdet;
00204   int j;
00205   for (j=ub;j>=lb+1;j--)
00206   {
00207     // part_prod(j)=part_prod(j-1)*b(j,j);
00208     dfpart_prod(j-1)+=dfpart_prod(j)*b(j,j);
00209     dfb(j,j)+=dfpart_prod(j)*part_prod(j-1);
00210     dfpart_prod(j)=0.;
00211   }
00212   //part_prod(1)=d*b(lb,lb);
00213   dfb(lb,lb)+=dfpart_prod(lb)*d;
00214   dfpart_prod(lb)=0.;
00215 
00216   double dfsum=0.;
00217   for (j=ub;j>=lb;j--)
00218   {
00219     for (int i=ub;i>=lb;i--)
00220     {
00221       if (i<=j)
00222       {
00223         // b(i,j)=sum;
00224         dfsum+=dfb(i,j);
00225         dfb(i,j)=0.;
00226       }
00227       else
00228       {
00229         // b(i,j)=sum/b(j,j);
00230         dfsum+=dfb(i,j)/b(j,j);
00231         dfb(j,j)-=dfb(i,j)*b(i,j)/b(j,j);
00232         dfb(i,j)=0.;
00233       }
00234 
00235       for (int k=min(i-1,j-1);k>=lb;k--)
00236       {
00237         // sum-=b(i,k)*b(k,j);
00238         dfb(i,k)-=dfsum*b(k,j);
00239         dfb(k,j)-=dfsum*b(i,k);
00240       }
00241       // sum=value(a(indx(i),j);
00242       save_dmatrix_derivatives(a_pos,dfsum,indx(i),j); // like this
00243       dfsum=0.;
00244     }
00245   }
00246 }
00247 
00248 #undef TINY