ADMB Documentation  11.2.2820
 All Classes Files Functions Variables Typedefs Friends Defines
f1b2sol2.cpp
Go to the documentation of this file.
00001 
00007 #include <df1b2fun.h>
00008 
00009 #ifdef __TURBOC__
00010   #pragma hdrstop
00011   #include <iostream.h>
00012 #endif
00013 
00014 #if defined (__WAT32__)
00015   #include <iostream.h>
00016   #include <strstrea.h>
00017 #endif
00018 
00019 #ifdef __ZTC__
00020   #include <iostream.hpp>
00021 #endif
00022 
00023 #define TINY 1.0e-20;
00024 df1b2vector solve(const df1b2matrix& aa,const df1b2vector& z,
00025   const df1b2variable & ln_unsigned_det,double& sign);
00026 
00027 df1b2vector solve(const df1b2matrix& aa,const df1b2vector& z,
00028   const df1b2variable& ld, df1b2variable& sign)
00029 {
00030   double sgn = 0;
00031   return solve(aa,z,ld,sgn);
00032 }
00033 
00034 df1b2vector csolve(const df1b2matrix& aa,const df1b2vector& z)
00035 {
00036   double ln_unsigned_det = 0;
00037   double sign = 0;
00038   df1b2vector sol=solve(aa,z,ln_unsigned_det,sign);
00039   return sol;
00040 }
00041 
00042 df1b2vector solve(const df1b2matrix& aa,const df1b2vector& z)
00043 {
00044   df1b2variable ln_unsigned_det;
00045   double sign;
00046   df1b2vector sol=solve(aa,z,ln_unsigned_det,sign);
00047   return sol;
00048 }
00049 
00055 df1b2vector solve(const df1b2matrix& aa,const df1b2vector& _z,
00056   const df1b2variable & _ln_unsigned_det,double& sign)
00057 {
00058   ADUNCONST(df1b2variable,ln_unsigned_det)
00059   ADUNCONST(df1b2vector,z)
00060   int k,n;
00061   n=aa.colsize();
00062   int lb=aa.colmin();
00063   int ub=aa.colmax();
00064   if (lb!=aa.rowmin()||ub!=aa.colmax())
00065   {
00066     cerr << "Error matrix not square in solve()"<<endl;
00067     ad_exit(1);
00068   }
00069   df1b2matrix bb(lb,ub,lb,ub);
00070   bb=aa;
00071   ivector indx(lb,ub);
00072   int One=1;
00073   indx.fill_seqadd(lb,One);
00074   double  d;
00075   df1b2variable big,dum,sum,temp;
00076   df1b2vector vv(lb,ub);
00077 
00078   d=1.0;
00079   for (int i=lb;i<=ub;i++)
00080   {
00081     big=0.0;
00082     for (int j=lb;j<=ub;j++)
00083     {
00084       temp=fabs(bb(i,j));
00085       if (value(temp) > value(big))
00086       {
00087         big=temp;
00088       }
00089     }
00090     if (value(big) == 0.0)
00091     {
00092       cerr <<
00093         "Error in matrix inverse -- matrix singular in inv(df1b2matrix)\n";
00094     }
00095     vv[i]=1.0/big;
00096   }
00097 
00098   for (int j=lb;j<=ub;j++)
00099   {
00100     for (int i=lb;i<j;i++)
00101     {
00102       sum=bb(i,j);
00103       for (k=lb;k<i;k++)
00104       {
00105         sum -= bb(i,k)*bb(k,j);
00106       }
00107       //a[i][j]=sum;
00108       bb(i,j)=sum;
00109     }
00110     int imax = j;
00111     big=0.0;
00112     for (int i=j;i<=ub;i++)
00113     {
00114       sum=bb(i,j);
00115       for (k=lb;k<j;k++)
00116       {
00117         sum -= bb(i,k)*bb(k,j);
00118       }
00119       bb(i,j)=sum;
00120       dum=vv[i]*fabs(sum);
00121       if ( value(dum) >= value(big))
00122       {
00123         big=dum;
00124         imax=i;
00125       }
00126     }
00127     if (j != imax)
00128     {
00129       for (k=lb;k<=ub;k++)
00130       {
00131         dum=bb(imax,k);
00132         bb(imax,k)=bb(j,k);
00133         bb(j,k)=dum;
00134       }
00135       d = -1.*d;
00136       vv[imax]=vv[j];
00137 
00138       //if (j<ub)
00139       {
00140         int itemp=indx(imax);
00141         indx(imax)=indx(j);
00142         indx(j)=itemp;
00143       }
00144       //cout << "indx= " <<indx<<endl;
00145     }
00146 
00147     if (value(bb(j,j)) == 0.0)
00148     {
00149       bb(j,j)=TINY;
00150     }
00151 
00152     if (j != n)
00153     {
00154       dum=1.0/bb(j,j);
00155       for (int i=j+1;i<=ub;i++)
00156       {
00157         bb(i,j) = bb(i,j) * dum;
00158       }
00159     }
00160   }
00161 
00162   // get the determinant
00163   sign=d;
00164   df1b2vector part_prod(lb,ub);
00165   part_prod(lb)=log(fabs(bb(lb,lb)));
00166   if (value(bb(lb,lb))<0) sign=-sign;
00167   for (int j=lb+1;j<=ub;j++)
00168   {
00169     if (value(bb(j,j))<0) sign=-sign;
00170     part_prod(j)=part_prod(j-1)+log(fabs(bb(j,j)));
00171   }
00172   ln_unsigned_det=part_prod(ub);
00173 
00174   df1b2vector x(lb,ub);
00175   df1b2vector y(lb,ub);
00176   //int lb=rowmin;
00177   //int ub=rowmax;
00178   df1b2matrix& b=bb;
00179   ivector indxinv(lb,ub);
00180   for (int i=lb;i<=ub;i++)
00181   {
00182     indxinv(indx(i))=i;
00183   }
00184 
00185   for (int i=lb;i<=ub;i++)
00186   {
00187     y(indxinv(i))=z(i);
00188   }
00189 
00190   for (int i=lb;i<=ub;i++)
00191   {
00192     sum=y(i);
00193     for (int j=lb;j<=i-1;j++)
00194     {
00195       sum-=b(i,j)*y(j);
00196     }
00197     y(i)=sum;
00198   }
00199   for (int i=ub;i>=lb;i--)
00200   {
00201     sum=y(i);
00202     for (int j=i+1;j<=ub;j++)
00203     {
00204       sum-=b(i,j)*x(j);
00205     }
00206     x(i)=sum/b(i,i);
00207   }
00208 
00209   return x;
00210 }
00211 #undef TINY