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