itersolv.c

OpenFVM-v1.1 open source cfd code

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", Dim, Normal, True);
    V_Constr(&p, "p", Dim, Normal, True);
    V_Constr(&q, "q", Dim, Normal, True);
    V_Constr(&z, "z", Dim, Normal, True);
    if (PrecondProc != NULL || Q_KerDefined(A))
        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 CGN */
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGNIterId)
                && Iter < MaxIter) {
                Iter++;
		if (PrecondProc != NULL) {
                    (*PrecondProc)(A, &z, &r, OmegaPrecond);
                    (*PrecondProc)(Transp_Q(A), &q, &z, OmegaPrecond);
                    Asgn_VV(&z, Mul_QV(Transp_Q(A), &q));
	        } else {
		    Asgn_VV(&z, &r);
		}
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&z, A);
                Rho = pow(l2Norm_V(&z), 2.0);
                if (Iter == 1) {
                    Asgn_VV(&p, &z);
                } else {
                    Beta = Rho / RhoOld;
                    Asgn_VV(&p, Add_VV(&z, Mul_SV(Beta, &p)));
                }
                Asgn_VV(&q, Mul_QV(A, &p));
		if (PrecondProc != NULL)
                    (*PrecondProc)(A, &s, &q, OmegaPrecond);
		else
		    Asgn_VV(&s, &q);
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&s, A);
                Alpha = Rho / pow(l2Norm_V(&s), 2.0);
                AddAsgn_VV(x, Mul_SV(Alpha, &p));
                SubAsgn_VV(&r, Mul_SV(Alpha, &q));
                RhoOld = Rho;
            }
        } else {
            /* plain CGN */
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGNIterId)
                && Iter < MaxIter) {
                Iter++;
                Asgn_VV(&z, Mul_QV(Transp_Q(A), &r));
                Rho = pow(l2Norm_V(&z), 2.0);
                if (Iter == 1) {
                    Asgn_VV(&p, &z);
                } else {
                    Beta = Rho / RhoOld;
                    Asgn_VV(&p, Add_VV(&z, Mul_SV(Beta, &p)));
                }
                Asgn_VV(&q, Mul_QV(A, &p));
                Alpha = Rho / pow(l2Norm_V(&q), 2.0);
                AddAsgn_VV(x, Mul_SV(Alpha, &p));
                SubAsgn_VV(&r, Mul_SV(Alpha, &q));
                RhoOld = Rho;
            }
        }
        if (Q_KerDefined(A))
	   OrthoRightKer_VQ(x, A);
    }
    
    V_Destr(&r);
    V_Destr(&p);
    V_Destr(&q);
    V_Destr(&z);
    if (PrecondProc != NULL || Q_KerDefined(A))
        V_Destr(&s);

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

    return(x);
}

Vector *GMRESIter(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
     *
     */

    char vName[10];
    int Iter, i, j, k;
    double h1, h2, r;
    double bNorm;
    double **h, *y, *s, *c1, *c2;
    size_t Dim;
    Boolean AllocOK;
    Vector *v;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    /* allocation of matrix H and vectors y, s, c1 and c2 */
    AllocOK = True;
    h = (double **)malloc((GMRESSteps + 1) * sizeof(double *));
    if (h == NULL) {
        AllocOK = False;
    } else {
        for (i = 1; i <= GMRESSteps; i++) {
            h[i] = (double *)malloc((GMRESSteps + 2) * sizeof(double));
            if (h[i] == NULL)
                AllocOK = False;
	}
    }
    y = (double *)malloc((GMRESSteps + 1) * sizeof(double));
    if (y == NULL)
        AllocOK = False;
    s = (double *)malloc((GMRESSteps + 2) * sizeof(double));
    if (s == NULL)
        AllocOK = False;
    c1 = (double *)malloc((GMRESSteps + 1) * sizeof(double));
    if (c1 == NULL)
        AllocOK = False;
    c2 = (double *)malloc((GMRESSteps + 1) * sizeof(double));
    if (c2 == NULL)
        AllocOK = False;

    /* ... and vectors u */    
    Dim = Q_GetDim(A);
    v = (Vector *)malloc((GMRESSteps + 2) * sizeof(Vector));
    if (v == NULL)
        AllocOK = False;
    else
        for (i = 1; i <= GMRESSteps + 1; i++) {
	    sprintf(vName, "v[%d]", i % 1000);
            V_Constr(&v[i], vName, Dim, Normal, True);
	}

    if (!AllocOK)
        LASError(LASMemAllocErr, "GMRESIter", NULL, NULL, NULL);
        
    if (LASResult() == LASOK) {
        bNorm = l2Norm_V(b);

        /* loop for 'MaxIter' GMRES cycles */
        Iter = 0;
        /* v[1] = r = b - A * x(i) */
        if (!IsZero(l1Norm_V(x) / Dim)) {
            if (Q_KerDefined(A))
	       OrthoRightKer_VQ(x, A);
            Asgn_VV(&v[1], Sub_VV(b, Mul_QV(A, x)));
        } else {
            Asgn_VV(&v[1], b);
	}

        while (!RTCResult(Iter, l2Norm_V(&v[1]), bNorm, GMRESIterId)
            && Iter < MaxIter) {
            if (PrecondProc != NULL)
	      {
	        (*PrecondProc)(A, &v[1], &v[1], OmegaPrecond);
	      }

            s[1] = l2Norm_V(&v[1]);
            MulAsgn_VS(&v[1], 1.0 / s[1]);
	    
            /* GMRES iteration */
            i = 0;
            while ((PrecondProc != NULL ? True : !RTCResult(Iter, fabs(s[i+1]),
                bNorm, GMRESIterId)) && i < GMRESSteps && Iter < MaxIter) {
 	        i++;
		Iter++;

                /* w = v[i+1] */
                if (PrecondProc != NULL) 
		    (*PrecondProc)(A, &v[i+1], Mul_QV(A, &v[i]), OmegaPrecond);
                else
                    Asgn_VV(&v[i+1], Mul_QV(A, &v[i]));

                /* modified Gram-Schmidt orthogonalization */
                for (k = 1; k <= i; k++) {
                    h[i][k] = Mul_VV(&v[i+1], &v[k]);
                    SubAsgn_VV(&v[i+1], Mul_SV(h[i][k], &v[k]));
                }

                h[i][i+1] = l2Norm_V(&v[i+1]);
                MulAsgn_VS(&v[i+1], 1.0 / h[i][i+1]);

                /* Q-R algorithm */
                for (k = 1; k < i; k++) {
                    h1 = c1[k] * h[i][k] + c2[k] * h[i][k+1];
                    h2 = - c2[k] * h[i][k] + c1[k] * h[i][k+1];
                    h[i][k] = h1;
                    h[i][k+1] = h2;
                }
                r = sqrt(h[i][i] * h[i][i] + h[i][i+1] * h[i][i+1]);

                c1[i] = h[i][i] / r;
                c2[i] = h[i][i+1] / r;
                h[i][i] = r;
                h[i][i+1] = 0.0;
                s[i+1] = - c2[i] * s[i];
                s[i] = c1[i] * s[i];

            }

            /* Solving of the system of equations : H y = s */
            for (j = i; j > 0; j--) {
                y[j] = s[j] / h[j][j];
		
                for (k = j - 1; k > 0; k--)
                    s[k] -= h[j][k] * y[j];
            }

            /* updating solution */
            for (j = i; j > 0; j--)
                AddAsgn_VV(x, Mul_SV(y[j], &v[j]));
            if (Q_KerDefined(A))
	       OrthoRightKer_VQ(x, A);
                
            /* computing new residual */
            Asgn_VV(&v[1], Sub_VV(b, Mul_QV(A, x)));
        }
    }

    /* release of vectors u, matrix H and vectors y, s, c1 and c2 */
    if (v != NULL) {
        for (i = 1; i <= GMRESSteps + 1; i++)
            V_Destr(&v[i]);
        free(v);
    }

    if (h != NULL) {
        for (i = 1; i <= GMRESSteps; i++)
	    if (h[i] != NULL)
                free(h[i]);
        free(h);
    }
    if (y != NULL)
        free(y);
    if (s != NULL)
        free(s);
    if (c1 != NULL)
        free(c1);
    if (c2 != NULL)
        free(c2);

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

    return(x);
}

Vector *BiCGIter(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, z, z_;

    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(&p_, "p_", Dim, Normal, True);
    V_Constr(&q, "q", Dim, Normal, True);
    if (PrecondProc != NULL || Q_KerDefined(A)) {
        V_Constr(&z, "z", Dim, Normal, True);
        V_Constr(&z_, "z_", 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 BiCG */
            Asgn_VV(&r_, &r);
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, BiCGIterId)
                && Iter < MaxIter) {
                Iter++;
		if (PrecondProc != NULL)
                    (*PrecondProc)(A, &z, &r, OmegaPrecond);
		else
		    Asgn_VV(&z, &r);
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&z, A);
		if (PrecondProc != NULL)
                    (*PrecondProc)(Transp_Q(A), &z_, &r_, OmegaPrecond);
		else
		    Asgn_VV(&z_, &r_);
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&z_, Transp_Q(A));
                Rho = Mul_VV(&z, &r_);
                if (IsZero(Rho)){
                    LASError(LASBreakdownErr, "BiCGIter", "Rho", NULL, NULL);
                    break;
                }
                if (Iter == 1) {
                    Asgn_VV(&p, &z);
                    Asgn_VV(&p_, &z_);
                } else {
                    Beta = Rho / RhoOld;
                    Asgn_VV(&p, Add_VV(&z, Mul_SV(Beta, &p)));
                    Asgn_VV(&p_, Add_VV(&z_, Mul_SV(Beta, &p_)));
                }
                Asgn_VV(&q, Mul_QV(A, &p));
                Alpha = Rho / Mul_VV(&p_, &q);
                AddAsgn_VV(x, Mul_SV(Alpha, &p));
                SubAsgn_VV(&r, Mul_SV(Alpha, &q));
                SubAsgn_VV(&r_, Mul_SV(Alpha, Mul_QV(Transp_Q(A), &p_)));
                RhoOld = Rho;
            }
        } else {
            /* plain BiCG (z = r, z_ = r_) */
            Asgn_VV(&r_, &r);
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, BiCGIterId)
                && Iter < MaxIter) {
                Iter++;
                Rho = Mul_VV(&r, &r_);
                if (IsZero(Rho)) {
                    LASError(LASBreakdownErr, "BiCGIter", "Rho", NULL, NULL);
                    break;
                }
                if (Iter == 1) {
                    Asgn_VV(&p, &r);
                    Asgn_VV(&p_, &r_);
                } else {
                    Beta = Rho / RhoOld;
                    Asgn_VV(&p, Add_VV(&r, Mul_SV(Beta, &p)));
                    Asgn_VV(&p_, Add_VV(&r_, Mul_SV(Beta, &p_)));
                }
                Asgn_VV(&q, Mul_QV(A, &p));
                Alpha = Rho / Mul_VV(&p_, &q);
                AddAsgn_VV(x, Mul_SV(Alpha, &p));
                SubAsgn_VV(&r, Mul_SV(Alpha, &q));
                SubAsgn_VV(&r_, Mul_SV(Alpha, Mul_QV(Transp_Q(A), &p_)));
                RhoOld = Rho;
            }
        }
        if (Q_KerDefined(A))
	   OrthoRightKer_VQ(x, A);
    }

    V_Destr(&r);
    V_Destr(&r_);
    V_Destr(&p);
    V_Destr(&p_);
    V_Destr(&q);
    if (PrecondProc != NULL || Q_KerDefined(A)) {
        V_Destr(&z);
        V_Destr(&z_);
    }

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

    return(x);
}

Vector *QMRIter(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 Rho, RhoNew, Xi, Gamma, GammaOld, Eta, Delta, Beta, Epsilon = 0.0,
        Theta, ThetaOld = 0.0;
    double bNorm;
    size_t Dim;
    Vector r, p, p_, q, y, y_, v, v_, w, w_, z, z_, s, d;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", Dim, Normal, True);
    V_Constr(&p, "p", Dim, Normal, True);
    V_Constr(&p_, "p_", Dim, Normal, True);
    V_Constr(&q, "q", Dim, Normal, True);
    V_Constr(&v, "v", Dim, Normal, True);
    V_Constr(&v_, "v_", Dim, Normal, True);
    V_Constr(&w, "w", Dim, Normal, True);
    V_Constr(&w_, "w_", Dim, Normal, True);
    V_Constr(&s, "s", Dim, Normal, True);
    V_Constr(&d, "d", Dim, Normal, True);
    if (PrecondProc != NULL || Q_KerDefined(A)) {
        V_Constr(&y, "y", Dim, Normal, True);
        V_Constr(&y_, "y_", Dim, Normal, True);
        V_Constr(&z, "z", Dim, Normal, True);
        V_Constr(&z_, "z_", 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 QMR (M1 ~ A, M2 = I) */
            Asgn_VV(&v_, &r);