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