systsolv.cpp 9.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439
  1. /*************LU DECOMPOSITION AND ITS APPLICATIONS***************/
  2. /*****Source - "NUMERICAL RECIPES IN C", www.nr.com*****/
  3. //#include "stdafx.h"
  4. #include <iostream.h>
  5. #include <math.h>
  6. #include <stdlib.h>
  7. #include "complex.h"
  8. #include "vector.h"
  9. #include "matrix.h"
  10. #include "cmplxvec.h"
  11. #include "cmplxmat.h"
  12. #include "systsolv.h"
  13. #define TINY 1.0e-30;
  14. /*********ludcmp - LU decomposition of a complex matrix***********/
  15. void ludcmp( CmplxMatrix& a, Vector& indx, double& d )
  16. {
  17. int const n = a.dim_i();
  18. int i, imax, j, k;
  19. double big, temp, dum;
  20. Vector vv(n), ind(n);
  21. Complex sum, dum1;
  22. d = 1.0;
  23. for( i=0; i < n; i++ ) {
  24. big = 0.0;
  25. for( j=0; j < n; j++ )
  26. if ( (temp = cabs(a(i,j))) > big ) big = temp;
  27. if( big == 0.0 ) {
  28. cerr << "Singular matrix in routine ludcmp !!!" << endl;
  29. exit(1);
  30. }
  31. vv[i] = 1.0/big;
  32. }
  33. for( j=0; j < n; j++ ) {
  34. for( i=0; i < j; i++ ) {
  35. sum = a(i,j);
  36. for( k=0; k < i; k++ ) sum -= a(i,k)*a(k,j);
  37. a(i,j) = sum;
  38. }
  39. big = 0.0;
  40. for( i=j; i < n; i++ ) {
  41. sum = a(i,j);
  42. for( k=0; k < j; k++ ) sum -= a(i,k)*a(k,j);
  43. a(i,j) = sum;
  44. if( (dum = vv[i]*cabs(sum)) >= big ) {
  45. big = dum;
  46. imax = i;
  47. }
  48. }
  49. if( j != imax ) {
  50. for( k=0; k < n; k++ ) {
  51. dum1 = a(imax,k);
  52. a(imax,k) = a(j,k);
  53. a(j,k) = dum1;
  54. }
  55. d = -d;
  56. vv[imax] = vv[j];
  57. }
  58. ind[j] = imax;
  59. if( a(j,j) == 0.0 ) a(j,j) = TINY;
  60. if( j != n-1 ) {
  61. dum1 = 1.0/a(j,j);
  62. for( i=j+1; i < n; i++ ) a(i,j) *= dum1; }
  63. }
  64. indx = ind;
  65. }
  66. /*******lubksb - solution of a system of linear equations*********/
  67. void lubksb( CmplxMatrix& a, CmplxVector& b, Vector& indx )
  68. {
  69. int const n = a.dim_i();
  70. int i, j, ip;
  71. int ii=-1;
  72. Complex sum(0.0);
  73. if (n != indx.dim()) {
  74. cerr << "Something is wrong in lubksb !!!" << endl;
  75. exit(1);
  76. }
  77. for(i=0; i<n; i++) {
  78. ip = int(indx[i]);
  79. sum = b[ip];
  80. b[ip] = b[i];
  81. if(ii!=-1)
  82. for(j=ii; j<=i-1; j++) sum -= a(i,j)*b[j];
  83. else if(sum!=0.0) ii = i;
  84. b[i] = sum;
  85. }
  86. for(i=n-1; i>=0; i--) {
  87. sum = b[i];
  88. for(j=i+1; j<n; j++) sum -= a(i,j)*b[j];
  89. b[i] = sum/a(i,i);
  90. }
  91. }
  92. /***********ludetr - determinant of a complex matrix**************/
  93. Complex ludetr( CmplxMatrix& a, double& d )
  94. {
  95. int const n = a.dim_i();
  96. Complex det(1.0);
  97. det *= d;
  98. for(int i=0; i<n; i++) det *= a(i,i);
  99. return det;
  100. }
  101. /*******************luinvm - inverse of a matrix*****************/
  102. void luinvm( CmplxMatrix& a, CmplxMatrix& mat, Vector& indx )
  103. {
  104. int const n = a.dim_i();
  105. int i, j;
  106. CmplxVector col(n);
  107. CmplxMatrix res(n,n);
  108. if (n != a.dim_j() || n != indx.dim()) {
  109. cerr << "Something is wrong in luinvm !!!" << endl;
  110. exit(1);
  111. }
  112. for(j=0; j<n; j++) {
  113. for(i=0; i<n; i++) col[i] = 0.0;
  114. col[j] = 1.0;
  115. lubksb( a, col, indx );
  116. for(i=0; i<n; i++) res(i,j) = col[i];
  117. }
  118. mat = res;
  119. }
  120. /*****improve - improvement of a solution to linear equations*****/
  121. void improve( CmplxMatrix& a, CmplxMatrix& alud, Vector& indx,
  122. CmplxVector& b, CmplxVector& x )
  123. {
  124. int const n = a.dim_i();
  125. //int i, j;
  126. CmplxVector r(n);
  127. r = a*x-b;
  128. lubksb( alud, r, indx );
  129. x -= r;
  130. }
  131. /*********ludcmp - LU decomposition of a real matrix**************/
  132. void ludcmp_real( Matrix& a, Vector& indx, double& d )
  133. {
  134. int const n = a.dim_i();
  135. int i, imax, j, k;
  136. double big, temp, dum;
  137. Vector vv(n), ind(n);
  138. double sum, dum1;
  139. d = 1.0;
  140. for( i=0; i < n; i++ ) {
  141. big = 0.0;
  142. for( j=0; j < n; j++ )
  143. if ( (temp = fabs(a(i,j))) > big ) big = temp;
  144. if( big == 0.0 ) {
  145. cerr << "Singular matrix in routine ludcmp !!!" << endl;
  146. exit(1);
  147. }
  148. vv[i] = 1.0/big;
  149. }
  150. for( j=0; j < n; j++ ) {
  151. for( i=0; i < j; i++ ) {
  152. sum = a(i,j);
  153. for( k=0; k < i; k++ ) sum -= a(i,k)*a(k,j);
  154. a(i,j) = sum;
  155. }
  156. big = 0.0;
  157. for( i=j; i < n; i++ ) {
  158. sum = a(i,j);
  159. for( k=0; k < j; k++ ) sum -= a(i,k)*a(k,j);
  160. a(i,j) = sum;
  161. if( (dum = vv[i]*fabs(sum)) >= big ) {
  162. big = dum;
  163. imax = i;
  164. }
  165. }
  166. if( j != imax ) {
  167. for( k=0; k < n; k++ ) {
  168. dum1 = a(imax,k);
  169. a(imax,k) = a(j,k);
  170. a(j,k) = dum1;
  171. }
  172. d = -d;
  173. vv[imax] = vv[j];
  174. }
  175. ind[j] = imax;
  176. if( a(j,j) == 0.0 ) a(j,j) = TINY;
  177. if( j != n-1 ) {
  178. dum1 = 1.0/a(j,j);
  179. for( i=j+1; i < n; i++ ) a(i,j) *= dum1; }
  180. }
  181. indx = ind;
  182. }
  183. /*******lubksb - solution of a system of linear equations*********/
  184. void lubksb_real( Matrix& a, Vector& b, Vector& indx )
  185. {
  186. int const n = a.dim_i();
  187. int i, j, ip;
  188. int ii=-1;
  189. double sum=0.0;
  190. if (n != indx.dim()) {
  191. cerr << "Something is wrong in lubksb !!!" << endl;
  192. exit(1);
  193. }
  194. for(i=0; i<n; i++) {
  195. ip = int(indx[i]);
  196. sum = b[ip];
  197. b[ip] = b[i];
  198. if(ii!=-1)
  199. for(j=ii; j<=i-1; j++) sum -= a(i,j)*b[j];
  200. else if(sum!=0.0) ii = i;
  201. b[i] = sum;
  202. }
  203. for(i=n-1; i>=0; i--) {
  204. sum = b[i];
  205. for(j=i+1; j<n; j++) sum -= a(i,j)*b[j];
  206. b[i] = sum/a(i,i);
  207. }
  208. }
  209. /***********ludetr - determinant of a complex matrix**************/
  210. double ludetr_real( Matrix& a, double& d )
  211. {
  212. int const n = a.dim_i();
  213. double det = 1.0;
  214. det *= d;
  215. for(int i=0; i<n; i++) det *= a(i,i);
  216. return det;
  217. }
  218. /*******************luinvm - inverse of a matrix*****************/
  219. void luinvm_real( Matrix& a, Matrix& mat, Vector& indx )
  220. {
  221. int const n = a.dim_i();
  222. int i, j;
  223. Vector col(n);
  224. Matrix res(n,n);
  225. if (n != a.dim_j() || n != indx.dim()) {
  226. cerr << "Something is wrong in luinvm !!!" << endl;
  227. exit(1);
  228. }
  229. for(j=0; j<n; j++) {
  230. for(i=0; i<n; i++) col[i] = 0.0;
  231. col[j] = 1.0;
  232. lubksb_real( a, col, indx );
  233. for(i=0; i<n; i++) res(i,j) = col[i];
  234. }
  235. mat = res;
  236. }
  237. /*****improve - improvement of a solution to linear equations*****/
  238. void improve_real( Matrix& a, Matrix& alud, Vector& indx,
  239. Vector& b, Vector& x )
  240. {
  241. int const n = a.dim_i();
  242. //int i, j;
  243. Vector r(n);
  244. r = a*x-b;
  245. lubksb_real( alud, r, indx );
  246. x -= r;
  247. }
  248. void get_l_u_real( Matrix& l, Matrix& u, const Matrix& a ) {
  249. const int n = a.dim_i();
  250. int i, j;
  251. double dlu;
  252. Vector Indx;
  253. Matrix res;
  254. if( n != a.dim_j() ) {
  255. cerr << " Something is wrong in getl_l_u_real !!! "
  256. << endl;
  257. exit(1);
  258. }
  259. res = a;
  260. ludcmp_real( res, Indx, dlu );
  261. l.resize(n,n);
  262. u.resize(n,n);
  263. for( i = 0; i < n; i++ ) {
  264. for( j = i; j < n; j++ )
  265. u(i,j) = res(i,j);
  266. // for( i = 0; i < n; i++ )
  267. l(i,i) = 1.0;
  268. // for( i = 1; i < n; i++ )
  269. for( j = 0; j < i; j++ )
  270. l(i,j) = res(i,j);
  271. }
  272. }
  273. void get_l_real( Matrix& l, const Matrix& u, const Matrix& a ) {
  274. const int n = a.dim_i();
  275. int i, j, k;
  276. double tmp;
  277. if( n != a.dim_j() || n != u.dim_i() || n != u.dim_j() ) {
  278. cerr << " Something is wrong in get_l_real !!! "
  279. << endl;
  280. exit(1);
  281. }
  282. l.resize(n,n);
  283. for( i = 0; i < n; i++ ) {
  284. l(i,0) = a(i,0) / u(0,0);
  285. // for( i = 1; i < n; i++ ) {
  286. for( j = 1; j <= i; j++ ) {
  287. tmp = 0.0;
  288. for( k = 0; k <= j-1; k++ )
  289. tmp += l(i,k) * u(k,j);
  290. l(i,j) = (a(i,j) - tmp) / u(j,j);
  291. }
  292. }
  293. }
  294. void solve_u_known_real( Matrix& l, const Matrix& u, const Matrix& a ) {
  295. const int n = a.dim_i();
  296. int i, j, k;
  297. double tmp;
  298. if( n != a.dim_j() || n != u.dim_i() || n != u.dim_j() ) {
  299. cerr << " Something is wrong in solve_u_known_real !!! "
  300. << endl;
  301. exit(1);
  302. }
  303. l.resize(n,n);
  304. for( i = 0; i < n; i++ ) {
  305. l(i,0) = a(i,0) / u(0,0);
  306. // for( i = 0; i < n; i++ ) {
  307. for( j = 1; j < n; j++ ) {
  308. tmp = 0.0;
  309. for( k = 0; k <= j-1; k++ )
  310. tmp += l(i,k) * u(k,j);
  311. l(i,j) = (a(i,j) - tmp) / u(j,j);
  312. }
  313. }
  314. }
  315. void solve_l_known_real( Matrix& u, const Matrix& l, const Matrix& a ) {
  316. const int n = a.dim_i();
  317. int i, j, k;
  318. double tmp;
  319. if( n != a.dim_j() || n != l.dim_i() || n != l.dim_j() ) {
  320. cerr << " Something is wrong in solve_l_known_real !!! "
  321. << endl;
  322. exit(1);
  323. }
  324. u.resize(n,n);
  325. for( i = 0; i < n; i++ ) {
  326. u(0,i) = a(0,i) / l(0,0);
  327. // for( i = 1; i < n; i++ ) {
  328. if ( i == 0 )
  329. continue;
  330. for( j = 0; j < n; j++ ) {
  331. tmp = 0.0;
  332. for( k = 0; k <= i-1; k++ )
  333. tmp += l(i,k) * u(k,j);
  334. u(i,j) = (a(i,j) - tmp) / l(i,i);
  335. }
  336. }
  337. }