operats.c

OpenFVM-v1.1 open source cfd code

                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, "Sub_QQ", Q_GetName(Q1), Q_GetName(Q2), NULL);
	            if (QRes != NULL)
	                free(QRes);
                }
	    
	        if (QResName != NULL)
	            free(QResName);
            } else {
                LASError(LASMatrCombErr, "Sub_QQ", Q_GetName(Q1), Q_GetName(Q2), NULL);
                QRes = NULL;
            }
        } else {
            LASError(LASDimErr, "Sub_QQ", Q_GetName(Q1), Q_GetName(Q2), NULL);
            QRes = NULL;
        }
    } else {
        QRes = NULL;
    }

    Q_Unlock(Q1);
    Q_Unlock(Q2);

    return(QRes);
}

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

    char *VResName;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
        VRes = (Vector *)malloc(sizeof(Vector));
	VResName = (char *)malloc((strlen(V_GetName(V)) + 20) * sizeof(char));
        if (VRes != NULL && VResName != NULL) {
	    sprintf(VResName, "%12.5e * (%s)", S, V_GetName(V));
            V_Constr(VRes, VResName, V->Dim, Tempor, False);
            if (LASResult() == LASOK) {
                VRes->Multipl = S * V->Multipl;
                if (V->Instance == Tempor && V->OwnData) {
                    V->OwnData = False;
                    VRes->OwnData = True;
                }
                VRes->Cmp = V->Cmp;
            }
        } else {
            LASError(LASMemAllocErr, "Mul_SV", V_GetName(V), NULL, NULL);
            if (VRes != NULL)
	        free(VRes);
        }
	    
	if (VResName != NULL)
	    free(VResName);
    } else {
        VRes = NULL;
    }

    V_Unlock(V);

    return(VRes);
}

Matrix *Mul_SM(double S, Matrix *M)
/* MRes = S * M */
{
    Matrix *MRes;

    char *MResName;

    M_Lock(M);

    if (LASResult() == LASOK) {
        MRes = (Matrix *)malloc(sizeof(Matrix));
	MResName = (char *)malloc((strlen(M_GetName(M)) + 20) * sizeof(char));
        if (MRes != NULL && MResName != NULL) {
	    sprintf(MResName, "%12.5e * (%s)", S, M_GetName(M));
            M_Constr(MRes, MResName,  M->RowDim, M->ClmDim, M->ElOrder, Tempor, False);
            if (LASResult() == LASOK) {
                MRes->Multipl = S * M->Multipl;
                MRes->Len = M->Len;
                MRes->El = M->El;
                MRes->ElSorted = M->ElSorted;
            }
        } else {
            LASError(LASMemAllocErr, "Mul_SM", M_GetName(M), NULL, NULL);
	    if (MRes != NULL)
	        free(MRes);
        }
	    
	if (MResName != NULL)
	    free(MResName);
    } else {
        MRes = NULL;
    }

    M_Unlock(M);

    return(MRes);
}

QMatrix *Mul_SQ(double S, QMatrix *Q)
/* QRes = S * Q */
{
    QMatrix *QRes;
    
    char *QResName;

    Q_Lock(Q);
    
    if (LASResult() == LASOK) {
        QRes = (QMatrix *)malloc(sizeof(QMatrix));
	QResName = (char *)malloc((strlen(Q_GetName(Q)) + 20) * sizeof(char));
        if (QRes != NULL && QResName != NULL) {
	    sprintf(QResName, "%12.5e * (%s)", S, Q_GetName(Q));
            Q_Constr(QRes, QResName, Q->Dim, Q->Symmetry, Q->ElOrder, Tempor, False);
            if (LASResult() == LASOK) {
                if (Q->Instance == Tempor && Q->OwnData) {
                    Q->OwnData = False;
                    QRes->OwnData = True;
                }
                QRes->MultiplD = S * Q->MultiplD;
                QRes->MultiplU = S * Q->MultiplU;
                QRes->MultiplL = S * Q->MultiplL;
                QRes->Len = Q->Len;
                QRes->El = Q->El;
                QRes->ElSorted = Q->ElSorted;
                QRes->DiagElAlloc = Q->DiagElAlloc;
                QRes->DiagEl = Q->DiagEl;
                QRes->ZeroInDiag = Q->ZeroInDiag;
                QRes->InvDiagEl = Q->InvDiagEl;
                QRes->UnitRightKer = Q->UnitRightKer;
                QRes->RightKerCmp = Q->RightKerCmp;
                QRes->UnitLeftKer = Q->UnitLeftKer;
                QRes->LeftKerCmp = Q->LeftKerCmp;
                QRes->ILUExists = Q->ILUExists;
                QRes->ILU = Q->ILU;
            }
        } else {
            LASError(LASMemAllocErr, "Mul_SQ", Q_GetName(Q), NULL, NULL);
	    if (QRes != NULL)
	        free(QRes);
        }
	    
	if (QResName != NULL)
	    free(QResName);
    } else {
        QRes = NULL;
    }

    Q_Unlock(Q);

    return(QRes);
}

double Mul_VV(Vector *V1, Vector *V2)
/* S = V1 * V2 */
{
    double SRes;

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

    V_Lock(V1);
    V_Lock(V2);
    
    if (LASResult() == LASOK) {
        if (V1->Dim == V2->Dim) {
            Dim = V1->Dim;
            V1Cmp = V1->Cmp;
            V2Cmp = V2->Cmp;
            SRes = 0.0;
            for_AllCmp
	        SRes += V1Cmp[Ind] * V2Cmp[Ind];
	    SRes *= V1->Multipl * V2->Multipl;
        } else {
            LASError(LASDimErr, "Mul_VV", V_GetName(V1),  V_GetName(V2), NULL);
            SRes = 1.0;
        }
    } else {
        SRes = 1.0;
    }

    V_Unlock(V1);
    V_Unlock(V2);

    return(SRes);
}

Vector *Mul_MV(Matrix *M, Vector *V)
/* VRes = M * V */
{
    Vector *VRes;

    char *VResName;
    double MultiplMV;
    double Sum, Cmp;
    size_t RowDim, ClmDim, Row, Clm, Len, ElCount;
    size_t *MLen;
    Boolean MultiplMVIsOne;
    ElType **MEl, *PtrEl;
    Real *VCmp, *VResCmp;

    M_Lock(M);
    V_Lock(V);
    
    if (LASResult() == LASOK) {
	if (M->ClmDim == V->Dim) {
            RowDim = M->RowDim;
            ClmDim = M->ClmDim;
	    VRes = (Vector *)malloc(sizeof(Vector));
	    VResName = (char *)malloc((strlen(M_GetName(M)) + strlen(V_GetName(V)) + 10)
		           * sizeof(char));
            if (VRes != NULL && VResName != NULL) {
	        sprintf(VResName, "(%s) * (%s)", M_GetName(M), V_GetName(V));
                V_Constr(VRes, VResName, RowDim, Tempor, True);
		if (LASResult() == LASOK) {
                    /* assignment of auxiliary lokal variables */
                    MLen = M->Len;
                    MEl = M->El;
                    VCmp = V->Cmp;
                    VResCmp = VRes->Cmp;

                    /* initialisation of the vector VRes */
		    if (M->ElOrder == Clmws)
                        V_SetAllCmp(VRes, 0.0);

                    /* analysis of multipliers of the matrix M and the vector V */
                    MultiplMV = M->Multipl * V->Multipl;
                    if (IsOne(MultiplMV)) {
                        MultiplMVIsOne = True;
                    } else {
                        MultiplMVIsOne = False;
                    }

                    /* multiplication of matrix elements by vector components */
		    if (M->ElOrder == Rowws) {
			for (Row = 1; Row <= RowDim; Row++) {
                            Len = MLen[Row];
                            PtrEl = MEl[Row];
                            Sum = 0.0;
			    for (ElCount = Len; ElCount > 0; ElCount--) {
                                Sum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                PtrEl++;
			    }
                            if (MultiplMVIsOne)
                                VResCmp[Row] = Sum;
                            else
                                VResCmp[Row] = MultiplMV * Sum;
			}
		    }
		    if (M->ElOrder == Clmws) {
			for (Clm = 1; Clm <= ClmDim; Clm++) {
                            Len = MLen[Clm];
                            PtrEl = MEl[Clm];
                            if (MultiplMVIsOne)
   			        Cmp = VCmp[Clm];
                            else
        		        Cmp = MultiplMV * VCmp[Clm];
			    for (ElCount = Len; ElCount > 0; ElCount--) {
				VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * Cmp;
                                PtrEl++;
			    }
			}
		    }
		}
	    } else {
		LASError(LASMemAllocErr, "Mul_MV", M_GetName(M), V_GetName(V), NULL);
		if (VRes != NULL)
		    free(VRes);
            }
	    
	    if (VResName != NULL)
	        free(VResName);
	} else {
	    LASError(LASDimErr, "Mul_MV", M_GetName(M), V_GetName(V), NULL);
            VRes = NULL;
	}
    } else {
        VRes = NULL;
    }

    M_Unlock(M);
    V_Unlock(V);
 
    return(VRes);
}

Vector *Mul_QV(QMatrix *Q, Vector *V)
/* VRes = Q * V */
{
    Vector *VRes;

    char *VResName;
    double MultiplDV, MultiplUV, MultiplLV;
    double Sum, PartSum, Cmp, PartCmp;
    size_t Dim, Row, Clm, RoC, Len, ElCount;
    size_t *QLen;
    Boolean MultiplDVIsZero, MultiplUVIsZero, MultiplLVIsZero;
    Boolean MultiplDVIsOne, MultiplUVIsOne, MultiplLVIsOne;
    Boolean MultiplDULVEquals;
    ElType **QEl, **QDiagEl, *PtrEl;
    Real *VCmp, *VResCmp;

    Q_Lock(Q);
    V_Lock(V);
    
    if (LASResult() == LASOK) {
	if (Q->Dim == V->Dim) {
	    Dim = V->Dim;
	    VRes = (Vector *)malloc(sizeof(Vector));
	    VResName = (char *)malloc((strlen(Q_GetName(Q)) + strlen(V_GetName(V)) + 10)
		           * sizeof(char));
            if (VRes != NULL && VResName != NULL) {
	        sprintf(VResName, "(%s) * (%s)", Q_GetName(Q), V_GetName(V));
                V_Constr(VRes, VResName, Dim, Tempor, True);
		
                /* sort of elements and allocation of diagonal elements
                   of the matrix Q */
                Q_SortEl(Q);
                Q_AllocInvDiagEl(Q);

		if (LASResult() == LASOK && Q->ElSorted) {
                    /* assignment of auxiliary lokal variables */
                    QLen = Q->Len;
                    QEl = Q->El;
                    QDiagEl = Q->DiagEl;
		    VCmp = V->Cmp;
		    VResCmp = VRes->Cmp;

                    /* initialisation of the vector VRes */
		    if (Q->Symmetry || Q->ElOrder == Clmws)
                        V_SetAllCmp(VRes, 0.0);

                    /* analysis of multipliers of the lower, diagonal and upper part
                       of the matrix Q and of the vector V */
                    MultiplDV = Q->MultiplD * V->Multipl;
                    MultiplUV = Q->MultiplU * V->Multipl;
                    MultiplLV = Q->MultiplL * V->Multipl;
                    MultiplDVIsZero = IsZero(MultiplDV);
                    MultiplUVIsZero = IsZero(MultiplUV);
                    MultiplLVIsZero = IsZero(MultiplLV);
                    MultiplDVIsOne = IsOne(MultiplDV);
                    MultiplUVIsOne = IsOne(MultiplUV);
                    MultiplLVIsOne = IsOne(MultiplLV);
                    if (!IsZero(MultiplDV) && IsOne(MultiplUV / MultiplDV)
			&& IsOne(MultiplLV / MultiplDV)) {
                        MultiplDULVEquals = True;
                    } else {
                        MultiplDULVEquals = False;
                    }

                    /* multiplication of the lower, diagonal and upper part
                       of the matrix Q by the vector V */
                    if (Q->Symmetry) {
                        if (MultiplDULVEquals) {
                            for (RoC = 1; RoC <= Dim; RoC++) {
                                Len = QLen[RoC];
                                PtrEl = QEl[RoC];
   	                        Sum = 0.0;
                                if (!MultiplDVIsOne)
        		            Cmp = MultiplDV * VCmp[RoC];
                                else
   			            Cmp = VCmp[RoC];
			        for (ElCount = Len; ElCount > 0; ElCount--) {
                                    if ((*PtrEl).Pos != RoC) {
			                VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * Cmp;
			                Sum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                    }
                                    PtrEl++;
			        }
	                        Sum += (*QDiagEl[RoC]).Val * VCmp[RoC];
                                if (MultiplDVIsOne)
			            VResCmp[RoC] += Sum;
                                else
			            VResCmp[RoC] += MultiplDV * Sum;
      			    }
                        } else {
                            for (RoC = 1; RoC <= Dim; RoC++) {
                                Len = QLen[RoC];
   	                        Sum = 0.0;
                                if ((!MultiplUVIsZero && Q->ElOrder == Rowws)
                                    || (!MultiplLVIsZero && Q->ElOrder == Clmws)) {
                                    PtrEl = QEl[RoC];
    			            for (ElCount = Len; ElCount > 0; ElCount--) {
                                        if ((*PtrEl).Pos != RoC) {
                                            Sum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                        }
                                        PtrEl++;
	                            }
                                    if (!MultiplUVIsZero && !MultiplUVIsOne)
   	                                Sum = MultiplUV * Sum;
                                    if (!MultiplLVIsZero && !MultiplLVIsOne)
   	                                Sum = MultiplLV * Sum;
                                }