ADMB Documentation  Fournier-pthread.1088
 All Classes Namespaces Files Functions Variables Typedefs Friends Defines
ludcmp.hpp
Go to the documentation of this file.
00001 /*
00002  * $Id: ludcmp.hpp 494 2012-06-13 20:41:16Z johnoel $
00003  *
00004  * Copyright (c) 2009-2012 ADMB Foundation
00005  */
00011 #ifndef __LUDCMP_HPP__
00012 #define __LUDCMP_HPP__
00013 #include <admodel.h>
00014 
00015 const double eps = 1.e-25;
00016 
00021 class cltudecomp
00022 {
00023    dmatrix L;
00024    dmatrix U;
00025    ivector indx;
00026    ivector indx2;
00027    int sign;
00028  public:
00029    void initialize(void)
00030    {
00031       indx.initialize();
00032       indx2.fill_seqadd(indexmin(), 1);
00033       sign = 1;
00034       L.initialize();
00035       U.initialize();
00036       for (int i = L.indexmin(); i <= L.indexmax(); i++)
00037       {
00038    L(i, i) = 1.0;
00039       }
00040    }
00041 
00042    cltudecomp(dmatrix & alpha, dmatrix & gamma);
00043 
00044    cltudecomp(void)
00045    {
00046    }
00047 
00048    void allocate(int lb, int ub)
00049    {
00050       indx.allocate(lb, ub);
00051       indx2.allocate(lb, ub);
00052       ivector iv(lb + 1, ub);
00053       iv.fill_seqadd(lb, 1);
00054       L.allocate(lb + 1, ub, lb, iv);
00055       ivector iv1(lb, ub);
00056       iv1.fill_seqadd(lb, 1);
00057       U.allocate(lb, ub, lb, iv1);
00058       indx2.fill_seqadd(lb, 1);
00059    }
00060    cltudecomp(int lb, int ub):indx(lb, ub), indx2(lb, ub)
00061    {
00062       ivector iv(lb + 1, ub);
00063       iv.fill_seqadd(lb, 1);
00064       L.allocate(lb + 1, ub, lb, iv);
00065       ivector iv1(lb, ub);
00066       iv1.fill_seqadd(lb, 1);
00067       U.allocate(lb, ub, lb, iv1);
00068       indx2.fill_seqadd(lb, 1);
00069    }
00070 
00071    cltudecomp & assign_value(const dvar_matrix & M)
00072    {
00073       int mmin = indexmin();
00074       int mmax = indexmax();
00075 
00076       if (mmin != M.indexmin() || mmax != M.indexmax())
00077       {
00078    cerr << "Shape error in =" << endl;
00079    ad_exit(1);
00080       }
00081       for (int i = mmin; i <= mmax; i++)
00082    for (int j = mmin; j <= mmax; j++)
00083       elem(i, j) = value(M(i, j));
00084       return *this;
00085    }
00086 
00087    dmatrix & get_L()
00088    {
00089       return L;
00090    }
00091    int indexmin()
00092    {
00093       return U.indexmin();
00094    }
00095    int indexmax()
00096    {
00097       return U.indexmax();
00098    }
00099    dmatrix & get_U()
00100    {
00101       return U;
00102    }
00103    ivector & get_index()
00104    {
00105       return indx;
00106    }
00107    ivector & get_index2()
00108    {
00109       return indx2;
00110    }
00111    int &get_sign()
00112    {
00113       return sign;
00114    }
00115    double &elem(int i, int j)
00116    {
00117       if (i > j)
00118    return L(i, j);
00119       else
00120    return U(j, i);
00121    }
00122    // overload () (int,int) to look like Numerical Recipes
00123    double &operator() (int i, int j)
00124    {
00125       if (i > j)
00126    return L(i, j);
00127       else
00128    return U(j, i);
00129    }
00130 };
00131 
00136 class dvector_for_adjoint
00137 {
00138    dmatrix D;
00139    ivector count;
00140  public:
00141    int indexmin()
00142    {
00143       return D.indexmin();
00144    }
00145    int indexmax()
00146    {
00147       return D.indexmax();
00148    }
00149    int indexmin() const
00150    {
00151       return D.indexmin();
00152    }
00153    int indexmax() const
00154    {
00155       return D.indexmax();
00156    }
00157    double &operator () (int i)
00158    {
00159       return D(i, count(i));
00160    }
00161    const double& operator () (int i) const
00162    {
00163       return D(i, count(i));
00164    }
00165    dvector_for_adjoint operator () (int mmin, int mmax)
00166    {
00167       return dvector_for_adjoint(D.sub(mmin, mmax), count(mmin, mmax));
00168    }
00169  dvector_for_adjoint(const dmatrix & _D, const ivector & _count):
00170    D(_D), count(_count)
00171    {
00172    }
00173 };
00174 
00175 double operator *
00176    (const dvector_for_adjoint & v, const dvector_for_adjoint & w);
00177 
00178 dvector operator *(double x, const dvector_for_adjoint & w);
00179 
00184 class dmatrix_for_adjoint
00185 {
00186    d3_array D;
00187    imatrix count;
00188  public:
00189    void initialize(void)
00190    {
00191       D.initialize();
00192       count = 1.0;
00193    }
00194    int indexmin()
00195    {
00196       return D.indexmin();
00197    }
00198    int indexmax()
00199    {
00200       return D.indexmax();
00201    }
00202    dvector_for_adjoint operator () (int i)
00203    {
00204       return dvector_for_adjoint(D(i), count(i));
00205    }
00206 
00207    void increment(int i, int j)
00208    {
00209       count(i, j)++;
00210       if (count(i, j) > 3)
00211    cout << count(i, j) << endl;
00212    }
00213    void decrement(int i, int j)
00214    {
00215       count(i, j)--;
00216    }
00217 
00218    double &operator () (int i, int j)
00219    {
00220       return D(i, j, count(i, j));
00221    }
00222    void allocate(int lb, int ub, int l1, ivector & iv, int l2, int u2)
00223    {
00224       D.allocate(lb, ub, l1, iv, l2, u2);
00225       count.allocate(lb, ub, l1, iv);
00226       count = 1;
00227    }
00228 };
00229 
00234 class cltudecomp_for_adjoint
00235 {
00236    dmatrix_for_adjoint L;
00237    dmatrix_for_adjoint U;
00238    ivector indx;
00239    ivector indx2;
00240    cltudecomp dfclu;
00241    dvar_matrix_position *pMpos;
00242    int sign;
00243  public:
00244    void initialize(void)
00245    {
00246       pMpos = 0;
00247       indx.initialize();
00248       indx2.fill_seqadd(indexmin(), 1);
00249       sign = 1;
00250       L.initialize();
00251       U.initialize();
00252       for (int i = L.indexmin(); i <= L.indexmax(); i++)
00253       {
00254    L(i, i) = 1.0;
00255       }
00256    }
00257 
00258    int indexmin()
00259    {
00260       return U.indexmin();
00261    }
00262    int indexmax()
00263    {
00264       return U.indexmax();
00265    }
00266 
00267    void ludecomp_pivot_for_adjoint_1(void);
00268 
00269    void ludecomp_pivot_for_adjoint_2(void);
00270    cltudecomp_for_adjoint(void)
00271    {
00272       pMpos = 0;
00273    }
00274 
00275    ~cltudecomp_for_adjoint()
00276    {
00277       if (pMpos)
00278    delete pMpos;
00279    }
00280 
00281    void allocate(int lb, int ub, int n, int m)
00282    {
00283       indx.allocate(lb, ub);
00284       indx2.allocate(lb, ub);
00285       dfclu.allocate(lb, ub);
00286       ivector iv(lb + 1, ub);
00287       iv.fill_seqadd(lb, 1);
00288       L.allocate(lb + 1, ub, lb, iv, 1, n);
00289       ivector iv1(lb, ub);
00290       iv1.fill_seqadd(lb, 1);
00291       U.allocate(lb, ub, lb, iv1, 1, m);
00292    }
00293 
00294    cltudecomp_for_adjoint(int lb, int ub, int n, int m):indx(lb, ub),
00295       dfclu(lb, ub), pMpos(0)
00296    {
00297       ivector iv(lb + 1, ub);
00298       iv.fill_seqadd(lb, 1);
00299       L.allocate(lb + 1, ub, lb, iv, 1, n);
00300       ivector iv1(lb, ub);
00301       iv1.fill_seqadd(lb, 1);
00302       U.allocate(lb, ub, lb, iv1, 1, m);
00303    }
00304 
00305    cltudecomp & get_dfclu()
00306    {
00307       return dfclu;
00308    }
00309    dmatrix_for_adjoint & get_L()
00310    {
00311       return L;
00312    }
00313    dmatrix_for_adjoint & get_U()
00314    {
00315       return U;
00316    }
00317    ivector & get_index()
00318    {
00319       return indx;
00320    }
00321    int &get_sign()
00322    {
00323       return sign;
00324    }
00325    // overload () (int,int) to look like Numerical Recipes
00326    void decrement(int i, int j)
00327    {
00328       if (i > j)
00329    L.decrement(i, j);
00330       else
00331    U.decrement(j, i);
00332    }
00333    void increment(int i, int j)
00334    {
00335       if (i > j)
00336    L.increment(i, j);
00337       else
00338    U.increment(j, i);
00339    }
00340    double &operator () (int i, int j)
00341    {
00342       if (i > j)
00343    return L(i, j);
00344       else
00345    return U(j, i);
00346    }
00347    friend class double_for_assign;
00348 };
00349 
00350 
00351 dmatrix get_dmatrix(cltudecomp & clu);
00352 dmatrix get_dmatrix(cltudecomp_for_adjoint & clu);
00353 dmatrix get_lower_matrix(cltudecomp_for_adjoint & clu);
00354 dmatrix get_upper_matrix(cltudecomp_for_adjoint & clu);
00355 dvar_vector solve(const dvar_matrix & aa, const dvar_vector & z,
00356       prevariable & ln_unsigned_det,
00357       const prevariable & _sign);
00358 #endif        //#ifndef __LUDCMP_HPP__