ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m15.cpp
Go to the documentation of this file.
00001 
00007 #include <fvar.hpp>
00008 
00009 #ifdef __TURBOC__
00010   #pragma hdrstop
00011   #include <iostream.h>
00012 #endif
00013 
00014 #ifdef __ZTC__
00015   #include <iostream.hpp>
00016 #endif
00017 
00018 #define TINY 1.0e-20;
00019 void dfinvpret(void);
00020 
00026 int min(const int a, const int b)
00027 {
00028   return a <= b ? a : b;
00029 }
00030 
00038 dvar_matrix inv(const dvar_matrix& aa)
00039 {
00040   int imax = 0;
00041   int n = aa.colsize();
00042   int lb=aa.colmin();
00043   int ub=aa.colmax();
00044   dvar_matrix vc(lb,ub,lb,ub);
00045   if (n==1)
00046   {
00047     if (aa(lb,lb)==0.0)
00048     {
00049       cerr << "Error in matrix inverse -- matrix singular in inv(dmatrix)\n";
00050       ad_exit(1);
00051     }
00052     else
00053     {
00054       vc(lb,lb)=1.0/aa(lb,lb);
00055       return vc;
00056     }
00057   }
00058   ivector indx(lb,ub);
00059   int One=1;
00060   indx.fill_seqadd(lb,One);
00061   double d;
00062   double big,dum,sum,temp;
00063   dvar_matrix_position dmp(aa,0);
00064   dmatrix bb=value(aa);
00065   dvector vv(lb,ub);
00066 
00067   d=1.0;
00068   for (int i=lb;i<=ub;i++)
00069   {
00070     big=0.0;
00071     for (int j=lb;j<=ub;j++)
00072     {
00073       temp=fabs(bb.elem(i,j));
00074       if (temp > big)
00075       {
00076         big=temp;
00077       }
00078     }
00079     if (big == 0.0)
00080     {
00081       cerr << "Error in matrix inverse -- matrix singular in inv(dmatrix)\n";
00082       ad_exit(1);
00083     }
00084     vv[i]=1.0/big;
00085   }
00086 
00087   for (int j=lb;j<=ub;j++)
00088   {
00089     for (int i=lb;i<j;i++)
00090     {
00091       sum=bb.elem(i,j);
00092       for (int k=lb;k<i;k++)
00093       {
00094         sum = sum - bb.elem(i,k)*bb.elem(k,j);
00095       }
00096       //a[i][j]=sum;
00097       bb.elem(i,j)=sum;
00098     }
00099     big=0.0;
00100     for (int i=j;i<=ub;i++)
00101     {
00102       sum=bb.elem(i,j);
00103       for (int k=lb;k<j;k++)
00104       {
00105         sum = sum - bb.elem(i,k)*bb.elem(k,j);
00106       }
00107       bb.elem(i,j)=sum;
00108       dum=vv[i]*fabs(sum);
00109       if ( dum >= big)
00110       {
00111         big=dum;
00112         imax=i;
00113       }
00114     }
00115     if (j != imax)
00116     {
00117       for (int k=lb;k<=ub;k++)
00118       {
00119         dum=bb.elem(imax,k);
00120         bb.elem(imax,k)=bb.elem(j,k);
00121         bb.elem(j,k)=dum;
00122       }
00123       d = -1.*d;
00124       vv[imax]=vv[j];
00125 
00126       //if (j<ub)
00127       {
00128         int itemp=indx.elem(imax);
00129         indx.elem(imax)=indx.elem(j);
00130         indx.elem(j)=itemp;
00131       }
00132       //cout << "indx= " <<indx<<endl;
00133     }
00134 
00135     if (bb.elem(j,j) == 0.0)
00136     {
00137       bb.elem(j,j)=TINY;
00138     }
00139 
00140     if (j != n)
00141     {
00142       dum=1.0/bb.elem(j,j);
00143       for (int i=j+1;i<=ub;i++)
00144       {
00145         bb.elem(i,j) = bb.elem(i,j) * dum;
00146       }
00147     }
00148   }
00149 
00150   dvector y(lb,ub);
00151   dvector x(lb,ub);
00152   //int lb=rowmin;
00153   //int ub=rowmax;
00154   dmatrix& b=bb;
00155   ivector indxinv(lb,ub);
00156   for (int i=lb;i<=ub;i++)
00157   {
00158     indxinv(indx.elem(i))=i;
00159   }
00160   for (int ii=lb;ii<=ub;ii++)
00161   {
00162     y.initialize();
00163     y(indxinv(ii))=1.;
00164     for (int i=indxinv(ii);i<=ub;i++)
00165     {
00166       // sum=y(ii,i);
00167       if (i==indxinv(ii))
00168       {
00169         sum=1.;
00170       }
00171       else
00172       {
00173         sum=0.;
00174       }
00175       for (int j=indxinv(ii);j<=i-1;j++)
00176       {
00177         sum-=b.elem(i,j)*y.elem(j);
00178       }
00179       y.elem(i)=sum;
00180     }
00181     for (int i=ub;i>=lb;i--)
00182     {
00183       sum=y.elem(i);
00184       for (int j=i+1;j<=ub;j++)
00185       {
00186         sum-=b.elem(i,j)*x.elem(j);
00187       }
00188       x.elem(i)=sum/b.elem(i,i);
00189     }
00190     y.save_dvector_value();
00191     x.save_dvector_value();
00192     nograd_assign_column(vc,x,ii);
00193   }
00194 
00195   save_identifier_string("P5");
00196   x.save_dvector_position();
00197   y.save_dvector_position();
00198   indx.save_ivector_value();
00199   indx.save_ivector_position();
00200   aa.save_dvar_matrix_position();
00201   vc.save_dvar_matrix_position();
00202   bb.save_dmatrix_value();
00203   bb.save_dmatrix_position();
00204   save_identifier_string("P1");
00205   gradient_structure::GRAD_STACK1->
00206       set_gradient_stack(dfinvpret);
00207   return vc;
00208 }
00209 
00212 void dfinvpret(void)
00213 {
00214   verify_identifier_string("P1");
00215   dmatrix_position bpos=restore_dmatrix_position();
00216   dmatrix b=restore_dmatrix_value(bpos);
00217   dvar_matrix_position v_pos=restore_dvar_matrix_position();
00218   dvar_matrix_position a_pos=restore_dvar_matrix_position();
00219   ivector_position indx_pos=restore_ivector_position();
00220   ivector indx=restore_ivector_value(indx_pos);
00221   dvector_position y_pos=restore_dvector_position();
00222   dvector_position x_pos=restore_dvector_position();
00223   verify_identifier_string("P5");
00224   int lb=b.colmin();
00225   int ub=b.colmax();
00226   dmatrix dfb(lb,ub,lb,ub);
00227   ivector indxinv(lb,ub);
00228   for (int i=lb;i<=ub;i++)
00229   {
00230     indxinv(indx.elem(i))=i;
00231   }
00232 
00233   double dfsum=0.;
00234   dvector dfy(lb,ub);
00235   #ifndef SAFE_INITIALIZE
00236     dfb.initialize();
00237     dfy.initialize();
00238   #endif
00239   for (int ii=ub;ii>=lb;ii--)
00240   {
00241     //x.save_dvector_value();
00242     dvector x=restore_dvector_value(x_pos);
00243     //y.save_dvector_value();
00244     dvector y=restore_dvector_value(y_pos);
00245     dvector dfx=restore_dvar_matrix_derivative_column(v_pos,ii);
00246     for (int i=lb;i<=ub;i++)
00247     {
00248       // x.elem(i)=sum/b.elem(i,i);
00249       dfsum+=dfx.elem(i)/b.elem(i,i);
00250       dfb.elem(i,i)-=dfx.elem(i)*x.elem(i)/b.elem(i,i);
00251       dfx.elem(i)=0.;
00252       for (int j=ub;j>=i+1;j--)
00253       {
00254         // sum -=b.elem(i,j)*x.elem(j);
00255         dfb.elem(i,j)-=dfsum*x.elem(j);
00256         dfx.elem(j)-=dfsum*b.elem(i,j);
00257       }
00258       // sum=y.elem(i);
00259       dfy.elem(i)+=dfsum;
00260       dfsum=0.;
00261     }
00262 
00263     //for (i=ub;i>=lb;i--)
00264     int i2;
00265     for (i2=ub;i2>=indxinv(ii);i2--)
00266     {
00267       // y.elem(i)=sum;
00268       dfsum+=dfy.elem(i2);
00269       dfy.elem(i2)=0.;
00270       // for (int j=i-1;j>=lb;j--)
00271       for (int j=i2-1;j>=indxinv(ii);j--)
00272       {
00273         // sum-=b.elem(i,j)*y.elem(j);
00274         dfb.elem(i2,j)-=dfsum*y.elem(j);
00275         dfy.elem(j)-=dfsum*b.elem(i2,j);
00276       }
00277       //sum=y.elem(i);
00278       dfy.elem(i2)=dfsum;
00279       dfsum=0.;
00280     }
00281     //x.initialize()
00282     //y.initialize()
00283     dfx.initialize();
00284     dfy.initialize();
00285   }
00286 
00287   for (int j=ub;j>=lb;j--)
00288   {
00289     for (int i=ub;i>=lb;i--)
00290     {
00291       if (i<=j)
00292       {
00293         // b.elem(i,j)=sum;
00294         dfsum+=dfb.elem(i,j);
00295         dfb.elem(i,j)=0.;
00296       }
00297       else
00298       {
00299         // b.elem(i,j)=sum/b.elem(j,j);
00300         dfsum+=dfb.elem(i,j)/b.elem(j,j);
00301         dfb.elem(j,j)-=dfb.elem(i,j)*b.elem(i,j)/b.elem(j,j);
00302         dfb.elem(i,j)=0.;
00303       }
00304 
00305       for (int k=min(i-1,j-1);k>=lb;k--)
00306       {
00307         // sum-=b.elem(i,k)*b.elem(k,j);
00308         dfb.elem(i,k)-=dfsum*b.elem(k,j);
00309         dfb.elem(k,j)-=dfsum*b.elem(i,k);
00310       }
00311       // sum=value(a(indx.elem(i),j);
00312       save_dmatrix_derivatives(a_pos,dfsum,indx.elem(i),j); // like this
00313       dfsum=0.;
00314     }
00315   }
00316 }
00317 
00318 #undef TINY