Revision 682

branches/davef-oct2012/src/linad99/dmat3.cpp (revision 682)
55 55
  return(y);
56 56
}
57 57

  
58
dmatrix inv_with_lu(const dmatrix& a,const ivector & indx,double d)
59
{
60
  if (a.rowmin()!=a.colmin() || a.rowmax() != a.colmax())
61
  {
62
    cerr << " Error in dmatrix inv(_CONST dmatrix&) -- matrix not square \n";
63
  }
64
 
65

  
66
  dmatrix y(a.rowmin(),a.rowmax(),a.rowmin(),a.rowmax());
67
  dvector col(a.rowmin(),a.rowmax());
68
  int i;
69

  
70
  for (int j=a.rowmin(); j<=a.rowmax(); j++)
71
  {
72
    for (i=a.rowmin(); i<=a.rowmax(); i++)
73
    {
74
      col[i]=0;
75
    }
76
    col[j]=1;
77

  
78
    lubksb(a,indx,col);
79
  
80
    for (i=a.rowmin(); i<=a.rowmax(); i++)
81
    {
82
      y[i][j]=col[i];
83
    }
84
  }
85
  return(y);
86
}
87

  
58 88
dmatrix inv(_CONST dmatrix& m1,const double& _ln_det, const int& _sgn)
59 89
{
60 90
  double d;
......
367 397
  }
368 398
}
369 399

  
400
dvector Lubksb(const dmatrix& a,const ivector&  indx,const dvector&  bb)
401
{
402
  int i,ii=0,ip,j,iiflag=0;
403
  dvector b(bb.indexmin(),bb.indexmax());
404
  b=bb;
405
  double sum;
406
  int lb=a.colmin();
407
  int ub=a.colmax();
408
  for (i=lb;i<=ub;i++)
409
  {
410
    ip=indx[i];
411
    sum=b[ip];
412
    b[ip]=b[i];
413
    if (iiflag)
414
    {
415
      for (j=ii;j<=i-1;j++)
416
      {
417
        sum -= a[i][j]*b[j];
418
      }
419
    }
420
    else if ( sum )
421
    {
422
      ii=i;
423
      iiflag=1;
424
    }
425
    b[i]=sum;
426
  }
427
 
428
  for (i=ub;i>=lb;i--) 
429
  {
430
    sum=b[i];
431
    for (j=i+1;j<=ub;j++) 
432
    {                        // !!! remove to show bug
433
      sum -= a[i][j]*b[j];
434
    }                        // !!! remove to show bug
435
    b[i]=sum/a[i][i];
436
  }
437
  return b;
438
}
439

  
370 440
double det(_CONST dmatrix& m1)
371 441
{
372 442
  double d;
......
553 623
  }
554 624
}
555 625
#undef TINY
626

  
627

  
628
dmatrix ludcmp(const dmatrix& _a,const ivector& _indx,
629
  const double& _det,const double& _sgn, const double& _d)
630
{
631
  const double TINY=1.e-60;
632
  int i=0;
633
  int imax=0;
634
  int j=0;
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff