ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
dmat37.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dmat37.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 dvector solve(const banded_symmetric_dmatrix& _m,const dvector&_v,
00018   const int& _ierr)
00019 {
00020   int & ierr=(int&)_ierr;
00021   ADUNCONST(dvector,v)
00022   ADUNCONST(banded_symmetric_dmatrix,m)
00023   int mmin=m.indexmin();
00024   m.shift(1);
00025   v.shift(1);
00026   const banded_lower_triangular_dmatrix& C=choleski_decomp(m,ierr);
00027   dvector w=solve_trans(C,solve(C,v));
00028   m.shift(mmin);
00029   w.shift(mmin);
00030   v.shift(mmin);
00031   return w;
00032 }
00033 
00038 dvector solve(const banded_symmetric_dmatrix& m,const dvector&v)
00039 {
00040   int ierr=0;
00041   return solve(m,v,ierr);
00042 }
00043 
00048 dmatrix solve(const banded_symmetric_dmatrix& m,const dmatrix& n,
00049   const int& _ierr)
00050 {
00051   int& ierr=(int&) _ierr;
00052   ierr=0;
00053   dmatrix tmp=trans(n);
00054   const banded_lower_triangular_dmatrix& C=choleski_decomp(m,ierr);
00055   int mmin=tmp.indexmin();
00056   int mmax=tmp.indexmax();
00057   dmatrix w(mmin,mmax);
00058   for (int i=mmin;i<=mmax;i++)
00059   {
00060     w(i)=solve_trans(C,solve(C,tmp(i)));
00061   }
00062   return trans(w);
00063 }
00064 
00069 dmatrix solve(const banded_symmetric_dmatrix& m,const dmatrix& n)
00070 {
00071   int ierr=0;
00072   return solve(m,n,ierr);
00073 }
00074 
00079 dvector solve(const banded_lower_triangular_dmatrix& m,const dvector&v)
00080 {
00081   int bw=m.bandwidth();
00082   int imin=m.indexmin();
00083   int imax=m.indexmax();
00084   if (v.indexmin() != imin || v.indexmax() != imax)
00085   {
00086     cerr << " Incompatible vector and matrix sizes in solve " << endl;
00087     ad_exit(1);
00088   }
00089   dvector x(imin,imax);
00090   x(imin)=v(imin)/m(imin,imin);
00091   for (int i=imin;i<=imax;i++)
00092   {
00093     int jmin=admax(imin,i-bw+1);
00094     double ssum=0.0;
00095     for (int j=jmin;j<=i-1;j++)
00096     {
00097       ssum+=m(i,j)*x(j);
00098     }
00099     x(i)=(v(i)-ssum)/m(i,i);
00100   }
00101   return x;
00102 }
00103 
00108 dvector solve_trans(const banded_lower_triangular_dmatrix& M,const dvector& y)
00109 {
00110   int mmin=M.indexmin();
00111   int mmax=M.indexmax();
00112   int bw=M.bandwidth();
00113 
00114   if (y.indexmin() !=mmin || y.indexmax() !=mmax)
00115   {
00116     cerr << "incompatible size in solve_trans" << endl;
00117     ad_exit(1);
00118   }
00119   dvector x(mmin,mmax);
00120   int i,j;
00121 
00122   for (i=mmax;i>=mmin;i--)
00123   {
00124     double sum=0.0;
00125     int jmax=admin(mmax,i+bw-1);
00126     for (j=i+1;j<=jmax;j++)
00127     {
00128       sum+=M(j,i)*x(j);
00129     }
00130     x(i)=(y(i)-sum)/M(i,i);
00131   }
00132 
00133   return x;
00134 }
00135 
00140 dvector lower_triangular_solve_trans(const dmatrix& M,const dvector& y)
00141 {
00142   int mmin=M.indexmin();
00143   int mmax=M.indexmax();
00144 
00145   if (y.indexmin() !=mmin || y.indexmax() !=mmax)
00146   {
00147     cerr << "incompatible size in solve_trans" << endl;
00148     ad_exit(1);
00149   }
00150   dvector x(mmin,mmax);
00151   int i,j;
00152 
00153   for (i=mmax;i>=mmin;i--)
00154   {
00155     double sum=0.0;
00156     int jmax=mmax;
00157     for (j=i+1;j<=jmax;j++)
00158     {
00159       sum+=M(j,i)*x(j);
00160     }
00161     x(i)=(y(i)-sum)/M(i,i);
00162   }
00163   return x;
00164 }
00165 
00170 dvector lower_triangular_solve(const dmatrix& m,const dvector&v)
00171 {
00172   int imin=m.indexmin();
00173   int imax=m.indexmax();
00174   if (v.indexmin() != imin || v.indexmax() != imax)
00175   {
00176     cerr << " Incompatible vector and matrix sizes in solve " << endl;
00177     ad_exit(1);
00178   }
00179   dvector x(imin,imax);
00180   x(imin)=v(imin)/m(imin,imin);
00181   for (int i=imin;i<=imax;i++)
00182   {
00183     int jmin=imin;
00184     double ssum=0.0;
00185     for (int j=jmin;j<=i-1;j++)
00186     {
00187       ssum+=m(i,j)*x(j);
00188     }
00189     x(i)=(v(i)-ssum)/m(i,i);
00190   }
00191   return x;
00192 }
00193 
00198 dvector choleski_solve_error(dmatrix M,dvector& v,int& ierror)
00199 {
00200   dmatrix C=choleski_decomp_error(M,ierror);
00201 
00202   if (ierror==1)
00203   {
00204     return v;
00205   }
00206   else
00207   {
00208     return lower_triangular_solve_trans(C,lower_triangular_solve(C,v));
00209   }
00210 }
00211 
00216 dvector choleski_solve_neghess_error(dmatrix M,dvector& v,int& ierror)
00217 {
00218   dmatrix C=choleski_decomp_neghess_error(M,ierror);
00219 
00220   if (ierror==1)
00221   {
00222     return v;
00223   }
00224   else
00225   {
00226     return lower_triangular_solve_trans(C,lower_triangular_solve(C,v));
00227   }
00228 }