Revision 682 branches/davefoct2012/src/linad99/dmat3.cpp
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<=i1;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.e60; 

632 
int i=0; 

633 
int imax=0; 

634 
int j=0; 
Also available in: Unified diff