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