ADMB Documentation  11.1.2397
 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,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     int imax = j;
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