itersolv.c

OpenFVM-v1.1 open source cfd code

	    if (PrecondProc != NULL)
		(*PrecondProc)(A, &y, &v_, OmegaPrecond);
	    else
		Asgn_VV(&y, &v_);
	    if (Q_KerDefined(A))
	        OrthoRightKer_VQ(&y, A);
            RhoNew = l2Norm_V(&y);
            Asgn_VV(&w_, &r);
            Asgn_VV(&z, &w_);
            Xi = l2Norm_V(&z);
            Gamma = 1.0;
            Eta = - 1.0;
            GammaOld = Gamma;
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, QMRIterId)
                && Iter < MaxIter) {
                Iter++;
                Rho = RhoNew;
                if (IsZero(Rho) || IsZero(Xi)) {
                    LASError(LASBreakdownErr, "QMRIter", "Rho", "Xi", NULL);
                    break;
                }
                Asgn_VV(&v, Mul_SV(1.0 / Rho, &v_));
                MulAsgn_VS(&y, 1.0 / Rho);
                Asgn_VV(&w, Mul_SV(1.0 / Xi, &w_));
                MulAsgn_VS(&z, 1.0 / Xi);
                Delta = Mul_VV(&z, &y);
                if (IsZero(Delta)) {
                    LASError(LASBreakdownErr, "QMRIter", "Delta", NULL, NULL);
                    break;
                }
                Asgn_VV(&y_, &y);
	        if (PrecondProc != NULL)
                    (*PrecondProc)(Transp_Q(A), &z_, &z, OmegaPrecond);
	        else
		    Asgn_VV(&z_, &z);
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&z_, Transp_Q(A));
                if (Iter == 1) {
                    Asgn_VV(&p, &y_);
                    Asgn_VV(&q, &z_);
                } else {
                    Asgn_VV(&p, Sub_VV(&y_, Mul_SV(Xi * Delta / Epsilon, &p)));
                    Asgn_VV(&q, Sub_VV(&z_, Mul_SV(Rho * Delta / Epsilon, &q)));
                }
                Asgn_VV(&p_, Mul_QV(A, &p));
                Epsilon = Mul_VV(&q, &p_);
                if (IsZero(Epsilon)) {
                    LASError(LASBreakdownErr, "QMRIter", "Epsilon", NULL, NULL);
                    break;
                }
                Beta = Epsilon / Delta;
                if (IsZero(Beta)) {
                    LASError(LASBreakdownErr, "QMRIter", "Beta", NULL, NULL);
                    break;
                }
                Asgn_VV(&v_, Sub_VV(&p_, Mul_SV(Beta, &v)));
	        if (PrecondProc != NULL)
                    (*PrecondProc)(A, &y, &v_, OmegaPrecond);
	        else
		    Asgn_VV(&y, &v_);
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&y, A);
                RhoNew = l2Norm_V(&y);
                Asgn_VV(&w_, Sub_VV(Mul_QV(Transp_Q(A), &q), Mul_SV(Beta, &w)));
                Asgn_VV(&z, &w_);
                Xi = l2Norm_V(&z);
                Theta = RhoNew / (GammaOld * fabs(Beta));
                Gamma = 1.0 / sqrt(1.0 + Theta * Theta);
                if (IsZero(Gamma)) {
                    LASError(LASBreakdownErr, "QMRIter", "Gamma", NULL, NULL);
                    break;
                }
                Eta = - Eta * Rho * Gamma * Gamma / (Beta * GammaOld * GammaOld);
                if (Iter == 1) {
                    Asgn_VV(&d, Mul_SV(Eta, &p));
                    Asgn_VV(&s, Mul_SV(Eta, &p_));
                } else {
                    Asgn_VV(&d, Add_VV(Mul_SV(Eta, &p), 
                        Mul_SV(ThetaOld * ThetaOld * Gamma * Gamma, &d)));
                    Asgn_VV(&s, Add_VV(Mul_SV(Eta, &p_), 
                        Mul_SV(ThetaOld * ThetaOld * Gamma * Gamma, &s)));
                }
                AddAsgn_VV(x, &d);
                SubAsgn_VV(&r, &s);
                GammaOld = Gamma;
                ThetaOld = Theta;
            }
        } else {
            /* plain QMR (M1 ~ A, M2 = I, y = y_ = v_, z = z_ = w_) */
            Asgn_VV(&v_, &r);
            RhoNew = l2Norm_V(&v_);
            Asgn_VV(&w_, &r);
            Xi = l2Norm_V(&w_);
            Gamma = 1.0;
            Eta = - 1.0;
            GammaOld = Gamma;
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, QMRIterId)
                && Iter < MaxIter) {
                Iter++;
                Rho = RhoNew;
                if (IsZero(Rho) || IsZero(Xi)) {
                    LASError(LASBreakdownErr, "QMRIter", "Rho", "Xi", NULL);
                    break;
                }
                Asgn_VV(&v, Mul_SV(1.0 / Rho, &v_));
                MulAsgn_VS(&v_, 1.0 / Rho);
                Asgn_VV(&w, Mul_SV(1.0 / Xi, &w_));
                MulAsgn_VS(&w_, 1.0 / Xi);
                Delta = Mul_VV(&w_, &v_);
                if (IsZero(Delta)) {
                    LASError(LASBreakdownErr, "QMRIter", "Rho", "Xi", NULL);
                    break;
                }
                if (Iter == 1) {
                    Asgn_VV(&p, &v_);
                    Asgn_VV(&q, &w_);
                } else {
                    Asgn_VV(&p, Sub_VV(&v_, Mul_SV(Xi * Delta / Epsilon, &p)));
                    Asgn_VV(&q, Sub_VV(&w_, Mul_SV(Rho * Delta / Epsilon, &q)));
                }
                Asgn_VV(&p_, Mul_QV(A, &p));
                Epsilon = Mul_VV(&q, &p_);
                if (IsZero(Epsilon)) {
                    LASError(LASBreakdownErr, "QMRIter", "Epsilon", NULL, NULL);
                    break;
                }
                Beta = Epsilon / Delta;
                if (IsZero(Beta)) {
                    LASError(LASBreakdownErr, "QMRIter", "Beta", NULL, NULL);
                    break;
                }
                Asgn_VV(&v_, Sub_VV(&p_, Mul_SV(Beta, &v)));
                RhoNew = l2Norm_V(&v);
                Asgn_VV(&w_, Sub_VV(Mul_QV(Transp_Q(A), &q), Mul_SV(Beta, &w)));
                Xi = l2Norm_V(&w_);
                Theta = RhoNew / (GammaOld * fabs(Beta));
                Gamma = 1.0 / sqrt(1.0 + Theta * Theta);
                if (IsZero(Gamma)) {
                    LASError(LASBreakdownErr, "QMRIter", "Gamma", NULL, NULL);
                    break;
                }
                Eta = - Eta * Rho * Gamma * Gamma / (Beta * GammaOld * GammaOld);
                if (Iter == 1) {
                    Asgn_VV(&d, Mul_SV(Eta, &p));
                    Asgn_VV(&s, Mul_SV(Eta, &p_));
                } else {
                    Asgn_VV(&d, Add_VV(Mul_SV(Eta, &p), 
                        Mul_SV(ThetaOld * ThetaOld * Gamma * Gamma, &d)));
                    Asgn_VV(&s, Add_VV(Mul_SV(Eta, &p_), 
                        Mul_SV(ThetaOld * ThetaOld * Gamma * Gamma, &s)));
                }
                AddAsgn_VV(x, &d);
                SubAsgn_VV(&r, &s);
		GammaOld = Gamma;
                ThetaOld = Theta;
            }
        }
        if (Q_KerDefined(A))
	   OrthoRightKer_VQ(x, A);
    }

    V_Destr(&r);
    V_Destr(&p);
    V_Destr(&p_);
    V_Destr(&q);
    V_Destr(&v);
    V_Destr(&v_);
    V_Destr(&w);
    V_Destr(&w_);
    V_Destr(&s);
    V_Destr(&d);
    if (PrecondProc != NULL || Q_KerDefined(A)) {
        V_Destr(&y);
        V_Destr(&y_);
        V_Destr(&z);
        V_Destr(&z_);
    }

    Q_Unlock(A);
    V_Unlock(x);
    V_Unlock(b);

    return(x);
}

Vector *CGSIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
            PrecondProcType PrecondProc, double OmegaPrecond)
{
    /*
     *  for details to the algorithm see
     *
     *  R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
     *  V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
     *  Templates for the Solution of Linear Systems: Building Blocks
     *  for Iterative Solvers;
     *  SIAM, Philadelphia, 1994
     *
     */

    int Iter;
    double Alpha, Beta, Rho, RhoOld = 0.0;
    double bNorm;
    size_t Dim;
    Vector r, r_, p, p_, q, u, u_, v_;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", Dim, Normal, True);
    V_Constr(&r_, "r_", Dim, Normal, True);
    V_Constr(&p, "p", Dim, Normal, True);
    V_Constr(&q, "q", Dim, Normal, True);
    V_Constr(&u, "u", Dim, Normal, True);
    V_Constr(&u_, "u_", Dim, Normal, True);
    V_Constr(&v_, "v_", Dim, Normal, True);
    if (PrecondProc != NULL || Q_KerDefined(A))
        V_Constr(&p_, "p_", Dim, Normal, True);

    if (LASResult() == LASOK) {
        bNorm = l2Norm_V(b);

        Iter = 0;
        /* r = b - A * x(i) */
        if (!IsZero(l1Norm_V(x) / Dim)) {
            if (Q_KerDefined(A))
	       OrthoRightKer_VQ(x, A);
            Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
        } else {
            Asgn_VV(&r, b);
	}
        if (PrecondProc != NULL || Q_KerDefined(A)) {
            /* preconditioned CGS */
            Asgn_VV(&r_, &r);
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGSIterId)
                && Iter < MaxIter) {
                Iter++;
                Rho = Mul_VV(&r_, &r);
                if (IsZero(Rho)) {
                    LASError(LASBreakdownErr, "CGSIter", "Rho", NULL, NULL);
                    break;
                }
                if (Iter == 1) {
                    Asgn_VV(&u, &r);
                    Asgn_VV(&p, &u);
                } else {
                    Beta = Rho / RhoOld;
                    Asgn_VV(&u, Add_VV(&r, Mul_SV(Beta, &q)));
                    Asgn_VV(&p, Add_VV(&u, Mul_SV(Beta, Add_VV(&q, Mul_SV(Beta, &p)))));
                }
	        if (PrecondProc != NULL)
                    (*PrecondProc)(A, &p_, &p, OmegaPrecond);
	        else
		    Asgn_VV(&p_, &p);
		if (Q_KerDefined(A))
		    OrthoRightKer_VQ(&p_, A);
                Asgn_VV(&v_, Mul_QV(A, &p_));
                Alpha = Rho / Mul_VV(&r_, &v_);
                Asgn_VV(&q, Sub_VV(&u, Mul_SV(Alpha, &v_)));
	        if (PrecondProc != NULL)
                    (*PrecondProc)(A, &u_, Add_VV(&u, &q), OmegaPrecond);
	        else
		    Asgn_VV(&u_, Add_VV(&u, &q));
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&u_, A);
                AddAsgn_VV(x, Mul_SV(Alpha, &u_));
                SubAsgn_VV(&r, Mul_SV(Alpha, Mul_QV(A, &u_)));
                RhoOld = Rho;
            }
        } else {
            /* plain CGS (p_ = p) */
            Asgn_VV(&r_, &r);
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGSIterId)
                && Iter < MaxIter) {
                Iter++;
                Rho = Mul_VV(&r_, &r);
                if (IsZero(Rho)) {
                    LASError(LASBreakdownErr, "CGSIter", "Rho", NULL, NULL);
                    break;
                }
                if (Iter == 1) {
                    Asgn_VV(&u, &r);
                    Asgn_VV(&p, &u);
                } else {
                    Beta = Rho / RhoOld;
                    Asgn_VV(&u, Add_VV(&r, Mul_SV(Beta, &q)));
                    Asgn_VV(&p, Add_VV(&u, Mul_SV(Beta, Add_VV(&q, Mul_SV(Beta, &p)))));
                }
                Asgn_VV(&v_, Mul_QV(A, &p));
                Alpha = Rho / Mul_VV(&r_, &v_);
                Asgn_VV(&q, Sub_VV(&u, Mul_SV(Alpha, &v_)));
                Asgn_VV(&u_, Add_VV(&u, &q));
                AddAsgn_VV(x, Mul_SV(Alpha, &u_));
                SubAsgn_VV(&r, Mul_SV(Alpha, Mul_QV(A, &u_)));
                RhoOld = Rho;
            }
        }
        if (Q_KerDefined(A))
	   OrthoRightKer_VQ(x, A);
    }
      
    V_Destr(&r);
    V_Destr(&r_);
    V_Destr(&p);
    V_Destr(&q);
    V_Destr(&u);
    V_Destr(&u_);
    V_Destr(&v_);
    if (PrecondProc != NULL || Q_KerDefined(A))
        V_Destr(&p_);

    Q_Unlock(A);
    V_Unlock(x);
    V_Unlock(b);

    return(x);
}

Vector *BiCGSTABIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
            PrecondProcType PrecondProc, double OmegaPrecond)
{
    /*
     *  for details to the algorithm see
     *
     *  R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
     *  V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
     *  Templates for the Solution of Linear Systems: Building Blocks
     *  for Iterative Solvers;
     *  SIAM, Philadelphia, 1994
     *
     */

    int Iter;
    double Alpha = 0.0, Beta, Rho, RhoOld = 0.0, Omega = 0.0;
    double bNorm;
    size_t Dim;
    Vector r, r_, p, p_, v, s, s_, t;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", Dim, Normal, True);
    V_Constr(&r_, "r_", Dim, Normal, True);
    V_Constr(&p, "p", Dim, Normal, True);
    V_Constr(&v, "v", Dim, Normal, True);
    V_Constr(&s, "s", Dim, Normal, True);
    V_Constr(&t, "t", Dim, Normal, True);
    
	if (PrecondProc != NULL || Q_KerDefined(A)) {
        V_Constr(&p_, "p_", Dim, Normal, True);
        V_Constr(&s_, "s_", Dim, Normal, True);
    }
    
    if (LASResult() == LASOK) {
        bNorm = l2Norm_V(b);

        Iter = 0;
        /* r = b - A * x(i) */
        if (!IsZero(l1Norm_V(x) / Dim)) {
            if (Q_KerDefined(A))
	       OrthoRightKer_VQ(x, A);
            Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
        } else {
            Asgn_VV(&r, b);
	}
        if (PrecondProc != NULL || Q_KerDefined(A)) {
            /* preconditioned BiCGStab */
            Asgn_VV(&r_, &r);
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, BiCGSTABIterId)
                && Iter < MaxIter) {
                Iter++;
                Rho = Mul_VV(&r_, &r);
                if (IsZero(Rho)) {
                    LASError(LASBreakdownErr, "BiCGSTABIter", "Rho", NULL, NULL);
                    break;
                }
                if (Iter == 1) {
                    Asgn_VV(&p, &r);
                } else {
                    Beta = (Rho / RhoOld) * (Alpha / Omega);
                    Asgn_VV(&p, Add_VV(&r, Mul_SV(Beta, Sub_VV(&p, Mul_SV(Omega, &v)))));
                }
	        if (PrecondProc != NULL)
                    (*PrecondProc)(A, &p_, &p, OmegaPrecond);
	        else
		    Asgn_VV(&p_, &p);
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&p_, A);
                Asgn_VV(&v, Mul_QV(A, &p_));
                Alpha = Rho / Mul_VV(&r_, &v);
                Asgn_VV(&s, Sub_VV(&r, Mul_SV(Alpha, &v)));
	        if (PrecondProc != NULL)
                    (*PrecondProc)(A, &s_, &s, OmegaPrecond);
	        else
		    Asgn_VV(&s_, &s);
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&s_, A);
                Asgn_VV(&t, Mul_QV(A, &s_));
                Omega = Mul_VV(&t, &s) / pow(l2Norm_V(&t), 2.0);
                AddAsgn_VV(x, Add_VV(Mul_SV(Alpha, &p_), Mul_SV(Omega, &s_)));
                Asgn_VV(&r, Sub_VV(&s, Mul_SV(Omega, &t)));
                RhoOld = Rho;
                if (IsZero(Omega))
                    break;
            }
        } else {
            /* plain BiCGStab (p_ = p) */
            Asgn_VV(&r_, &r);
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, BiCGSTABIterId)
                && Iter < MaxIter) {
                Iter++;
                Rho = Mul_VV(&r_, &r);
                if (IsZero(Rho)) {
                    LASError(LASBreakdownErr, "BiCGSTABIter", "Rho", NULL, NULL);
                    break;
                }
                if (Iter == 1) {
                    Asgn_VV(&p, &r);
                } else {
                    Beta = (Rho / RhoOld) * (Alpha / Omega);
                    Asgn_VV(&p, Add_VV(&r, Mul_SV(Beta, Sub_VV(&p, Mul_SV(Omega, &v)))));
                }
                Asgn_VV(&v, Mul_QV(A, &p));
                Alpha = Rho / Mul_VV(&r_, &v);
                Asgn_VV(&s, Sub_VV(&r, Mul_SV(Alpha, &v)));
                Asgn_VV(&t, Mul_QV(A, &s));
                Omega = Mul_VV(&t, &s) / pow(l2Norm_V(&t), 2.0);
                AddAsgn_VV(x, Add_VV(Mul_SV(Alpha, &p), Mul_SV(Omega, &s)));
                Asgn_VV(&r, Sub_VV(&s, Mul_SV(Omega, &t)));
                RhoOld = Rho;
                if (IsZero(Omega))
                    break;
            }
        }
        if (Q_KerDefined(A))
	   OrthoRightKer_VQ(x, A);
    }

    V_Destr(&r);
    V_Destr(&r_);
    V_Destr(&p);
    V_Destr(&v);
    V_Destr(&s);
    V_Destr(&t);
    if (PrecondProc != NULL || Q_KerDefined(A)) {
        V_Destr(&p_);
        V_Destr(&s_);
    }
    
    Q_Unlock(A);
    V_Unlock(x);
    V_Unlock(b);

    return(x);
}

void SetGMRESRestart(int Steps)
{
    GMRESSteps = Steps;
}