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.0e-20;
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);
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff