Revision 795 branches/replacement/src/linad99/dmat38.cpp

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