ADMB Documentation  11.1.1916
 All Classes Files Functions Variables Typedefs Friends Defines
f1b2v12.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: f1b2v12.cpp 1107 2013-07-11 21:56:15Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 
00017 df1b2vector lower_triangular_solve_trans(const df1b2matrix& M,
00018   const df1b2vector& y)
00019 {
00020   int mmin=M.indexmin();
00021   int mmax=M.indexmax();
00022 
00023   if (y.indexmin() !=mmin || y.indexmax() !=mmax)
00024   {
00025     cerr << "incompatible size in solve_trans" << endl;
00026     ad_exit(1);
00027   }
00028   df1b2vector x(mmin,mmax);
00029   int i,j;
00030 
00031   for (i=mmax;i>=mmin;i--)
00032   {
00033     df1b2variable ssum=0.0;
00034     int jmax=mmax;
00035     for (j=i+1;j<=jmax;j++)
00036     {
00037       ssum+=M(j,i)*x(j);
00038     }
00039     x(i)=(y(i)-ssum)/M(i,i);
00040   }
00041   return x;
00042 }
00043 
00048 df1b2vector lower_triangular_solve(const df1b2matrix& m,const df1b2vector& v)
00049 {
00050   int imin=m.indexmin();
00051   int imax=m.indexmax();
00052   if (v.indexmin() != imin || v.indexmax() != imax)
00053   {
00054     cerr << " Incompatible vector and matrix sizes in solve " << endl;
00055     ad_exit(1);
00056   }
00057   df1b2vector x(imin,imax);
00058   x(imin)=v(imin)/m(imin,imin);
00059   for (int i=imin+1;i<=imax;i++)
00060   {
00061     int jmin=imin;
00062     df1b2variable ssum=0.0;
00063     for (int j=jmin;j<=i-1;j++)
00064     {
00065       ssum+=m(i,j)*x(j);
00066     }
00067     x(i)=(v(i)-ssum)/m(i,i);
00068   }
00069   return x;
00070 }
00071 
00076 df1b2variable lower_triangular_ln_det(const df1b2matrix& m)
00077 {
00078   int sgn;
00079   return lower_triangular_ln_det(m,sgn);
00080 }
00081 
00086 df1b2variable lower_triangular_ln_det(const df1b2matrix& m,int& sgn)
00087 {
00088   sgn=1;
00089   int imin=m.indexmin();
00090   int imax=m.indexmax();
00091   df1b2variable ssum=0.0;
00092   for (int i=imin;i<=imax;i++)
00093   {
00094     if (value(m(i,i))>0.0)
00095     {
00096       ssum+=log(m(i,i));
00097     }
00098     else
00099     {
00100       sgn=-sgn;
00101       ssum+=log(-m(i,i));
00102     }
00103   }
00104   return ssum;
00105 }