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

1  1 
// For AD Model Builder by Anders Nielsen <anders@nielsensweb.org> and Casper Berg <cbe@aqua.dtu.dk> 
2  2 
// Aug. 2010 
3  3 
#include <fvar.hpp> 
4 
#include "ludcmp.hpp" 

5 
dvar_matrix solve(const dvar_matrix& aa, const dvar_matrix& zz); 

6 
cltudecomp xludecomp_pivot(const dvar_matrix & M); 

7 
static void df_mm_solve(void); 

8  
4  9 
#define TINY 1.0e20; 
5  10  
6  11 
dmatrix fabs(const dmatrix & X){ 
...  ...  
18  23 
} 
19  24  
20  25  
21 
dvar_matrix solve(const dvar_matrix& aa, const dvar_matrix& tz, dvariable ln_unsigned_det, dvariable& sign); 

22  26  
27  
28 
/*dvar_matrix solve(const dvar_matrix& aa, const dvar_matrix& tz, dvariable ln_unsigned_det, dvariable& sign); 

29  
23  30 
dvar_matrix solve(const dvar_matrix& aa, const dvar_matrix& tz) 
24  31 
{ 
25  32 
dvariable ln; 
...  ...  
188  195 
} 
189  196 
RETURN_ARRAYS_DECREMENT(); 
190  197 
return trans(x); 
191 
} 

198 
}*/


192  199  
193  200 
/** 
194  201 
\ingroup PDF 
...  ...  
247  254 
return E; 
248  255 
} 
249  256  
257 
#undef TINY 

250  258  
251 
#undef TINY 

259  
260 
/** Solve a linear system using LU decomposition. 

261 
\param aa A variable matrix. \f$A\f$. 

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

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

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

265 
*/ 

266 
dvar_matrix solve(const dvar_matrix& aa, const dvar_matrix& zz) 

267 
{ 

268 
int n=aa.colsize(); 

269 
int lb=aa.colmin(); 

270 
int ub=aa.colmax(); 

271 
if (lb!=aa.rowmin()ub!=aa.colmax()) 

272 
{ 

273 
cerr << "Error matrix not square in " 

274 
<< "solve(const dvar_matrix& aa, const dvar_matrix& zz)"<<endl; 

275 
ad_exit(1); 

276 
} 

277 
if (zz.rowmin()!=aa.rowmin()zz.colmax()!=aa.colmax()) 

278 
{ 

279 
cerr << "Error matrices not compatible in " 

280 
<< "solve(const dvar_matrix& aa, const dvar_matrix& zz)"<<endl; 

281 
ad_exit(1); 

282 
} 

283  
284 
dvar_matrix xx(lb,ub,lb,ub); 

285  
286 
if(ub==lb) 

287 
{ 

288 
xx=zz*inv(aa); 

289 
return(xx); 

290 
} 

291 


292 


293  
294 
cltudecomp clu1=xludecomp_pivot(aa); 

295 
ivector index2=clu1.get_index2(); 

296 
dmatrix & gamma=clu1.get_U(); 

297 
dmatrix & alpha=clu1.get_L(); 

298  
299 
//check if invertable 

300 
int i; 

301 
double det=1.0; 

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

303 
{ 

304 
det*=clu1(i,i); 

305 
} 

306 
if (det==0.0) 

307 
{ 

308 
cerr << "Error in matrix inverse  matrix singular in solve(dmatrix)\n"; 

309 
ad_exit(1); 

310 
} 

311  
312  
313 
dmatrix yy(lb,ub,lb,ub); 

314 
dmatrix ttmp1(lb,ub,lb,ub); 

315 
dmatrix ttmp2(lb,ub,lb,ub); 
Also available in: Unified diff