123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566 |
- /*****************************************************************/
- /*** ***/
- /*** Bi_CGSTAB iterative method ***/
- /*** ***/
- /*** A - system matrix, xi - initial guess and output vector, ***/
- /*** b - right hand side ***/
- /*** ***/
- /*****************************************************************/
- //#include "stdafx.h"
- #include <iostream.h>
- #include <math.h>
- #include "vector.h"
- #include "matrix.h"
- #include "complex.h"
- #include "cmplxvec.h"
- #include "cmplxmat.h"
- #include "ivectorl.h"
- #include "sparse.h"
- #include "bicgstab.h"
- int bi_cgstab_d( Matrix& A, Vector& xi, Vector& b, double EPS ) {
- const int ns = A.dim_i();
- int maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- double bet, roi;
- double roim1 = 1.0, al = 1.0, wi = 1.0;
- Vector ri, roc;
- Vector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
-
- ri = b - A * xi;
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = roc * ri;
- bet = ( roi / roim1 ) * ( al / wi );
- pi = ri + ( pi - vi * wi ) * bet;
- vi = A * pi;
- al = roi / ( roc * vi );
- s = ri - vi * al;
- t = A * s;
- wi = ( t * s ) / ( t * t );
- xi += pi * al + s * wi;
- ri = s - t * wi;
- deltai = norm2( ri );
- // cerr << " " << icount << " "
- // << setprecision(10)
- // << deltai << " "
- // << delta0 << " "
- // << deltai/delta0 << endl;
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cout << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cgstab_c( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS ) {
- const int ns = A.dim_i();
- int maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- Complex bet, roi;
- Complex roim1 = 1.0, al = 1.0, wi = 1.0;
- CmplxVector ri, roc;
- CmplxVector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
-
- ri = b - A * xi;
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = inner( roc, ri );
- bet = ( roi / roim1 ) * ( al / wi );
- pi = ri + ( pi - vi * wi ) * bet;
- vi = A * pi;
- al = roi / inner( roc, vi );
- s = ri - vi * al;
- t = A * s;
- wi = inner( t, s ) / inner( t, t );
- xi += pi * al + s * wi;
- ri = s - t * wi;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cgstab_d_sparse( Matrix& A, Vector& xi, Vector& b, double EPS, double thresh ) {
- const int ns = A.dim_i();
- int maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- double bet, roi;
- double roim1 = 1.0, al = 1.0, wi = 1.0;
- Vector ri, roc;
- Vector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
- Vector sa;
- IVectorl ija;
-
- sprsin_d( A, thresh, sa, ija );
- ri = b - sprsax_d_v( sa, ija, xi );
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = roc * ri;
- bet = ( roi / roim1 ) * ( al / wi );
- pi = ri + ( pi - vi * wi ) * bet;
- vi = sprsax_d_v( sa, ija, pi );
- al = roi / ( roc * vi );
- s = ri - vi * al;
- t = sprsax_d_v( sa, ija, s );
- wi = ( t * s ) / ( t * t );
- xi += pi * al + s * wi;
- ri = s - t * wi;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cgstab_c_sparse( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS, double thresh ) {
- const int ns = A.dim_i();
- int maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- Complex bet, roi;
- Complex roim1 = 1.0, al = 1.0, wi = 1.0;
- CmplxVector ri, roc;
- CmplxVector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
- CmplxVector sa;
- IVectorl ija;
-
- sprsin_c( A, thresh, sa, ija );
-
- ri = b - sprsax_c_v( sa, ija, xi );
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = inner( roc, ri );
- bet = ( roi / roim1 ) * ( al / wi );
- pi = ri + ( pi - vi * wi ) * bet;
- vi = sprsax_c_v( sa, ija, pi );
- al = roi / inner( roc, vi );
- s = ri - vi * al;
- t = sprsax_c_v( sa, ija, s );
- wi = inner( t, s ) / inner( t, t );
- xi += pi * al + s * wi;
- ri = s - t * wi;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cg_d( Matrix& A, Vector& xi, Vector& b, double EPS ) {
- const int ns = A.dim_i();
- int maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- double bet, roi;
- double roim1 = 1.0, al;
- Vector ri, roc;
- Vector vi(ns,0.0), pi(ns,0.0), pih(ns,0.0);
- Matrix At = A.trans();
-
- ri = b - A * xi;
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
-
- icount += 1;
-
- roi = roc * ri;
- bet = roi / roim1;
- pi = ri + pi * bet;
- pih = roc + pih * bet;
- vi = A * pi;
- al = roi / ( pih * vi );
- xi += pi * al;
- ri -= vi * al;
- roc -= ( At * pih ) * al;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cg_c( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS ) {
- const int ns = A.dim_i();
- int maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- Complex bet, roi;
- Complex roim1 = 1.0, al;
- CmplxVector ri, roc;
- CmplxVector vi(ns,0.0), pi(ns,0.0), pih(ns,0.0);
- CmplxMatrix At = A.trans();
-
- ri = b - A * xi;
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = inner( roc, ri );
- bet = roi / roim1;
- pi = ri + pi * bet;
- pih = roc + pih * bet;
- vi = A * pi;
- al = roi / inner( pih, vi );
- xi += pi * al;
- ri -= vi * al;
- roc -= ( At * pih ) * al;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cgstab_d_precond( Matrix& A, Vector& xi, Vector& b, double EPS ) {
- const int ns = A.dim_i();
- int i, maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- double bet, roi;
- double roim1 = 1.0, al = 1.0, wi = 1.0;
- Vector mp(ns), y(ns), z(ns), ri, roc;
- Vector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
-
- for( i = 0; i < ns; i++ )
- mp[i] = 1.0 / A(i,i);
-
- ri = b - A * xi;
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = roc * ri;
- bet = ( roi / roim1 ) * ( al / wi );
- pi = ri + ( pi - vi * wi ) * bet;
- y = mp ^ pi;
- vi = A * y;
- al = roi / ( roc * vi );
- s = ri - vi * al;
- z = mp ^ s;
- t = A * z;
- wi = ( t * s ) / ( t * t );
- xi += y * al + z * wi;
- ri = s - t * wi;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cgstab_c_precond( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS ) {
- const int ns = A.dim_i();
- int i, maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- Complex bet, roi;
- Complex roim1 = 1.0, al = 1.0, wi = 1.0;
- CmplxVector mp(ns), y(ns), z(ns), ri, roc;
- CmplxVector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
-
- for( i = 0; i < ns; i++ )
- mp[i] = 1.0 / A(i,i);
-
- ri = b - A * xi;
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = inner( roc, ri );
- bet = ( roi / roim1 ) * ( al / wi );
- pi = ri + ( pi - vi * wi ) * bet;
- y = mp ^ pi;
- vi = A * y;
- al = roi / inner( roc, vi );
- s = ri - vi * al;
- z = mp ^ s;
- t = A * z;
- wi = inner( t, s ) / inner( t, t );
- xi += y * al + z * wi;
- ri = s - t * wi;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cgstab_d_sparse_precond( Matrix& A, Vector& xi, Vector& b, double EPS, double thresh ) {
- const int ns = A.dim_i();
- int i, maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- double bet, roi;
- double roim1 = 1.0, al = 1.0, wi = 1.0;
- Vector mp(ns), y(ns), z(ns), ri, roc;
- Vector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
-
- for( i = 0; i < ns; i++ )
- mp[i] = 1.0 / A(i,i);
- Vector sa;
- IVectorl ija;
-
- sprsin_d( A, thresh, sa, ija );
- ri = b - sprsax_d_v( sa, ija, xi ); ;
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = roc * ri;
- bet = ( roi / roim1 ) * ( al / wi );
- pi = ri + ( pi - vi * wi ) * bet;
- y = mp ^ pi;
- vi = sprsax_d_v( sa, ija, y );
- al = roi / ( roc * vi );
- s = ri - vi * al;
- z = mp ^ s;
- t = sprsax_d_v( sa, ija, z ); ;
- wi = ( t * s ) / ( t * t );
- xi += y * al + z * wi;
- ri = s - t * wi;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
- int bi_cgstab_c_sparse_precond( CmplxMatrix& A, CmplxVector& xi, CmplxVector& b, double EPS, double thresh ) {
- const int ns = A.dim_i();
- int i, maxit = 2000;
- int iflag = 1, icount = 0;
- double delta0, deltai;
- Complex bet, roi;
- Complex roim1 = 1.0, al = 1.0, wi = 1.0;
- CmplxVector mp(ns), y(ns), z(ns), ri, roc;
- CmplxVector s(ns,0.0), t(ns,0.0), vi(ns,0.0), pi(ns,0.0);
- for( i = 0; i < ns; i++ )
- mp[i] = 1.0 / A(i,i);
- CmplxVector sa;
- IVectorl ija;
-
- sprsin_c( A, thresh, sa, ija );
-
- ri = b - sprsax_c_v( sa, ija, xi );
- roc = ri;
- delta0 = norm2( ri );
-
- while( iflag != 0 && icount < maxit ) {
- icount += 1;
-
- roi = inner( roc, ri );
- bet = ( roi / roim1 ) * ( al / wi );
- pi = ri + ( pi - vi * wi ) * bet;
- y = mp ^ pi;
- vi = sprsax_c_v( sa, ija, y );
- al = roi / inner( roc, vi );
- s = ri - vi * al;
- z = mp ^ s;
- t = sprsax_c_v( sa, ija, z );
- wi = inner( t, s ) / inner( t, t );
- xi += y * al + z * wi;
- ri = s - t * wi;
- deltai = norm2( ri );
- if( deltai < EPS * delta0 )
- iflag = 0;
- else
- roim1 = roi;
- }
- cerr << endl
- << endl
- << " Number of iterations = "
- << icount
- << endl << endl;
- return( 0 );
- }
|