bicgstab.cpp 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566
  1. /*****************************************************************/
  2. /*** ***/
  3. /*** Bi_CGSTAB iterative method ***/
  4. /*** ***/
  5. /*** A - system matrix, xi - initial guess and output vector, ***/
  6. /*** b - right hand side ***/
  7. /*** ***/
  8. /*****************************************************************/
  9. //#include "stdafx.h"
  10. #include <iostream.h>
  11. #include <math.h>
  12. #include "vector.h"
  13. #include "matrix.h"
  14. #include "complex.h"
  15. #include "cmplxvec.h"
  16. #include "cmplxmat.h"
  17. #include "ivectorl.h"
  18. #include "sparse.h"
  19. #include "bicgstab.h"
  20. int bi_cgstab_d( Matrix& A, Vector& xi, Vector& b, double EPS ) {
  21. const int ns = A.dim_i();
  22. int maxit = 2000;
  23. int iflag = 1, icount = 0;
  24. double delta0, deltai;
  25. double bet, roi;
  26. double roim1 = 1.0, al = 1.0, wi = 1.0;
  27. Vector ri, roc;
  28. Vector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
  29. ri = b - A * xi;
  30. roc = ri;
  31. delta0 = norm2( ri );
  32. while( iflag != 0 && icount < maxit ) {
  33. icount += 1;
  34. roi = roc * ri;
  35. bet = ( roi / roim1 ) * ( al / wi );
  36. pi = ri + ( pi - vi * wi ) * bet;
  37. vi = A * pi;
  38. al = roi / ( roc * vi );
  39. s = ri - vi * al;
  40. t = A * s;
  41. wi = ( t * s ) / ( t * t );
  42. xi += pi * al + s * wi;
  43. ri = s - t * wi;
  44. deltai = norm2( ri );
  45. // cerr << " " << icount << " "
  46. // << setprecision(10)
  47. // << deltai << " "
  48. // << delta0 << " "
  49. // << deltai/delta0 << endl;
  50. if( deltai < EPS * delta0 )
  51. iflag = 0;
  52. else
  53. roim1 = roi;
  54. }
  55. cout << endl
  56. << endl
  57. << " Number of iterations = "
  58. << icount
  59. << endl << endl;
  60. return( 0 );
  61. }
  62. int bi_cgstab_c( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS ) {
  63. const int ns = A.dim_i();
  64. int maxit = 2000;
  65. int iflag = 1, icount = 0;
  66. double delta0, deltai;
  67. Complex bet, roi;
  68. Complex roim1 = 1.0, al = 1.0, wi = 1.0;
  69. CmplxVector ri, roc;
  70. CmplxVector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
  71. ri = b - A * xi;
  72. roc = ri;
  73. delta0 = norm2( ri );
  74. while( iflag != 0 && icount < maxit ) {
  75. icount += 1;
  76. roi = inner( roc, ri );
  77. bet = ( roi / roim1 ) * ( al / wi );
  78. pi = ri + ( pi - vi * wi ) * bet;
  79. vi = A * pi;
  80. al = roi / inner( roc, vi );
  81. s = ri - vi * al;
  82. t = A * s;
  83. wi = inner( t, s ) / inner( t, t );
  84. xi += pi * al + s * wi;
  85. ri = s - t * wi;
  86. deltai = norm2( ri );
  87. if( deltai < EPS * delta0 )
  88. iflag = 0;
  89. else
  90. roim1 = roi;
  91. }
  92. cerr << endl
  93. << endl
  94. << " Number of iterations = "
  95. << icount
  96. << endl << endl;
  97. return( 0 );
  98. }
  99. int bi_cgstab_d_sparse( Matrix& A, Vector& xi, Vector& b, double EPS, double thresh ) {
  100. const int ns = A.dim_i();
  101. int maxit = 2000;
  102. int iflag = 1, icount = 0;
  103. double delta0, deltai;
  104. double bet, roi;
  105. double roim1 = 1.0, al = 1.0, wi = 1.0;
  106. Vector ri, roc;
  107. Vector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
  108. Vector sa;
  109. IVectorl ija;
  110. sprsin_d( A, thresh, sa, ija );
  111. ri = b - sprsax_d_v( sa, ija, xi );
  112. roc = ri;
  113. delta0 = norm2( ri );
  114. while( iflag != 0 && icount < maxit ) {
  115. icount += 1;
  116. roi = roc * ri;
  117. bet = ( roi / roim1 ) * ( al / wi );
  118. pi = ri + ( pi - vi * wi ) * bet;
  119. vi = sprsax_d_v( sa, ija, pi );
  120. al = roi / ( roc * vi );
  121. s = ri - vi * al;
  122. t = sprsax_d_v( sa, ija, s );
  123. wi = ( t * s ) / ( t * t );
  124. xi += pi * al + s * wi;
  125. ri = s - t * wi;
  126. deltai = norm2( ri );
  127. if( deltai < EPS * delta0 )
  128. iflag = 0;
  129. else
  130. roim1 = roi;
  131. }
  132. cerr << endl
  133. << endl
  134. << " Number of iterations = "
  135. << icount
  136. << endl << endl;
  137. return( 0 );
  138. }
  139. int bi_cgstab_c_sparse( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS, double thresh ) {
  140. const int ns = A.dim_i();
  141. int maxit = 2000;
  142. int iflag = 1, icount = 0;
  143. double delta0, deltai;
  144. Complex bet, roi;
  145. Complex roim1 = 1.0, al = 1.0, wi = 1.0;
  146. CmplxVector ri, roc;
  147. CmplxVector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
  148. CmplxVector sa;
  149. IVectorl ija;
  150. sprsin_c( A, thresh, sa, ija );
  151. ri = b - sprsax_c_v( sa, ija, xi );
  152. roc = ri;
  153. delta0 = norm2( ri );
  154. while( iflag != 0 && icount < maxit ) {
  155. icount += 1;
  156. roi = inner( roc, ri );
  157. bet = ( roi / roim1 ) * ( al / wi );
  158. pi = ri + ( pi - vi * wi ) * bet;
  159. vi = sprsax_c_v( sa, ija, pi );
  160. al = roi / inner( roc, vi );
  161. s = ri - vi * al;
  162. t = sprsax_c_v( sa, ija, s );
  163. wi = inner( t, s ) / inner( t, t );
  164. xi += pi * al + s * wi;
  165. ri = s - t * wi;
  166. deltai = norm2( ri );
  167. if( deltai < EPS * delta0 )
  168. iflag = 0;
  169. else
  170. roim1 = roi;
  171. }
  172. cerr << endl
  173. << endl
  174. << " Number of iterations = "
  175. << icount
  176. << endl << endl;
  177. return( 0 );
  178. }
  179. int bi_cg_d( Matrix& A, Vector& xi, Vector& b, double EPS ) {
  180. const int ns = A.dim_i();
  181. int maxit = 2000;
  182. int iflag = 1, icount = 0;
  183. double delta0, deltai;
  184. double bet, roi;
  185. double roim1 = 1.0, al;
  186. Vector ri, roc;
  187. Vector vi(ns,0.0), pi(ns,0.0), pih(ns,0.0);
  188. Matrix At = A.trans();
  189. ri = b - A * xi;
  190. roc = ri;
  191. delta0 = norm2( ri );
  192. while( iflag != 0 && icount < maxit ) {
  193. icount += 1;
  194. roi = roc * ri;
  195. bet = roi / roim1;
  196. pi = ri + pi * bet;
  197. pih = roc + pih * bet;
  198. vi = A * pi;
  199. al = roi / ( pih * vi );
  200. xi += pi * al;
  201. ri -= vi * al;
  202. roc -= ( At * pih ) * al;
  203. deltai = norm2( ri );
  204. if( deltai < EPS * delta0 )
  205. iflag = 0;
  206. else
  207. roim1 = roi;
  208. }
  209. cerr << endl
  210. << endl
  211. << " Number of iterations = "
  212. << icount
  213. << endl << endl;
  214. return( 0 );
  215. }
  216. int bi_cg_c( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS ) {
  217. const int ns = A.dim_i();
  218. int maxit = 2000;
  219. int iflag = 1, icount = 0;
  220. double delta0, deltai;
  221. Complex bet, roi;
  222. Complex roim1 = 1.0, al;
  223. CmplxVector ri, roc;
  224. CmplxVector vi(ns,0.0), pi(ns,0.0), pih(ns,0.0);
  225. CmplxMatrix At = A.trans();
  226. ri = b - A * xi;
  227. roc = ri;
  228. delta0 = norm2( ri );
  229. while( iflag != 0 && icount < maxit ) {
  230. icount += 1;
  231. roi = inner( roc, ri );
  232. bet = roi / roim1;
  233. pi = ri + pi * bet;
  234. pih = roc + pih * bet;
  235. vi = A * pi;
  236. al = roi / inner( pih, vi );
  237. xi += pi * al;
  238. ri -= vi * al;
  239. roc -= ( At * pih ) * al;
  240. deltai = norm2( ri );
  241. if( deltai < EPS * delta0 )
  242. iflag = 0;
  243. else
  244. roim1 = roi;
  245. }
  246. cerr << endl
  247. << endl
  248. << " Number of iterations = "
  249. << icount
  250. << endl << endl;
  251. return( 0 );
  252. }
  253. int bi_cgstab_d_precond( Matrix& A, Vector& xi, Vector& b, double EPS ) {
  254. const int ns = A.dim_i();
  255. int i, maxit = 2000;
  256. int iflag = 1, icount = 0;
  257. double delta0, deltai;
  258. double bet, roi;
  259. double roim1 = 1.0, al = 1.0, wi = 1.0;
  260. Vector mp(ns), y(ns), z(ns), ri, roc;
  261. Vector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
  262. for( i = 0; i < ns; i++ )
  263. mp[i] = 1.0 / A(i,i);
  264. ri = b - A * xi;
  265. roc = ri;
  266. delta0 = norm2( ri );
  267. while( iflag != 0 && icount < maxit ) {
  268. icount += 1;
  269. roi = roc * ri;
  270. bet = ( roi / roim1 ) * ( al / wi );
  271. pi = ri + ( pi - vi * wi ) * bet;
  272. y = mp ^ pi;
  273. vi = A * y;
  274. al = roi / ( roc * vi );
  275. s = ri - vi * al;
  276. z = mp ^ s;
  277. t = A * z;
  278. wi = ( t * s ) / ( t * t );
  279. xi += y * al + z * wi;
  280. ri = s - t * wi;
  281. deltai = norm2( ri );
  282. if( deltai < EPS * delta0 )
  283. iflag = 0;
  284. else
  285. roim1 = roi;
  286. }
  287. cerr << endl
  288. << endl
  289. << " Number of iterations = "
  290. << icount
  291. << endl << endl;
  292. return( 0 );
  293. }
  294. int bi_cgstab_c_precond( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS ) {
  295. const int ns = A.dim_i();
  296. int i, maxit = 2000;
  297. int iflag = 1, icount = 0;
  298. double delta0, deltai;
  299. Complex bet, roi;
  300. Complex roim1 = 1.0, al = 1.0, wi = 1.0;
  301. CmplxVector mp(ns), y(ns), z(ns), ri, roc;
  302. CmplxVector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
  303. for( i = 0; i < ns; i++ )
  304. mp[i] = 1.0 / A(i,i);
  305. ri = b - A * xi;
  306. roc = ri;
  307. delta0 = norm2( ri );
  308. while( iflag != 0 && icount < maxit ) {
  309. icount += 1;
  310. roi = inner( roc, ri );
  311. bet = ( roi / roim1 ) * ( al / wi );
  312. pi = ri + ( pi - vi * wi ) * bet;
  313. y = mp ^ pi;
  314. vi = A * y;
  315. al = roi / inner( roc, vi );
  316. s = ri - vi * al;
  317. z = mp ^ s;
  318. t = A * z;
  319. wi = inner( t, s ) / inner( t, t );
  320. xi += y * al + z * wi;
  321. ri = s - t * wi;
  322. deltai = norm2( ri );
  323. if( deltai < EPS * delta0 )
  324. iflag = 0;
  325. else
  326. roim1 = roi;
  327. }
  328. cerr << endl
  329. << endl
  330. << " Number of iterations = "
  331. << icount
  332. << endl << endl;
  333. return( 0 );
  334. }
  335. int bi_cgstab_d_sparse_precond( Matrix& A, Vector& xi, Vector& b, double EPS, double thresh ) {
  336. const int ns = A.dim_i();
  337. int i, maxit = 2000;
  338. int iflag = 1, icount = 0;
  339. double delta0, deltai;
  340. double bet, roi;
  341. double roim1 = 1.0, al = 1.0, wi = 1.0;
  342. Vector mp(ns), y(ns), z(ns), ri, roc;
  343. Vector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
  344. for( i = 0; i < ns; i++ )
  345. mp[i] = 1.0 / A(i,i);
  346. Vector sa;
  347. IVectorl ija;
  348. sprsin_d( A, thresh, sa, ija );
  349. ri = b - sprsax_d_v( sa, ija, xi ); ;
  350. roc = ri;
  351. delta0 = norm2( ri );
  352. while( iflag != 0 && icount < maxit ) {
  353. icount += 1;
  354. roi = roc * ri;
  355. bet = ( roi / roim1 ) * ( al / wi );
  356. pi = ri + ( pi - vi * wi ) * bet;
  357. y = mp ^ pi;
  358. vi = sprsax_d_v( sa, ija, y );
  359. al = roi / ( roc * vi );
  360. s = ri - vi * al;
  361. z = mp ^ s;
  362. t = sprsax_d_v( sa, ija, z ); ;
  363. wi = ( t * s ) / ( t * t );
  364. xi += y * al + z * wi;
  365. ri = s - t * wi;
  366. deltai = norm2( ri );
  367. if( deltai < EPS * delta0 )
  368. iflag = 0;
  369. else
  370. roim1 = roi;
  371. }
  372. cerr << endl
  373. << endl
  374. << " Number of iterations = "
  375. << icount
  376. << endl << endl;
  377. return( 0 );
  378. }
  379. int bi_cgstab_c_sparse_precond( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS, double thresh ) {
  380. const int ns = A.dim_i();
  381. int i, maxit = 2000;
  382. int iflag = 1, icount = 0;
  383. double delta0, deltai;
  384. Complex bet, roi;
  385. Complex roim1 = 1.0, al = 1.0, wi = 1.0;
  386. CmplxVector mp(ns), y(ns), z(ns), ri, roc;
  387. CmplxVector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
  388. for( i = 0; i < ns; i++ )
  389. mp[i] = 1.0 / A(i,i);
  390. CmplxVector sa;
  391. IVectorl ija;
  392. sprsin_c( A, thresh, sa, ija );
  393. ri = b - sprsax_c_v( sa, ija, xi );
  394. roc = ri;
  395. delta0 = norm2( ri );
  396. while( iflag != 0 && icount < maxit ) {
  397. icount += 1;
  398. roi = inner( roc, ri );
  399. bet = ( roi / roim1 ) * ( al / wi );
  400. pi = ri + ( pi - vi * wi ) * bet;
  401. y = mp ^ pi;
  402. vi = sprsax_c_v( sa, ija, y );
  403. al = roi / inner( roc, vi );
  404. s = ri - vi * al;
  405. z = mp ^ s;
  406. t = sprsax_c_v( sa, ija, z );
  407. wi = inner( t, s ) / inner( t, t );
  408. xi += y * al + z * wi;
  409. ri = s - t * wi;
  410. deltai = norm2( ri );
  411. if( deltai < EPS * delta0 )
  412. iflag = 0;
  413. else
  414. roim1 = roi;
  415. }
  416. cerr << endl
  417. << endl
  418. << " Number of iterations = "
  419. << icount
  420. << endl << endl;
  421. return( 0 );
  422. }