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