Revision 752 branches/threaded2/src/linad99/dmat3.cpp
dmat3.cpp (revision 752)  

4  4 
* Author: David Fournier 
5  5 
* Copyright (c) 20092012 ADMB Foundation 
6  6 
*/ 
7 
#include "fvar.hpp" 

8  7 
#include <math.h> 
8 
#include <cassert> 

9 
#include "fvar.hpp" 

9  10 
#ifndef __ZTC__ 
10  11 
//#include <iomanip.h> 
11  12 
#endif 
...  ...  
19  20 
void lubksb(dmatrix a, const ivector& indx,dvector b); 
20  21 
void ludcmp(const dmatrix& a, const ivector& indx, const double& d); 
21  22  
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 

23 
/** 

24 
Inverse of a constant matrix by LU decomposition. 

25 
\ingroup matop 

26 
\param m1 A dmatrix, \f$M\f$, for which the inverse is to be computed. 

27 
\return A dmatrix containing \f$M^{1}\f$. 

28 
\n\n The implementation of this algorithm was inspired by 

29 
"Numerical Recipes in C", 2nd edition, 

30 
Press, Teukolsky, Vetterling, Flannery, chapter 2 

29  31 
*/ 
30  32 
dmatrix inv(const dmatrix& m1) 
31  33 
{ 
32 
double d; 

34  
33  35 
if (m1.rowmin()!=m1.colmin()  m1.rowmax() != m1.colmax()) 
34  36 
{ 
35  37 
cerr << " Error in dmatrix inv(const dmatrix&)  matrix not square \n"; 
36  38 
} 
39 
if (true) return m1; 

37  40 

38  41 
dmatrix a(m1.rowmin(),m1.rowmax(),m1.rowmin(),m1.rowmax()); 
39  
40 
int i; 

41 
for (i=m1.rowmin(); i<=m1.rowmax(); i++) 

42 
for (int i = m1.rowmin(); i <= m1.rowmax(); i++) 

42  43 
{ 
43 
for (int j=m1.rowmin(); j<=m1.rowmax(); j++)


44 
for (int j = m1.rowmin(); j <= m1.rowmax(); j++)


44  45 
{ 
45  46 
a[i][j]=m1[i][j]; 
46  47 
} 
47  48 
} 
48  49 
ivector indx(m1.rowmin(),m1.rowmax()); 
49 
//int indx[30]; 

50  
50 


51 
double d = 0; 

51  52 
ludcmp(a,indx,d); 
52  53  
53  54 
dmatrix y(m1.rowmin(),m1.rowmax(),m1.rowmin(),m1.rowmax()); 
...  ...  
55  56  
56  57 
for (int j=m1.rowmin(); j<=m1.rowmax(); j++) 
57  58 
{ 
58 
for (i=m1.rowmin(); i<=m1.rowmax(); i++)


59 
for (int i = m1.rowmin(); i <= m1.rowmax(); i++)


59  60 
{ 
60 
col[i]=0;


61 
col[i] = 0;


61  62 
} 
62 
col[j]=1;


63 
col[j] = 1;


63  64  
64  65 
lubksb(a,indx,col); 
65  66 

66 
for (i=m1.rowmin(); i<=m1.rowmax(); i++)


67 
for (int i = m1.rowmin(); i <= m1.rowmax(); i++)


67  68 
{ 
68 
y[i][j]=col[i];


69 
y[i][j] = col[i];


69  70 
} 
70  71 
} 
71 
return(y);


72 
return y;


72  73 
} 
73  74  
74  75 
/** Inverse of a constant matrix by LU decomposition. 
...  ...  
183  184 
double big,dum,sum,temp; 
184  185  
185  186 
dvector vv(lb,ub); 
186  
187  
188  187 
d=1.0; 
189  
190  188 
for (i=lb;i<=ub;i++) 
Also available in: Unified diff