Revision 795 branches/replacement/src/linad99/dmat38.cpp
dmat38.cpp (revision 795)  

25  25  
26  26 
#define TINY 1.0e20; 
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 forwardsubstitution (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