ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
dmat36.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dmat36.cpp 1112 2013-07-12 21:41:41Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00017 lower_triangular_dmatrix::lower_triangular_dmatrix(int min,int max):
00018   dmatrix(min,max)
00019 {
00020   for (int i=min;i<=max;i++)
00021   {
00022     (*this)(i).allocate(min,i);
00023   }
00024 }
00025 
00030 dvector solve(const lower_triangular_dmatrix& m,const dvector&v)
00031 {
00032   int imin=m.indexmin();
00033   int imax=m.indexmax();
00034   if (v.indexmin() != imin || v.indexmax() != imax)
00035   {
00036     cerr << " Incompatible vector and matrix sizes in solve " << endl;
00037     ad_exit(1);
00038   }
00039   dvector x(imin,imax);
00040   x(imin)=v(imin)/m(imin,imin);
00041   for (int i=2;i<=imax;i++)
00042   {
00043     int jmin=imin;
00044     double ssum=0.0;
00045     for (int j=jmin;j<=i-1;j++)
00046     {
00047       ssum+=m(i,j)*x(j);
00048     }
00049     x(i)=(v(i)-ssum)/m(i,i);
00050   }
00051   return x;
00052 }
00053 
00058 dmatrix symmetrize(const lower_triangular_dmatrix & T)
00059 {
00060   int min=T.indexmin();
00061   int max=T.indexmax();
00062   dmatrix tmp(min,max,min,max);
00063   int i,j;
00064   for (i=min;i<=max;i++)
00065   {
00066     tmp(i,i)=T(i,i);
00067     for (j=i+1;j<=max;j++)
00068     {
00069       tmp(i,j)=T(j,i);
00070       tmp(j,i)=T(j,i);
00071     }
00072   }
00073   return tmp;
00074 }
00075 
00080 dvector solve_trans(const lower_triangular_dmatrix& M,const dvector& y)
00081 {
00082   int mmin=M.indexmin();
00083   int mmax=M.indexmax();
00084 
00085   if (y.indexmin() !=mmin || y.indexmax() !=mmax)
00086   {
00087     cerr << "incompatible size in solve_trans" << endl;
00088     ad_exit(1);
00089   }
00090   dvector x(mmin,mmax);
00091   int i,j;
00092 
00093   for (i=mmax;i>=mmin;i--)
00094   {
00095     double sum=0.0;
00096     for (j=i+1;j<=mmax;j++)
00097     {
00098       sum+=M(j,i)*x(j);
00099     }
00100     x(i)=(y(i)-sum)/M(i,i);
00101   }
00102 
00103   return x;
00104 }
00105 
00110 dmatrix fillout(const lower_triangular_dmatrix& M)
00111 {
00112   int mmin=M.indexmin();
00113   int mmax=M.indexmax();
00114   int i,j;
00115   dmatrix tmp(mmin,mmax,mmin,mmax);
00116   for (i=mmin;i<=mmax;i++)
00117   {
00118     for (j=1;j<i;j++)
00119       tmp(j,i)=0;
00120     for (j=i;j<=mmax;j++)
00121       tmp(j,i)=M(j,i);
00122   }
00123   return tmp;
00124 }
00125 
00130 dmatrix fillout_trans(const lower_triangular_dmatrix& M)
00131 {
00132   int mmin=M.indexmin();
00133   int mmax=M.indexmax();
00134   int i,j;
00135   dmatrix tmp(mmin,mmax,mmin,mmax);
00136   for (i=mmin;i<=mmax;i++)
00137   {
00138     for (j=1;j<i;j++)
00139       tmp(i,j)=0;
00140     for (j=i;j<=mmax;j++)
00141       tmp(i,j)=M(j,i);
00142   }
00143   return tmp;
00144 }