ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
dmat15.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dmat15.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 
00013 #ifdef __TURBOC__
00014 #pragma hdrstop
00015 #include <iostream.h>
00016 #endif
00017 
00018 #ifdef __ZTC__
00019 #include <iostream.hpp>
00020 #endif
00021 
00022 #ifdef __SUN__
00023 #include <iostream.h>
00024 #endif
00025 #ifdef __NDPX__
00026 #include <iostream.h>
00027 #endif
00028 
00033 dmatrix choleski_decomp(const dmatrix& MM)
00034 {
00035   // kludge to deal with constantness
00036   dmatrix & M= * (dmatrix *) &MM;
00037   if (M.colsize() != M.rowsize())
00038   {
00039     cerr << "Error in chol_decomp. Matrix not square" << endl;
00040     ad_exit(1);
00041   }
00042   int rowsave=M.rowmin();
00043   int colsave=M.colmin();
00044   M.rowshift(1);
00045   M.colshift(1);
00046   int n=M.rowmax();
00047 
00048   dmatrix L(1,n,1,n);
00049 #ifndef SAFE_INITIALIZE
00050     L.initialize();
00051 #endif
00052 
00053   int i,j,k;
00054   double tmp;
00055     if (M(1,1)<=0)
00056     {
00057       cerr << "Error matrix not positive definite in choleski_decomp"
00058         <<endl;
00059       ad_exit(1);
00060     }
00061   L(1,1)=sqrt(M(1,1));
00062   for (i=2;i<=n;i++)
00063   {
00064     L(i,1)=M(i,1)/L(1,1);
00065   }
00066 
00067   for (i=2;i<=n;i++)
00068   {
00069     for (j=2;j<=i-1;j++)
00070     {
00071       tmp=M(i,j);
00072       for (k=1;k<=j-1;k++)
00073       {
00074         tmp-=L(i,k)*L(j,k);
00075       }
00076       L(i,j)=tmp/L(j,j);
00077     }
00078     tmp=M(i,i);
00079     for (k=1;k<=i-1;k++)
00080     {
00081       tmp-=L(i,k)*L(i,k);
00082     }
00083     if (tmp<=0)
00084     {
00085       cerr << "Error matrix not positive definite in choleski_decomp"
00086         <<endl;
00087       ad_exit(1);
00088     }
00089     L(i,i)=sqrt(tmp);
00090   }
00091   L.rowshift(rowsave);
00092   L.colshift(colsave);
00093 
00094   return L;
00095 }
00096 
00101 static dmatrix onerror(dmatrix& L,int & ierror)
00102 {
00103   ierror=1;
00104   return L;
00105 }
00106 
00111 dmatrix choleski_decomp_error(const dmatrix& MM,int& ierror)
00112 {
00113   ierror=0;
00114   // kludge to deal with constantness
00115   dmatrix & M= * (dmatrix *) &MM;
00116   if (M.colsize() != M.rowsize())
00117   {
00118     cerr << "Error in chol_decomp. Matrix not square" << endl;
00119     ad_exit(1);
00120   }
00121   int rowsave=M.rowmin();
00122   int colsave=M.colmin();
00123   M.rowshift(1);
00124   M.colshift(1);
00125   int n=M.rowmax();
00126 
00127   dmatrix L(1,n,1,n);
00128 #ifndef SAFE_INITIALIZE
00129     L.initialize();
00130 #endif
00131 
00132   int i,j,k;
00133   double tmp;
00134     if (M(1,1)<=0)
00135     {
00136       return onerror(L,ierror);
00137     }
00138   L(1,1)=sqrt(M(1,1));
00139   for (i=2;i<=n;i++)
00140   {
00141     L(i,1)=M(i,1)/L(1,1);
00142   }
00143 
00144   for (i=2;i<=n;i++)
00145   {
00146     for (j=2;j<=i-1;j++)
00147     {
00148       tmp=M(i,j);
00149       for (k=1;k<=j-1;k++)
00150       {
00151         tmp-=L(i,k)*L(j,k);
00152       }
00153       L(i,j)=tmp/L(j,j);
00154     }
00155     tmp=M(i,i);
00156     for (k=1;k<=i-1;k++)
00157     {
00158       tmp-=L(i,k)*L(i,k);
00159     }
00160     if (tmp<=0)
00161     {
00162       return onerror(L,ierror);
00163     }
00164     L(i,i)=sqrt(tmp);
00165   }
00166   L.rowshift(rowsave);
00167   L.colshift(colsave);
00168 
00169   return L;
00170 }
00171 
00176 dmatrix choleski_decomp_neghess_error(const dmatrix& MM, int& ierror)
00177 {
00178   ierror=0;
00179   // kludge to deal with constantness
00180   dmatrix & M= * (dmatrix *) &MM;
00181   if (M.colsize() != M.rowsize())
00182   {
00183     cerr << "Error in chol_decomp. Matrix not square" << endl;
00184     ad_exit(1);
00185   }
00186   int rowsave=M.rowmin();
00187   int colsave=M.colmin();
00188   M.rowshift(1);
00189   M.colshift(1);
00190   int n=M.rowmax();
00191 
00192   dmatrix L(1,n,1,n);
00193 #ifndef SAFE_INITIALIZE
00194     L.initialize();
00195 #endif
00196 
00197   int i,j,k;
00198   double tmp;
00199     if (M(1,1)>=0)
00200     {
00201       return onerror(L,ierror);
00202     }
00203   L(1,1)=sqrt(-M(1,1));
00204   for (i=2;i<=n;i++)
00205   {
00206     L(i,1)=-M(i,1)/L(1,1);
00207   }
00208 
00209   for (i=2;i<=n;i++)
00210   {
00211     for (j=2;j<=i-1;j++)
00212     {
00213       tmp=-M(i,j);
00214       for (k=1;k<=j-1;k++)
00215       {
00216         tmp-=L(i,k)*L(j,k);
00217       }
00218       L(i,j)=tmp/L(j,j);
00219     }
00220     tmp=-M(i,i);
00221     for (k=1;k<=i-1;k++)
00222     {
00223       tmp-=L(i,k)*L(i,k);
00224     }
00225     if (tmp<=0)
00226     {
00227       return onerror(L,ierror);
00228     }
00229     L(i,i)=sqrt(tmp);
00230   }
00231   L.rowshift(rowsave);
00232   L.colshift(colsave);
00233 
00234   return L;
00235 }