ADMB Documentation  11.1.2397
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m53.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar_m53.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00017 dvar_vector lower_triangular_solve(const dvar_matrix& m,const dvar_vector&v)
00018 {
00019   int imin=m.indexmin();
00020   int imax=m.indexmax();
00021   if (v.indexmin() != imin || v.indexmax() != imax)
00022   {
00023     cerr << " Incompatible vector and matrix sizes in solve " << endl;
00024     ad_exit(1);
00025   }
00026   dvar_vector x(imin,imax);
00027   x(imin)=v(imin)/m(imin,imin);
00028   for (int i=2;i<=imax;i++)
00029   {
00030     int jmin=imin;
00031     dvariable ssum=0.0;
00032     for (int j=jmin;j<=i-1;j++)
00033     {
00034       ssum+=m(i,j)*x(j);
00035     }
00036     x(i)=(v(i)-ssum)/m(i,i);
00037   }
00038   return x;
00039 }
00040 
00045 dvariable lower_triangular_ln_det(const dvar_matrix& m)
00046 {
00047   int sgn;
00048   return lower_triangular_ln_det(m,sgn);
00049 }
00050 
00055 dvariable lower_triangular_ln_det(const dvar_matrix& m,int& sgn)
00056 {
00057   sgn=1;
00058   int imin=m.indexmin();
00059   int imax=m.indexmax();
00060   dvariable ssum=0.0;
00061   for (int i=imin;i<=imax;i++)
00062   {
00063     if (value(m(i,i))>0.0)
00064     {
00065       ssum+=log(m(i,i));
00066     }
00067     else
00068     {
00069       sgn=-sgn;
00070       ssum+=log(-m(i,i));
00071     }
00072   }
00073   return ssum;
00074 }
00075 
00076 
00077 
00078 /*
00079 dvar_vector lower_triangular_solve(const dvar_matrix& a, const dvar_vector & b)
00080 {
00081   int i,k;
00082   dvariable sum;
00083 
00084   int mmin=b.indexmin();
00085   int mmax=b.indexmax();
00086   dvar_vector x(mmin,mmax);
00087 
00088   for (i=mmin;i<=mmax;i++)
00089   {
00090     sum=b[i];
00091     for (k=i-1;k>=mmin;k--)
00092     {
00093       sum -= a[i][k]*x[k];
00094     }
00095     x[i]=sum/a[i][i];
00096   }
00097   return x;
00098   for (i=mmax;i>=mmin;i--)
00099   {
00100     for (sum=x[i],k=i+1;k<=mmax;k++)
00101     {
00102       sum -= a[k][i]*x[k];
00103     }
00104     x[i]=sum/a[i][i];
00105   }
00106 }
00107  */
00108