Revision 692 branches/mergetrunkdavef/src/linad99/dmat3.cpp
dmat3.cpp (revision 692)  

1 
/** 

2 
* $Id: dmat3.cpp 789 20101005 01:01:09Z johnoel $ 

3 
* 

4 
* Author: David Fournier 

5 
* Copyright (c) 20092012 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.0e20; 
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 

Also available in: Unified diff