Revision 692 branches/merge-trunk-davef/src/linad99/dmat3.cpp

dmat3.cpp (revision 692)
1
/**
2
 * $Id: dmat3.cpp 789 2010-10-05 01:01:09Z johnoel $
3
 *
4
 * Author: David Fournier
5
 * Copyright (c) 2009-2012 ADMB Foundation
6
 */
1
#define HOME_VERSION
2
//COPYRIGHT (c) 1991 OTTER RESEARCH LTD
3

  
7 4
#include "fvar.hpp"
8 5
#include <math.h>
9 6
#ifndef __ZTC__
10 7
//#include <iomanip.h>
11 8
#endif
12
/**
13
\def TINY
14
A small number. Used to avoid divide by zero in the LU decomposition. Locally defined,
15
undefined, redefined and undefined in this file.
16
*/
17 9
#define TINY 1.0e-20;
18 10

  
19
void lubksb(dmatrix a, const ivector&  indx,dvector b);
20
void ludcmp(const dmatrix& a, const ivector& indx, const double& d);
11
void lubksb(dmatrix a,_CONST ivector&  indx,dvector b);
12
void ludcmp(BOR_CONST dmatrix& a,BOR_CONST ivector& indx,BOR_CONST double& d);
21 13

  
22
/** Inverse of a constant matrix by LU decomposition.
23
    \ingroup matop
24
    \param m1 A dmatrix, \f$M\f$, for which the inverse is to be computed.
25
    \return A dmatrix containing \f$M^{-1}\f$.
26
    \n\n The implementation of this algorithm was inspired by
27
    "Numerical Recipes in C", 2nd edition,
28
    Press, Teukolsky, Vetterling, Flannery, chapter 2
29
*/
30
dmatrix inv(const dmatrix& m1)
14
dmatrix inv(_CONST dmatrix& m1)
31 15
{
32 16
  double d;
33 17
  if (m1.rowmin()!=m1.colmin() || m1.rowmax() != m1.colmax())
34 18
  {
35
    cerr << " Error in dmatrix inv(const dmatrix&) -- matrix not square \n";
19
    cerr << " Error in dmatrix inv(_CONST dmatrix&) -- matrix not square \n";
36 20
  }
37 21
 
38 22
  dmatrix a(m1.rowmin(),m1.rowmax(),m1.rowmin(),m1.rowmax());
......
71 55
  return(y);
72 56
}
73 57

  
74
/** Inverse of a constant matrix by LU decomposition.
75
    \param m1 A dmatrix, \f$M\f$, for which the inverse is to be computed.
76
    \param _ln_det On return contains \f$|\log M|\f$
77
    \param _sign
78
    \return A dmatrix containing \f$M^{-1}\f$.
79
    \n\n The implementation of this algorithm was inspired by
80
    "Numerical Recipes in C", 2nd edition,
81
    Press, Teukolsky, Vetterling, Flannery, chapter 2
82
*/
83
dmatrix inv(const dmatrix& m1,const double& _ln_det, const int& _sgn)
58
dmatrix inv_with_lu(const dmatrix& a,const ivector & indx,double d)
84 59
{
60
  if (a.rowmin()!=a.colmin() || a.rowmax() != a.colmax())
61
  {
62
    cerr << " Error in dmatrix inv(_CONST dmatrix&) -- matrix not square \n";
63
  }
64
 
65

  
66
  dmatrix y(a.rowmin(),a.rowmax(),a.rowmin(),a.rowmax());
67
  dvector col(a.rowmin(),a.rowmax());
68
  int i;
69

  
70
  for (int j=a.rowmin(); j<=a.rowmax(); j++)
71
  {
72
    for (i=a.rowmin(); i<=a.rowmax(); i++)
73
    {
74
      col[i]=0;
75
    }
76
    col[j]=1;
77

  
78
    lubksb(a,indx,col);
79
  
80
    for (i=a.rowmin(); i<=a.rowmax(); i++)
81
    {
82
      y[i][j]=col[i];
83
    }
84
  }
85
  return(y);
86
}
87

  
88
dmatrix inv(_CONST dmatrix& m1,const double& _ln_det, const int& _sgn)
89
{
85 90
  double d;
86 91
  double& ln_det=(double&)(_ln_det);
87 92
  ln_det=0.0;
......
89 94
  
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff