dmat3.cpp (revision 752)
4 4
``` * Author: David Fournier
```
5 5
``` * Copyright (c) 2009-2012 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++)
```
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff