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