ADMB Documentation  11.1.2541
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2mat.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2mat.cpp 1967 2014-04-30 23:36:04Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <df1b2fun.h>
00017 df1b2matrix choleski_decomp(const df1b2matrix& MM)
00018 {
00019   // kludge to deal with constantness
00020   df1b2matrix & M= * (df1b2matrix *) &MM;
00021   int rmin=M.indexmin();
00022   int cmin=M(rmin).indexmin();
00023   int rmax=M.indexmax();
00024   int cmax=M(rmin).indexmax();
00025 
00026   if (rmin != cmin )
00027   {
00028     cerr << "minimum row and column inidices must equal 1 in "
00029       "df1b2matrix choleski_decomp(const df1b2matrix& MM)"
00030          << endl;
00031     ad_exit(1);
00032   }
00033   if (rmax !=cmax)
00034   {
00035     cerr << "Error in df1b2matrix choleski_decomp(const df1b2matrix& MM)"
00036       " Matrix not square" << endl;
00037     ad_exit(1);
00038   }
00039 
00040   //int n=rmax;
00041   df1b2matrix L(rmin,rmax,rmin,rmax);
00042 #ifndef SAFE_INITIALIZE
00043     L.initialize();
00044 #endif
00045 
00046   int i,j,k;
00047   df1b2variable tmp;
00048     if (value(M(rmin,rmin))<=0)
00049     {
00050       cerr << "Error matrix not positive definite in choleski_decomp"
00051         <<endl;
00052       ad_exit(1);
00053     }
00054 
00055   L(rmin,rmin)=sqrt(M(rmin,rmin));
00056   for (i=rmin+1;i<=rmax;i++)
00057   {
00058     L(i,rmin)=M(i,rmin)/L(rmin,rmin);
00059   }
00060 
00061   for (i=rmin+1;i<=rmax;i++)
00062   {
00063     for (j=rmin+1;j<=i-1;j++)
00064     {
00065       tmp=M(i,j);
00066       for (k=rmin;k<=j-1;k++)
00067       {
00068         tmp-=L(i,k)*L(j,k);
00069       }
00070       L(i,j)=tmp/L(j,j);
00071     }
00072     tmp=M(i,i);
00073     for (k=rmin;k<=i-1;k++)
00074     {
00075       tmp-=L(i,k)*L(i,k);
00076     }
00077     if (value(tmp)<=0)
00078     {
00079       cerr << "Error matrix not positive definite in choleski_decomp"
00080         <<endl;
00081       ad_exit(1);
00082     }
00083     L(i,i)=sqrt(tmp);
00084   }
00085   return L;
00086 }
00087 
00088 #undef HOME_VERSION