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