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