Revision 752 branches/threaded2/src/linad99/dmat3.cpp

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