ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
dmat38.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 
00025 dmatrix solve(const dmatrix& aa,const dmatrix& tz,
00026   const double& ln_unsigned_det,double& sign);
00027 
00032 dmatrix solve(const dmatrix& aa, const dmatrix& tz)
00033 {
00034   double ln = 0;
00035   double sgn = 0;
00036   return solve(aa,tz,ln,sgn);
00037 }
00038 
00045 dmatrix solve(const dmatrix& aa,const dmatrix& tz,
00046   const double& _ln_unsigned_det,double& sign)
00047 {
00048   double& ln_unsigned_det = (double&)_ln_unsigned_det;
00049   int n = aa.colsize();
00050   int lb=aa.colmin();
00051   int ub=aa.colmax();
00052   if (lb!=aa.rowmin()||ub!=aa.colmax())
00053   {
00054     cerr << "Error matrix not square in solve()"<<endl;
00055     ad_exit(1);
00056   }
00057   dmatrix bb(lb,ub,lb,ub);
00058   bb=aa;
00059   ivector indx(lb,ub);
00060   int One=1;
00061   indx.fill_seqadd(lb,One);
00062   dvector vv(lb,ub);
00063 
00064   double d = 1.0;
00065   for (int i=lb;i<=ub;i++)
00066   {
00067     double big=0.0;
00068     for (int j=lb;j<=ub;j++)
00069     {
00070       double temp=fabs(bb(i,j));
00071       if (temp > big)
00072       {
00073         big=temp;
00074       }
00075     }
00076     if (big == 0.0)
00077     {
00078       cerr << "Error in matrix inverse -- matrix singular in inv(dmatrix)\n";
00079     }
00080     vv[i]=1.0/big;
00081   }
00082 
00083   for (int j=lb;j<=ub;j++)
00084   {
00085     for (int i=lb;i<j;i++)
00086     {
00087       double sum=bb(i,j);
00088       for (int k=lb;k<i;k++)
00089       {
00090         sum -= bb(i,k)*bb(k,j);
00091       }
00092       //a[i][j]=sum;
00093       bb(i,j)=sum;
00094     }
00095     int imax = j;
00096     double big=0.0;
00097     for (int i=j;i<=ub;i++)
00098     {
00099       double sum=bb(i,j);
00100       for (int k=lb;k<j;k++)
00101       {
00102         sum -= bb(i,k)*bb(k,j);
00103       }
00104       bb(i,j)=sum;
00105       double dum=vv[i]*fabs(sum);
00106       if (dum >= big)
00107       {
00108         big=dum;
00109         imax=i;
00110       }
00111     }
00112     if (j != imax)
00113     {
00114       for (int k=lb;k<=ub;k++)
00115       {
00116         double dum=bb(imax,k);
00117         bb(imax,k)=bb(j,k);
00118         bb(j,k)=dum;
00119       }
00120       d = -1.*d;
00121       vv[imax]=vv[j];
00122 
00123       //if (j<ub)
00124       {
00125         int itemp=indx(imax);
00126         indx(imax)=indx(j);
00127         indx(j)=itemp;
00128       }
00129       //cout << "indx= " <<indx<<endl;
00130     }
00131 
00132     if (bb(j,j) == 0.0)
00133     {
00134       bb(j,j)=TINY;
00135     }
00136 
00137     if (j != n)
00138     {
00139       double dum=1.0/bb(j,j);
00140       for (int i=j+1;i<=ub;i++)
00141       {
00142         bb(i,j) = bb(i,j) * dum;
00143       }
00144     }
00145   }
00146 
00147   // get the determinant
00148   sign=d;
00149   dvector part_prod(lb,ub);
00150   part_prod(lb)=log(fabs(bb(lb,lb)));
00151   if (bb(lb,lb)<0) sign=-sign;
00152   for (int j=lb+1;j<=ub;j++)
00153   {
00154     if (bb(j,j)<0) sign=-sign;
00155     part_prod(j)=part_prod(j-1)+log(fabs(bb(j,j)));
00156   }
00157   ln_unsigned_det=part_prod(ub);
00158 
00159   dmatrix z=trans(tz);
00160   int mmin=z.indexmin();
00161   int mmax=z.indexmax();
00162   dmatrix x(mmin,mmax,lb,ub);
00163   //dvector x(lb,ub);
00164 
00165   dvector y(lb,ub);
00166   //int lb=rowmin;
00167   //int ub=rowmax;
00168   dmatrix& b=bb;
00169   ivector indxinv(lb,ub);
00170   for (int i=lb;i<=ub;i++)
00171   {
00172     indxinv(indx(i))=i;
00173   }
00174   for (int kk=mmin;kk<=mmax;kk++)
00175   {
00176     for (int i=lb;i<=ub;i++)
00177     {
00178       y(indxinv(i))=z(kk)(i);
00179     }
00180 
00181     for (int i=lb;i<=ub;i++)
00182     {
00183       double sum=y(i);
00184       for (int j=lb;j<=i-1;j++)
00185       {
00186         sum-=b(i,j)*y(j);
00187       }
00188       y(i)=sum;
00189     }
00190     for (int i=ub;i>=lb;i--)
00191     {
00192       double sum=y(i);
00193       for (int j=i+1;j<=ub;j++)
00194       {
00195         sum-=b(i,j)*x(kk)(j);
00196       }
00197       x(kk)(i)=sum/b(i,i);
00198     }
00199   }
00200   return trans(x);
00201 }
00202 
00203 double ln_det_choleski(
00204   const banded_symmetric_dmatrix& MM, const int& _ierr)
00205 {
00206   banded_lower_triangular_dmatrix tmp=choleski_decomp(MM,_ierr);
00207 
00208   int mmin=tmp.indexmin();
00209   int mmax=tmp.indexmax();
00210   double ld=0.0;
00211   for (int i=mmin;i<=mmax;i++)
00212   {
00213     ld+=log(tmp(i,i));
00214   }
00215   return 2.0*ld;
00216 }
00217 
00218 double norm(const banded_symmetric_dmatrix& B)
00219 {
00220   return sqrt(norm2(B));
00221 }
00222 
00223 double norm2(const banded_symmetric_dmatrix& B)
00224 {
00225   double nm=0.0;
00226   for (int i=1;i<=B.bw-1;i++)
00227   {
00228     nm+=norm2(B.d(i));
00229   }
00230   nm*=2;
00231   nm+=norm2(B.d(0));
00232   return nm;
00233 }
00234 
00235 #undef TINY
00236