dmat38.cpp (revision 795)
25 25

26 26
#define TINY 1.0e-20;

27 27

28
dmatrix solve(const dmatrix& aa,const dmatrix& tz,

28
#include "ludcmp.hpp"

29

30
cltudecomp ludecomp_pivot(const dmatrix & M);

31

32
/** Solve a linear system using LU decomposition.

33
    \param aa A constant matrix. \f$A\f$.

34
    \param z A matrix containing the RHS, \f$B\f$ of the linear equation

35
    \f$A\cdot X = B\f$, to be solved.

36
    \return A matrix containing solution \f$X\f$.

37
*/

38
dmatrix solve(const dmatrix & aa, const dmatrix & zz)

39
{

40
   int n = aa.colsize();

41
   int lb = aa.colmin();

42
   int ub = aa.colmax();

43
   if (lb != aa.rowmin() || ub != aa.colmax())

44
   {

45
      cerr << "Error matrix not square in"

46
           << "solve(const dmatrix & aa, const dmatrix & z)" << endl;

47
      ad_exit(1);

48
   }

49
   dmatrix bb(lb, ub, lb, ub);

50
   bb = aa;

51
   cltudecomp dcmp;

52
   dcmp = ludecomp_pivot(bb);

53
   ivector index2 = dcmp.get_index2();

54

55
   //check if invertable

56
   int i;

57
   double det = 1.0;

58
   for (i = lb; i <= ub; i++)

59
   {

60
      det *= dcmp(i, i);

61
   }

62
   if (det == 0.0)

63
   {

64
      cerr <<

65
	 "Error in matrix inverse -- matrix singular in solve(dmatrix)\n";

66
      ad_exit(1);

67
   }

68

69
   //for each column of X and B solve A*x_i = b_i

70
   dmatrix xx(lb,ub,lb,ub);

71
   for(int k=lb;k<=ub;k++)

72
   {

73
     dvector z = column(zz,k);

74

75
     //Solve L*y=b with forward-substitution (before solving Ux=y)

76
     dvector y(lb, ub);

77
     for (i = lb; i <= ub; i++)

78
     {

79
        double tmp = 0.0;

80
        for (int j = lb; j < i; j++)

81
        {

82
	   tmp += dcmp(i, j) * y(j);

83
        }

84
        y(i) = z(index2(i)) - tmp;

85
     }

86

87
     //Now solve U*x=y with back substitution

88
     dvector x(lb, ub);

89
     for (i = ub; i >= lb; i--)

90
     {

91
        double tmp = 0.0;

92
        for (int j = ub; j > i; j--)

93
        {

94
	   tmp += dcmp(i, j) * x(j);

95
        }

96
        x(i) = (y(i) - tmp) / dcmp(i, i);

97
     }

98

99
     xx.colfill(k,x);

100
   }

101

102
   return xx;

103
}

104

105

106
/*dmatrix solve(const dmatrix& aa,const dmatrix& tz,

29 107
  double ln_unsigned_det,double& sign);

30 108

31 109
dmatrix solve(const dmatrix& aa,const dmatrix& tz)

......
195 273


196 274
  }

197 275
  return trans(x);

198
}

276
}*/

199 277

200 278
double ln_det_choleski(

201 279
  const banded_symmetric_dmatrix& MM, const int& _ierr)


Also available in: Unified diff