itersolv.c

OpenFVM-v1.1 open source cfd code

/****************************************************************************/
/*                                itersolv.c                                */
/****************************************************************************/
/*                                                                          */
/* ITERative SOLVers for systems of linear equations                        */
/*                                                                          */
/* Copyright (C) 1992-1996 Tomas Skalicky. All rights reserved.             */
/*                                                                          */
/****************************************************************************/
/*                                                                          */
/*        ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS          */
/*              OF THE COPYRIGHT NOTICE (SEE FILE COPYRGHT.H)               */
/*                                                                          */
/****************************************************************************/

#include <stddef.h>
#include <math.h>

#include "itersolv.h"
#include "elcmp.h"
#include "errhandl.h"
#include "operats.h"
#include "rtc.h"
#include "copyrght.h"

/* number of GMRES steps before restart */
static int GMRESSteps = 10;

Vector *JacobiIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
            PrecondProcType Dummy, double Omega)
{
    int Iter;
    double bNorm;
    size_t Dim;
    Vector r;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", 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);
	}
        while (!RTCResult(Iter, l2Norm_V(&r), bNorm, JacobiIterId)
            && Iter < MaxIter) {
            Iter++;
            /* x(i+1) = x(i) + Omega * D^(-1) * r */
            AddAsgn_VV(x, Mul_SV(Omega, MulInv_QV(Diag_Q(A), &r)));
            if (Q_KerDefined(A))
	       OrthoRightKer_VQ(x, A);
            /* r = b - A * x(i) */
            if (Iter < MaxIter)
                Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
        }
    }

    V_Destr(&r);

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

    return(x);
}

Vector *SORForwIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
            PrecondProcType Dummy, double Omega)
{
    int Iter;
    double bNorm;
    size_t Dim;
    Vector r;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", 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);
	}
        while (!RTCResult(Iter, l2Norm_V(&r), bNorm, SORForwIterId)
            && Iter < MaxIter) {
            Iter++;
            /* x(i+1) = x(i) + (D / Omega + L)^(-1) r */
            AddAsgn_VV(x, MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Lower_Q(A)), &r));
            if (Q_KerDefined(A))
	       OrthoRightKer_VQ(x, A);
            /* r = b - A * x(i) */
            if (Iter < MaxIter)
                Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
        }
    }

    V_Destr(&r);

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

    return(x);
}

Vector *SORBackwIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
            PrecondProcType Dummy, double Omega)
{
    int Iter;
    double bNorm;
    size_t Dim;
    Vector r;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", 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);
	}
        while (!RTCResult(Iter, l2Norm_V(&r), bNorm, SORBackwIterId)
            && Iter < MaxIter) {
            Iter++;
            /* x(i+1) = x(i) + (D / Omega + U)^(-1) r */
            AddAsgn_VV(x, MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Upper_Q(A)), &r));
            if (Q_KerDefined(A))
	       OrthoRightKer_VQ(x, A);
            /* r = b - A * x(i) */
            if (Iter < MaxIter)
                Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
        }
    }

    V_Destr(&r);

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

    return(x);
}

Vector *SSORIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
            PrecondProcType Dummy, double Omega)
{
    int Iter;
    double bNorm;
    size_t Dim;
    Vector r;

    Q_Lock(A);
    V_Lock(x);
    V_Lock(b);
    
    Dim = Q_GetDim(A);
    V_Constr(&r, "r", 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);
	}
        while (!RTCResult(Iter, l2Norm_V(&r), bNorm, SSORIterId)
            && Iter < MaxIter) {
            Iter++;
            /* x(i+1) = x(i) + (2 - Omega) * (Diag(A) / Omega + Upper(A))^(-1)
                   * Diag(A) / Omega * (Diag(A) / Omega + Lower(A))^(-1) r */
            AddAsgn_VV(x, Mul_SV((2.0 - Omega) / Omega,
                MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Upper_Q(A)),
                Mul_QV(Diag_Q(A),
                MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Lower_Q(A)), &r)))));
            if (Q_KerDefined(A))
	       OrthoRightKer_VQ(x, A);
            /* r = b - A * x(i) */
            if (Iter < MaxIter)
                Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
        }
    }

    V_Destr(&r);

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

    return(x);
}

Vector *ChebyshevIter(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
     *
     *  The original algorithm doesn't seem to work correctly.
     *  The three term recursion was therefore in the curent version
     *  adopted after
     *
     *  W. Hackbusch: 
     *  Iterative Solution of Large Sparse Systems of Equations,
     *  Springer-Verlag, Berlin, 1994
     *
     */
     
    int Iter;
    double MaxEigenval, MinEigenval;
    double c, d, Alpha, Beta;
    double T = 0.0, TOld = 0.0, TNew; /* values of Chebyshev polynomials */
    double bNorm;
    size_t Dim;
    Vector r, p, z;
    

    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);
    if (PrecondProc != NULL || Q_KerDefined(A))
        V_Constr(&z, "z", Dim, Normal, True);

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

        MaxEigenval = GetMaxEigenval(A, PrecondProc, OmegaPrecond);
        MinEigenval = GetMinEigenval(A, PrecondProc, OmegaPrecond);

        c = (MaxEigenval - MinEigenval) / 2.0;
        d = (MaxEigenval + MinEigenval) / 2.0;
	
        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 Chebyshev method */
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, ChebyshevIterId)
                && 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 (Iter == 1) {
                    TOld = 1.0;
	            T = d / c;
                    Alpha = 1.0 / d;
                    Asgn_VV(&p, Mul_SV(Alpha, &z));
                } else {
		    TNew = 2.0 * d / c * T - TOld;
		    TOld = T;
		    T = TNew;
                    Alpha = 2.0 / c * TOld / T;
                    Beta = 2.0 * d / c * TOld / T - 1.0;
                    Asgn_VV(&p, Add_VV(Mul_SV(Alpha, &z), Mul_SV(Beta, &p)));
                }
                AddAsgn_VV(x, &p);
                if (Iter < MaxIter)
                    Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
            }
        } else {
            /* plain Chebyshev method (z = r) */
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, ChebyshevIterId)
                && Iter < MaxIter) {
                Iter++;
                if (Iter == 1) {
                    TOld = 1.0;
	            T = d / c;
                    Alpha = 1.0 / d;
                    Asgn_VV(&p, Mul_SV(Alpha, &r));
                } else {
		    TNew = 2.0 * d / c * T - TOld;
		    TOld = T;
		    T = TNew;
                    Alpha = 2.0 / c * TOld / T;
                    Beta = 2.0 * d / c * TOld / T - 1.0;
                    Asgn_VV(&p, Add_VV(Mul_SV(Alpha, &r), Mul_SV(Beta, &p)));
                }
                AddAsgn_VV(x, &p);
                if (Iter < MaxIter)
                    Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
            }
        }
        if (Q_KerDefined(A))
	   OrthoRightKer_VQ(x, A);
    }

    V_Destr(&r);
    V_Destr(&p);
    if (PrecondProc != NULL || Q_KerDefined(A))
        V_Destr(&z);

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

    return(x);
}

Vector *CGIter(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, p, q, z;

    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);
    if (PrecondProc != NULL || Q_KerDefined(A))
        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 CG */
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGIterId)
                && Iter < MaxIter) {
                Iter++;
		if (PrecondProc != NULL)
                    (*PrecondProc)(A, &z, &r, OmegaPrecond);
		else
		    Asgn_VV(&z, &r);
		if (Q_KerDefined(A))
		   OrthoRightKer_VQ(&z, A);
                Rho = Mul_VV(&r, &z);
                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 / Mul_VV(&p, &q);
                AddAsgn_VV(x, Mul_SV(Alpha, &p));
                SubAsgn_VV(&r, Mul_SV(Alpha, &q));
                RhoOld = Rho;
            }
        } else {
            /* plain CG (z = r) */
            while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGIterId)
                && Iter < MaxIter) {
                Iter++;
                Rho = pow(l2Norm_V(&r), 2.0);
                if (Iter == 1) {
                    Asgn_VV(&p, &r);
                } else {
                    Beta = Rho / RhoOld;
                    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));
                RhoOld = Rho;
            }
        }
        if (Q_KerDefined(A))
	   OrthoRightKer_VQ(x, A);
    }
    
    V_Destr(&r);
    V_Destr(&p);
    V_Destr(&q);
    if (PrecondProc != NULL || Q_KerDefined(A))
        V_Destr(&z);

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

    return(x);
}

Vector *CGNIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
            PrecondProcType PrecondProc, double OmegaPrecond)
{
    int Iter;
    double Alpha, Beta, Rho, RhoOld = 0.0;
    double bNorm;
    size_t Dim;
    Vector r, p, q, s, z;