update clcg

example/ex7.cpp

@@ -123,6 +123,14 @@ int main(int argc, char const *argv[])
 
     m.assign_all(std::complex<double>(0.0, 0.0));
     test.CLCG_Minimize(m, B, gctl::CLCG_BICG);
+    std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
+
+    m.assign_all(std::complex<double>(0.0, 0.0));
+    test.CLCG_Minimize(m, B, gctl::CLCG_CGS);
+    std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
+
+    m.assign_all(std::complex<double>(0.0, 0.0));
+    test.CLCG_Minimize(m, B, gctl::CLCG_BICGSTAB);
     std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
 	return 0;
 }

lib/optimization/clcg.cpp

@@ -30,7 +30,7 @@
 /**
  * Default parameter for conjugate gradient methods
  */
-static const gctl::clcg_para clcg_defparam = {0, 1e-8, 0};
+static const gctl::clcg_para clcg_defparam = {0, 1e-8, 0, 1e-8};
 
 gctl::clcg_solver::clcg_solver()
 {
@@ -259,13 +259,13 @@ void gctl::clcg_solver::CLCG_Minimize(array<std::complex<double> > &m, const arr
         switch (solver_id)
 		{
 			case CLCG_BICG:
-				ss << "Solver: " << std::setw(9) << "Bi-CG, Times cost: " << costime << " ms" << std::endl; break;
+				ss << "Solver: " << std::setw(9) << "BiCG, Times cost: " << costime << " ms" << std::endl; break;
 			case CLCG_BICG_SYM:
-				ss << "Solver: " << std::setw(9) << "Bi-CG (symmetrically accelerated), Times cost: " << costime << " ms" << std::endl; break;
+				ss << "Solver: " << std::setw(9) << "BiCG-Sym, Times cost: " << costime << " ms" << std::endl; break;
 			case CLCG_CGS:
 				ss << "Solver: " << std::setw(9) << "CGS, Times cost: " << costime << " ms" << std::endl; break;
             case CLCG_BICGSTAB:
-				ss << "Solver: " << std::setw(9) << "CGS, Times cost: " << costime << " ms" << std::endl; break;
+				ss << "Solver: " << std::setw(9) << "BiCG-Stab, Times cost: " << costime << " ms" << std::endl; break;
 			case CLCG_TFQMR:
 				ss << "Solver: " << std::setw(9) << "TFQMR, Times cost: " << costime << " ms" << std::endl; break;
 			default:
@@ -462,12 +462,219 @@ void gctl::clcg_solver::clbicg_symmetric(array<std::complex<double> > &m, const
 
 void gctl::clcg_solver::clcgs(array<std::complex<double> > &m, const array<std::complex<double> > &B, std::ostream &ss)
 {
-	return;
+    size_t n_size = B.size();
+	//check parameters
+	if (n_size <= 0) return clcg_error_str(CLCG_INVILAD_VARIABLE_SIZE, ss);
+	if (clcg_param_.max_iterations < 0) return clcg_error_str(CLCG_INVILAD_MAX_ITERATIONS, ss);
+	if (clcg_param_.epsilon <= 0.0 || clcg_param_.epsilon >= 1.0) return clcg_error_str(CLCG_INVILAD_EPSILON, ss);
+
+    r1k.resize(n_size); r2k.resize(n_size); d1k.resize(n_size);
+    uk.resize(n_size); qk.resize(n_size); wk.resize(n_size);
+	Ax.resize(n_size);
+
+	CLCG_Ax(m, Ax, gctl::NoTrans, gctl::NoConj);
+
+    std::complex<double> one_z(1.0, 0.0);
+    std::complex<double> one_one(1.0, 1.0);
+    vecdiff(r1k, B, Ax, one_z, one_z);
+    veccpy(uk, r1k, one_z);
+    veccpy(d1k, r1k, one_z);
+
+    srand(time(0));
+	std::complex<double> rhok;
+	do
+	{
+        for (size_t i = 0; i < n_size; i++)
+        {
+            r2k[i] = random(one_z, 2.0*one_one);
+        }
+
+        rhok = vecinner(r2k, r1k);
+	} while (std::norm(rhok) < clcg_param_.lamda);
+
+	double r0_square, rk_square;
+	std::complex<double> r0_mod, rk_mod;
+
+	rk_mod = vecinner(r1k, r1k);
+	r0_square = rk_square = std::norm(rk_mod);
+	if (r0_square < 1.0) r0_square = 1.0;
+
+	if (clcg_param_.abs_diff && sqrt(rk_square)/n_size <= clcg_param_.epsilon)
+	{
+        CLCG_Progress(m, sqrt(rk_square)/n_size, clcg_param_, 0, ss);
+		return clcg_error_str(CLCG_ALREADY_OPTIMIZIED, ss);
+	}	
+	else if (rk_square/r0_square <= clcg_param_.epsilon)
+	{
+        CLCG_Progress(m, rk_square/r0_square, clcg_param_, 0, ss);
+		return clcg_error_str(CLCG_ALREADY_OPTIMIZIED, ss);
+	}
+
+	double residual;
+    std::complex<double> ak, rhok2, sigma, betak;
+
+    size_t t = 0;
+	while(1)
+	{
+		if (clcg_param_.abs_diff) residual = sqrt(rk_square)/n_size;
+		else residual = rk_square/r0_square;
+
+        if (CLCG_Progress(m, residual, clcg_param_, t, ss))
+        {
+            return clcg_error_str(CLCG_STOP, ss);
+        }
+
+		if (residual <= clcg_param_.epsilon)
+		{
+			return clcg_error_str(CLCG_CONVERGENCE, ss);
+		}
+
+		if (clcg_param_.max_iterations > 0 && t+1 > clcg_param_.max_iterations)
+		{
+			return clcg_error_str(CLCG_REACHED_MAX_ITERATIONS, ss);
+		}
+		
+		t++;
+
+		CLCG_Ax(d1k, Ax, gctl::NoTrans, gctl::NoConj); // vk = Apk
+		sigma = vecinner(r2k, Ax);
+		ak = rhok/sigma;
+
+        vecdiff(qk, uk, Ax, one_z, ak);
+        vecadd(wk, uk, qk, one_z, one_z);
+
+		CLCG_Ax(wk, Ax, gctl::NoTrans, gctl::NoConj);
+
+        vecapp(m, wk, ak);
+        vecsub(r1k, Ax, ak);
+
+		rk_mod = vecinner(r1k, r1k);
+		rk_square = std::norm(rk_mod);
+
+		if (!vecvalid(m))
+        {
+            return clcg_error_str(CLCG_NAN_VALUE, ss);
+        }
+
+		rhok2 = vecinner(r2k, r1k);
+		betak = rhok2/rhok;
+		rhok = rhok2;
+
+        vecadd(uk, r1k, qk, one_z, betak);
+        vecadd(d1k, qk, d1k, one_z, betak);
+        vecadd(d1k, uk, d1k, one_z, betak);
+	}
+
+	return clcg_error_str(CLCG_UNKNOWN_ERROR, ss);
 }
 
 void gctl::clcg_solver::clbicgstab(array<std::complex<double> > &m, const array<std::complex<double> > &B, std::ostream &ss)
 {
-	return;
+	size_t n_size = B.size();
+	//check parameters
+	if (n_size <= 0) return clcg_error_str(CLCG_INVILAD_VARIABLE_SIZE, ss);
+	if (clcg_param_.max_iterations < 0) return clcg_error_str(CLCG_INVILAD_MAX_ITERATIONS, ss);
+	if (clcg_param_.epsilon <= 0.0 || clcg_param_.epsilon >= 1.0) return clcg_error_str(CLCG_INVILAD_EPSILON, ss);
+	
+    r1k.resize(n_size); r2k.resize(n_size);
+    d1k.resize(n_size); d2k.resize(n_size);
+	Ax.resize(n_size); Ad.resize(n_size);
+
+	CLCG_Ax(m, Ax, gctl::NoTrans, gctl::NoConj);
+
+    std::complex<double> one_z(1.0, 0.0);
+    std::complex<double> one_one(1.0, 1.0);
+    vecdiff(r1k, B, Ax, one_z, one_z);
+    veccpy(d1k, r1k, one_z);
+
+    srand(time(0));
+	std::complex<double> rhok;
+	do
+	{
+        for (size_t i = 0; i < n_size; i++)
+        {
+            r2k[i] = random(one_z, 2.0*one_one);
+        }
+
+        rhok = vecinner(r2k, r1k);
+	} while (std::norm(rhok) < clcg_param_.lamda);
+
+	double r0_square, rk_square;
+	std::complex<double> r0_mod, rk_mod;
+	rk_mod = vecinner(r1k, r1k);
+	r0_square = rk_square = std::norm(rk_mod);
+
+	if (r0_square < 1.0) r0_square = 1.0;
+
+	if (clcg_param_.abs_diff && sqrt(rk_square)/n_size <= clcg_param_.epsilon)
+	{
+        CLCG_Progress(m, sqrt(rk_square)/n_size, clcg_param_, 0, ss);
+		return clcg_error_str(CLCG_ALREADY_OPTIMIZIED, ss);
+	}	
+	else if (rk_square/r0_square <= clcg_param_.epsilon)
+	{
+        CLCG_Progress(m, rk_square/r0_square, clcg_param_, 0, ss);
+		return clcg_error_str(CLCG_ALREADY_OPTIMIZIED, ss);
+	}
+
+    double residual;
+    std::complex<double> ak, rhok2, sigma, omega, betak, Ass, AsAs;
+
+	size_t t = 0;
+	while(1)
+	{
+		if (clcg_param_.abs_diff) residual = sqrt(rk_square)/n_size;
+		else residual = rk_square/r0_square;
+
+        if (CLCG_Progress(m, residual, clcg_param_, t, ss))
+        {
+            return clcg_error_str(CLCG_STOP, ss);
+        }
+
+		if (residual <= clcg_param_.epsilon)
+		{
+			return clcg_error_str(CLCG_CONVERGENCE, ss);
+		}
+
+		if (clcg_param_.max_iterations > 0 && t+1 > clcg_param_.max_iterations)
+		{
+			return clcg_error_str(CLCG_REACHED_MAX_ITERATIONS, ss);
+		}
+		
+		t++;
+
+		CLCG_Ax(d1k, Ax, gctl::NoTrans, gctl::NoConj);
+		sigma = vecinner(r2k, Ax);
+		ak = rhok/sigma;
+
+        vecdiff(d2k, r1k, Ax, one_z, ak);
+
+		CLCG_Ax(d2k, Ad, gctl::NoTrans, gctl::NoConj);
+		Ass = vecinner(Ad, d2k);
+		AsAs = vecinner(Ad, Ad);
+		omega = Ass/AsAs;
+
+        vecdiff(r1k, d2k, Ad, one_z, omega);
+        vecadd(d2k, d1k, d2k, ak, omega);
+        vecapp(m, d2k, one_z);
+
+		rk_mod = vecinner(r1k, r1k);
+		rk_square = std::norm(rk_mod);
+
+		if (!vecvalid(m))
+        {
+            return clcg_error_str(CLCG_NAN_VALUE, ss);
+        }
+
+		rhok2 = vecinner(r2k, r1k);
+		betak = rhok2*ak/(rhok*omega);
+		rhok = rhok2;
+
+        vecdiff(d1k, d1k, Ax, one_z, omega);
+        vecadd(d1k, r1k, d1k, one_z, betak);
+	}
+
+	return clcg_error_str(CLCG_UNKNOWN_ERROR, ss);
 }
 
 void gctl::clcg_solver::cltfqmr(array<std::complex<double> > &m, const array<std::complex<double> > &B, std::ostream &ss)

lib/optimization/clcg.h

@@ -132,6 +132,13 @@ namespace gctl
          * applied to the non-constrained methods.
          */
         int abs_diff;
+
+        /**
+         * Minimal value for testing if a vector's module is bigger than the threshold.
+         * The default value is 1e-8 
+         * 
+         */
+        double lamda;
     };
 
     class clcg_solver
@@ -143,7 +150,8 @@ namespace gctl
 
         // make them class variables are more suitable for repetitively usages
         array<std::complex<double> > r1k, r2k, d1k, d2k;
-	    array<std::complex<double> > Ax;
+        array<std::complex<double> > uk, qk, wk;
+	    array<std::complex<double> > Ax, Ad;
 
         void clcg_error_str(clcg_return_code err_code, std::ostream &ss);