operats.c

OpenFVM-v1.1 open source cfd code

/****************************************************************************/
/*                                operats.c                                 */
/****************************************************************************/
/*                                                                          */
/* basic OPERATionS for the types vector, matrix and qmatrix                */
/*                                                                          */
/* 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 <string.h>

#include "operats.h"
#include "errhandl.h"
#include "copyrght.h"

/* the following macro allows to optimize vector operations
   for different computer architecture */
#if !defined(__hppa) && !defined(sparc)
#define for_AllCmp for(Ind = 1; Ind <= Dim; Ind++)
#else
#define for_AllCmp for(Ind = Dim; Ind > 0; Ind--)
#endif   

Vector *Asgn_VV(Vector *V1, Vector *V2)
/* VRes = V1 = V2 */
{
    Vector *VRes;

    double Multipl;
    size_t Dim, Ind;
    Real *V1Cmp, *V2Cmp;

    V_Lock(V1);
    V_Lock(V2);
    
    if (LASResult() == LASOK) {
        if (V1->Instance == Normal && V1->Dim == V2->Dim) {
            V1Cmp = V1->Cmp;
            V2Cmp = V2->Cmp;
            if (V2->Instance == Normal) {
                Dim = V1->Dim;
                for_AllCmp
		    V1Cmp[Ind] = V2Cmp[Ind];
            } else {
                if (V2->OwnData) {
                    V1->Cmp = V2Cmp;
                    V2->Cmp = V1Cmp;
                } else {
                    if (IsOne(V2->Multipl)) {
                        Dim = V1->Dim;
                        for_AllCmp
			    V1Cmp[Ind] = V2Cmp[Ind];
                    } else {
                        Dim = V1->Dim;
                        Multipl = V2->Multipl;
                        for_AllCmp
			    V1Cmp[Ind] = Multipl * V2Cmp[Ind];
                    }
                }
            }
            VRes = V1;
        } else {
            if (V1->Instance != Normal)
                LASError(LASLValErr, "Asgn_VV", V_GetName(V1), V_GetName(V2), NULL);
            if (V1->Dim != V2->Dim)
                LASError(LASDimErr, "Asgn_VV", V_GetName(V1), V_GetName(V2), NULL);
            VRes = NULL;
        }
    } else {
        VRes = NULL;
    }
    
    V_Unlock(V1);
    V_Unlock(V2);

    return(VRes);
}

Vector *AddAsgn_VV(Vector *V1, Vector *V2)
/* VRes = V1 += V2 */
{
    Vector *VRes;

    double Multipl;
    size_t Dim, Ind;
    Real *V1Cmp, *V2Cmp;

    V_Lock(V1);
    V_Lock(V2);
    
    if (LASResult() == LASOK) {
        if (V1->Instance == Normal && V1->Dim == V2->Dim) {
            Dim = V1->Dim;
            V1Cmp = V1->Cmp;
            V2Cmp = V2->Cmp;
            if (IsOne(V2->Multipl)) {
	       for_AllCmp
		    V1Cmp[Ind] += V2Cmp[Ind];
            } else {
                Multipl = V2->Multipl;
                for_AllCmp
		    V1Cmp[Ind] += Multipl * V2Cmp[Ind];
            }
            VRes = V1;
        } else {
            if (V1->Instance != Normal)
                LASError(LASLValErr, "AddAsgn_VV", V_GetName(V1), V_GetName(V2), NULL);
            if (V1->Dim != V2->Dim) 
                LASError(LASDimErr, "AddAsgn_VV", V_GetName(V1), V_GetName(V2), NULL);
            VRes = NULL;
        }
    } else {
        VRes = NULL;
    }

    V_Unlock(V1);
    V_Unlock(V2);

    return(VRes);
}

Vector *SubAsgn_VV(Vector *V1, Vector *V2)
/* VRes = V1 -= V2 */
{
    Vector *VRes;

    double Multipl;
    size_t Dim, Ind;
    Real *V1Cmp, *V2Cmp;

    V_Lock(V1);
    V_Lock(V2);
    
    if (LASResult() == LASOK) {
        if (V1->Instance == Normal && V1->Dim == V2->Dim) {
            Dim = V1->Dim;
            V1Cmp = V1->Cmp;
            V2Cmp = V2->Cmp;
            if (IsOne(V2->Multipl)) {
	       for_AllCmp
		    V1Cmp[Ind] -= V2Cmp[Ind];
            } else {
                Multipl = V2->Multipl;
                for_AllCmp
		    V1Cmp[Ind] -= Multipl * V2Cmp[Ind];
            }
            VRes = V1;
        } else {
            if (V1->Instance != Normal)
                LASError(LASLValErr, "SubAsgn_VV", V_GetName(V1), V_GetName(V2), NULL);
            if (V1->Dim != V2->Dim) 
                LASError(LASDimErr, "SubAsgn_VV", V_GetName(V1), V_GetName(V2), NULL);
            VRes = NULL;
        }
    } else {
        VRes = NULL;
    }

    V_Unlock(V1);
    V_Unlock(V2);

    return(VRes);
}

Vector *MulAsgn_VS(Vector *V, double S)
/* VRes = V *= S */
{
    Vector *VRes;

    size_t Dim, Ind;
    Real *VCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
        if (V->Instance == Normal) {
            Dim = V->Dim;
            VCmp = V->Cmp;
            for_AllCmp
	        VCmp[Ind] *= S;
            VRes = V;
        } else {
            LASError(LASLValErr, "MulAsgn_VS", V_GetName(V), NULL, NULL);
            VRes = NULL;
        }
    } else {
        VRes = NULL;
    }

    V_Unlock(V);

    return(VRes);
}

Vector *Add_VV(Vector *V1, Vector *V2)
/* VRes = V1 + V2 */
{
    Vector *VRes;

    char *VResName;
    double Multipl1, Multipl2;
    size_t Dim, Ind;
    Real *V1Cmp, *V2Cmp, *VResCmp;

    V_Lock(V1);
    V_Lock(V2);
    
    if (LASResult() == LASOK) {
        if (V1->Dim == V2->Dim) {
            Dim = V1->Dim;
            VRes = (Vector *)malloc(sizeof(Vector));
	    VResName = (char *)malloc((strlen(V_GetName(V1)) + strlen(V_GetName(V2)) + 10)
		           * sizeof(char));
            if (VRes != NULL && VResName != NULL) {
	        sprintf(VResName, "%s + %s", V_GetName(V1), V_GetName(V2));
                V_Constr(VRes, VResName, Dim, Tempor, True);
                if (LASResult() == LASOK) {
                    V1Cmp = V1->Cmp;
                    V2Cmp = V2->Cmp;
                    VResCmp = VRes->Cmp;
                    if (IsOne(V1->Multipl)) {
                        if (IsOne(V2->Multipl)) {
			   for_AllCmp
			        VResCmp[Ind] = V1Cmp[Ind] + V2Cmp[Ind];
                        } else {
                            Multipl2 = V2->Multipl;
                            for_AllCmp
			        VResCmp[Ind] = V1Cmp[Ind] + Multipl2 * V2Cmp[Ind];
                        }
                    } else {
                        Multipl1 = V1->Multipl;
                        if (IsOne(V2->Multipl)) {
			   for_AllCmp
			       VResCmp[Ind] = Multipl1 * V1Cmp[Ind] + V2Cmp[Ind];
                        } else {
                            Multipl2 = V2->Multipl;
                            for_AllCmp
			        VResCmp[Ind] = Multipl1 * V1Cmp[Ind] + Multipl2 * V2Cmp[Ind];
                        }
                    }
                }
            } else {
                LASError(LASMemAllocErr, "Add_VV", V_GetName(V1), V_GetName(V2), NULL);
		if (VRes != NULL)
		    free(VRes);
            }
	    
	    if (VResName != NULL)
	        free(VResName);
        } else {
            LASError(LASDimErr, "Add_VV", V_GetName(V1), V_GetName(V2), NULL);
            VRes = NULL;
        }
    } else {
        VRes = NULL;
    }

    V_Unlock(V1);
    V_Unlock(V2);

    return(VRes);
}

QMatrix *Add_QQ(QMatrix *Q1, QMatrix *Q2)
/* QRes = Q1 + Q2, this operation is limited to matrices, which are
   derived from the same matrix */
{
    QMatrix *QRes;

    char *QResName;

    Q_Lock(Q1);
    Q_Lock(Q2);
    
    if (LASResult() == LASOK) {
        if (Q1->Dim == Q2->Dim) {
            if (Q1->Symmetry == Q2->Symmetry &&
                Q1->ElOrder == Q2->ElOrder &&
                Q1->Len == Q2->Len &&
                Q1->El == Q2->El) {
                QRes = (QMatrix *)malloc(sizeof(QMatrix));
	        QResName = (char *)malloc((strlen(Q_GetName(Q1)) + strlen(Q_GetName(Q1)) + 10)
		           * sizeof(char));
                if (QRes != NULL && QResName != NULL) {
	            sprintf(QResName, "%s + %s", Q_GetName(Q1), Q_GetName(Q2));
                    Q_Constr(QRes, QResName, Q1->Dim, Q1->Symmetry, Q1->ElOrder, Tempor, False);
                    if (LASResult() == LASOK) {
                        QRes->MultiplD = Q1->MultiplD + Q2->MultiplD;
                        QRes->MultiplU = Q1->MultiplU + Q2->MultiplU;
                        QRes->MultiplL = Q1->MultiplL + Q2->MultiplL;
                        QRes->Len = Q1->Len;
                        QRes->El = Q1->El;
                        QRes->ElSorted = Q1->ElSorted;
                        QRes->DiagElAlloc = Q1->DiagElAlloc;
                        QRes->DiagEl = Q1->DiagEl;
                        QRes->ZeroInDiag = Q1->ZeroInDiag;
                        QRes->InvDiagEl = Q1->InvDiagEl;
                        QRes->UnitRightKer = Q1->UnitRightKer;
                        QRes->RightKerCmp = Q1->RightKerCmp;
                        QRes->UnitLeftKer = Q1->UnitLeftKer;
                        QRes->LeftKerCmp = Q1->LeftKerCmp;
                        QRes->ILUExists = Q1->ILUExists;
                        QRes->ILU = Q1->ILU;
                    }
                } else {
                    LASError(LASMemAllocErr, "Add_QQ", Q_GetName(Q1), Q_GetName(Q2), NULL);
	            if (QRes != NULL)
	                free(QRes);
                }
	    
	        if (QResName != NULL)
	            free(QResName);
            } else {
                LASError(LASMatrCombErr, "Add_QQ", Q_GetName(Q1), Q_GetName(Q2), NULL);
                QRes = NULL;
            }
        } else {
            LASError(LASDimErr, "Add_QQ", Q_GetName(Q1), Q_GetName(Q2), NULL);
            QRes = NULL;
        }
    } else {
        QRes = NULL;
    }

    Q_Unlock(Q1);
    Q_Unlock(Q2);

    return(QRes);
}

Vector *Sub_VV(Vector *V1, Vector *V2)
/* VRes = V1 - V2 */
{
    Vector *VRes;

    char *VResName;
    double Multipl1, Multipl2;
    size_t Dim, Ind;
    Real *V1Cmp, *V2Cmp, *VResCmp;

    V_Lock(V1);
    V_Lock(V2);
    
    if (LASResult() == LASOK) {
        if (V1->Dim == V2->Dim) {
            Dim = V1->Dim;
            VRes = (Vector *)malloc(sizeof(Vector));
	    VResName = (char *)malloc((strlen(V_GetName(V1)) + strlen(V_GetName(V2)) + 10)
		           * sizeof(char));
            if (VRes != NULL && VResName != NULL) {
	        sprintf(VResName, "%s - %s", V_GetName(V1), V_GetName(V2));
                V_Constr(VRes, VResName, Dim, Tempor, True);
                if (LASResult() == LASOK) {
                    V1Cmp = V1->Cmp;
                    V2Cmp = V2->Cmp;
                    VResCmp = VRes->Cmp;
                    if (IsOne(V1->Multipl)) {
                        if (IsOne(V2->Multipl)) {
			   for_AllCmp
			        VResCmp[Ind] = V1Cmp[Ind] - V2Cmp[Ind];
                        } else {
                            Multipl2 = V2->Multipl;
                            for_AllCmp
			        VResCmp[Ind] = V1Cmp[Ind] - Multipl2 * V2Cmp[Ind];
                        }
                    } else {
                        Multipl1 = V1->Multipl;
                        if (IsOne(V2->Multipl)) {
			   for_AllCmp
			        VResCmp[Ind] = Multipl1 * V1Cmp[Ind] - V2Cmp[Ind];
                        } else {
                            Multipl2 = V2->Multipl;
                            for_AllCmp
			        VResCmp[Ind] = Multipl1 * V1Cmp[Ind] - Multipl2 * V2Cmp[Ind];
                        }
                    }
                }
            } else {
                LASError(LASMemAllocErr, "Sub_VV", V_GetName(V1), V_GetName(V2), NULL);
		if (VRes != NULL)
		    free(VRes);
            }
	    
	    if (VResName != NULL)
	        free(VResName);
        } else {
            LASError(LASDimErr, "Sub_VV", V_GetName(V1), V_GetName(V2), NULL);
            VRes = NULL;
        }
    } else {
        VRes = NULL;
    }

    V_Unlock(V1);
    V_Unlock(V2);

    return(VRes);
}

QMatrix *Sub_QQ(QMatrix *Q1, QMatrix *Q2)
/* QRes = Q1 - Q2, this operation ist limited to matricies, which are
   derived from the same matrix */
{
    QMatrix *QRes;

    char *QResName;

    Q_Lock(Q1);
    Q_Lock(Q2);
    
    if (LASResult() == LASOK) {
        if (Q1->Dim == Q2->Dim) {
            if (Q1->Symmetry == Q2->Symmetry &&
                Q1->ElOrder == Q2->ElOrder &&
                Q1->Len == Q2->Len &&