ADMB Documentation  11.1.2493
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_ma4.cpp
Go to the documentation of this file.
00001 
00007 #include "fvar.hpp"
00008 #include <math.h>
00009 
00010 static double eps0=1.e-50;
00011 
00012 void lubksb(dvar_matrix a, const ivector&  indx,dvar_vector b);
00013 void ludcmp(const dvar_matrix& a, const ivector& indx, const prevariable& d);
00014 
00023 void ludcmp(const dvar_matrix& _a, const ivector& _indx, const prevariable& _d)
00024 {
00025   ADUNCONST(dvar_matrix,a)
00026   ADUNCONST(prevariable,d)
00027   ivector& indx= (ivector&) _indx;
00028   int i,j,k,n;
00029 
00030   n=a.colsize();
00031   int lb=a.colmin();
00032   int ub=a.colmax();
00033 
00034   dvariable big,dum,sum,temp;
00035 
00036   dvar_vector vv(lb,ub);
00037 
00038 
00039   d=1.0;
00040 
00041   for (i=lb;i<=ub;i++)
00042   {
00043     big=0.0;
00044     for (j=lb;j<=ub;j++)
00045     {
00046       temp=fabs(a(i,j));
00047       if (temp > big)
00048       {
00049         big=temp;
00050       }
00051     }
00052     if (big == 0.0)
00053     {
00054       cerr << "Error in matrix inverse -- matrix singular in inv(dmatrix)\n";
00055     }
00056     vv(i)=1.0/big;
00057   }
00058 
00059   for (j=lb;j<=ub;j++)
00060   {
00061     for (i=lb;i<j;i++)
00062     {
00063       sum=a(i,j);
00064       for (k=lb;k<i;k++)
00065       {
00066         sum = sum - a(i,k)*a(k,j);
00067       }
00068       a(i,j)=sum;
00069     }
00070     int imax = j;
00071     big=0.0;
00072     for (i=j;i<=ub;i++)
00073     {
00074       sum=a(i,j);
00075       for (k=lb;k<j;k++)
00076       {
00077         sum = sum - a(i,k)*a(k,j);
00078       }
00079       a(i,j)=sum;
00080       dum=vv(i)*fabs(sum);
00081       if ( dum >= big)
00082       {
00083         big=dum;
00084         imax=i;
00085       }
00086     }
00087     if (j != imax)
00088     {
00089       for (k=lb;k<=ub;k++)
00090       {
00091         dum=a(imax,k);
00092         a(imax,k)=a(j,k);
00093         a(j,k)=dum;
00094       }
00095       d = -1.*d;
00096       vv(imax)=vv(j);
00097     }
00098     indx(j)=imax;
00099 
00100     if (a(j,j) == 0.0)
00101     {
00102       a(j,j)=eps0;
00103     }
00104 
00105     if (j != n)
00106     {
00107       dum=1.0/(a(j,j));
00108       for (i=j+1;i<=ub;i++)
00109       {
00110         a(i,j) = a(i,j) * dum;
00111       }
00112     }
00113   }
00114 }
00115 
00125 void lubksb(dvar_matrix a, const ivector& indx,dvar_vector b)
00126 {
00127   int i,ii=0,ip,j,iiflag=0;
00128   dvariable sum;
00129   int lb=a.colmin();
00130   int ub=a.colmax();
00131   for (i=lb;i<=ub;i++)
00132   {
00133     ip=indx(i);
00134     sum=b(ip);
00135     b(ip)=b(i);
00136     if (iiflag)
00137     {
00138       for (j=ii;j<=i-1;j++)
00139       {
00140         sum -= a.elem(i,j)*b.elem(j);
00141       }
00142     }
00143     else if ( value(sum) )
00144     {
00145       ii=i;
00146       iiflag=1;
00147     }
00148     b(i)=sum;
00149   }
00150 
00151   for (i=ub;i>=lb;i--)
00152   {
00153     sum=b(i);
00154     for (j=i+1;j<=ub;j++)
00155     {                        // !!! remove to show bug
00156       sum -= a.elem(i,j)*b.elem(j);
00157     }                        // !!! remove to show bug
00158     b.elem(i)=sum/a.elem(i,i);
00159   }
00160 }
00161 
00162