operats.c

OpenFVM-v1.1 open source cfd code

                                if (!MultiplDVIsZero)
                                    if (MultiplDVIsOne)
			                Sum += (*QDiagEl[RoC]).Val * VCmp[RoC];
                                    else
			                Sum += MultiplDV * (*QDiagEl[RoC]).Val * VCmp[RoC];
			        VResCmp[RoC] += Sum;
                                if ((!MultiplUVIsZero && Q->ElOrder == Clmws)
                                    || (!MultiplLVIsZero && Q->ElOrder == Rowws)) {
  	                            Cmp = VCmp[RoC];
                                    if (!MultiplUVIsZero && !MultiplUVIsOne)
         		                Cmp *= MultiplUV;
                                    if (!MultiplLVIsZero && !MultiplLVIsOne)
             		                Cmp *= MultiplLV;
                                    PtrEl = QEl[RoC];
    			            for (ElCount = Len; ElCount > 0; ElCount--) {
                                        if ((*PtrEl).Pos != RoC)
 			                    VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * Cmp;
                                        PtrEl++;
	                            }
                                }
  			    }
                        }
                    }
                    if (!Q->Symmetry && Q->ElOrder == Rowws) {
                        if (MultiplDULVEquals) {
                            for (Row = 1; Row <= Dim; Row++) {
                                Len = QLen[Row];
                                PtrEl = QEl[Row];
		    	        Sum = 0.0;
			        for (ElCount = Len; ElCount > 0; ElCount--) {
			            Sum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                    PtrEl++;
			        }
                                if (MultiplDVIsOne)
			            VResCmp[Row] = Sum;
                                else
			            VResCmp[Row] = MultiplDV * Sum;
  			    }
  			} else {
                            for (Row = 1; Row <= Dim; Row++) {
                                Len = QLen[Row];
                                Sum = 0.0;
                                if (!MultiplLVIsZero) {
                                    PtrEl = QEl[Row];
		    	            PartSum = 0.0;
  			            for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos < Row;
                                        ElCount--) {
			                PartSum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                        PtrEl++;
			            }
                                    if (MultiplLVIsOne)
			                Sum += PartSum;
                                    else
			                Sum += MultiplLV * PartSum;
                                }
                                if (!MultiplDVIsZero) {
                                    if (MultiplDVIsOne)
			                Sum += (*QDiagEl[Row]).Val * VCmp[Row];
                                    else
			                Sum += MultiplDV * (*QDiagEl[Row]).Val * VCmp[Row];
                                }
                                if (!MultiplUVIsZero) {
                                    PtrEl = QEl[Row] + Len - 1;
		    	            PartSum = 0.0;
  			            for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos > Row;
                                        ElCount--) {
			                PartSum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                        PtrEl--;
			            }
                                    if (MultiplUVIsOne)
			                Sum += PartSum;
                                    else
			                Sum += MultiplUV * PartSum;
                                }
			        VResCmp[Row] = Sum;
  			    }
                        }
                    }
                    if (!Q->Symmetry && Q->ElOrder == Clmws) {
                        if (MultiplDULVEquals) {
 			    for (Clm = 1; Clm <= Dim; Clm++) {
                                Len = QLen[Clm];
                                PtrEl = QEl[Clm];
                                if (MultiplDVIsOne)
   			            Cmp = VCmp[Clm];
                                else
        		            Cmp = MultiplDV * VCmp[Clm];
			        for (ElCount = Len; ElCount > 0; ElCount--) {
				    VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * Cmp;
                                    PtrEl++;
			        }
			    }
  			} else {
                            for (Clm = 1; Clm <= Dim; Clm++) {
                                Len = QLen[Clm];
		                Cmp = VCmp[Clm];
                                if (!MultiplUVIsZero) {
                                    PtrEl = QEl[Clm];
                                    if (MultiplUVIsOne)
                                        PartCmp = Cmp;
                                    else
                                        PartCmp = MultiplUV * Cmp;
  			            for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos < Clm;
                                        ElCount--) {
 				        VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * PartCmp;
                                        PtrEl++;
			            }
                                }
                                if (!MultiplDVIsZero) {
                                    if (MultiplDVIsOne)
 				        VResCmp[Clm] += (*QDiagEl[Clm]).Val * Cmp;
                                    else
 				        VResCmp[Clm] += MultiplDV * (*QDiagEl[Clm]).Val * Cmp;
                                }
                                if (!MultiplLVIsZero) {
                                    PtrEl = QEl[Clm] + Len - 1;
                                    if (MultiplLVIsOne)
                                        PartCmp = Cmp;
                                    else
                                        PartCmp = MultiplLV * Cmp;
  			            for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos > Clm;
                                        ElCount--) {
 				        VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * Cmp;
                                        PtrEl--;
			            }
                                }
  			    }
                        }
                    }
                } else {
                    if (LASResult() == LASOK && !(*Q->ElSorted))
                        LASError(LASElNotSortedErr, "Mul_QV", Q_GetName(Q), V_GetName(V), NULL);
                }
	    } else {
		LASError(LASMemAllocErr, "Mul_QV", Q_GetName(Q), V_GetName(V), NULL);
		if (VRes != NULL)
		    free(VRes);
            }
	    
	    if (VResName != NULL)
	        free(VResName);
	} else {
	    LASError(LASDimErr, "Mul_QV", Q_GetName(Q), V_GetName(V), NULL);
	    VRes = NULL;
	}
    } else {
        VRes = NULL;
    }

    Q_Unlock(Q);
    V_Unlock(V);

    return(VRes);
}

Vector *MulInv_QV(QMatrix *Q, Vector *V)
/* VRes = Q^(-1) * V, this operation is limited to diagonal or triangular
   matrices */
{
    Vector *VRes;

    char *VResName;
    double MultiplD, MultiplU, MultiplL, MultiplV, MultiplDV;
    double Sum, Cmp;
    size_t Dim, Row, Clm, Ind, Len, ElCount;
    size_t *QLen;
    Boolean MultiplDIsZero, MultiplUIsZero, MultiplLIsZero;
    Boolean MultiplDIsOne, MultiplUIsOne, MultiplLIsOne, MultiplVIsOne,
        MultiplDVIsOne;
    ElType **QEl, *PtrEl;
    Real *QInvDiagEl;
    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)) + 20)
		           * sizeof(char));
            if (VRes != NULL && VResName != NULL) {
	        sprintf(VResName, "(%s)^(-1) * (%s)", Q_GetName(Q), V_GetName(V));
                V_Constr(VRes, VResName, Dim, Tempor, True);

                /* sort of elements and allocation of the inverse of diagonal elements
                   of the matrix Q */
                Q_SortEl(Q);
                Q_AllocInvDiagEl(Q);

		if (LASResult() == LASOK && *Q->ElSorted && !(*Q->ZeroInDiag)
		    && !IsZero(Q->MultiplD)) {
                    /* assignment of auxiliary lokal variables */
                    QLen = Q->Len;
                    QEl = Q->El;
                    QInvDiagEl = Q->InvDiagEl;
                    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 */
                    MultiplD = 1.0 / Q->MultiplD;  /* attention here !!! */
                    MultiplU = Q->MultiplU;
                    MultiplL = Q->MultiplL;
                    MultiplV = V->Multipl;
                    MultiplDV = V->Multipl / Q->MultiplD;  /* attention here !!! */
                    MultiplDIsZero = IsZero(MultiplD);
                    MultiplUIsZero = IsZero(MultiplU);
                    MultiplLIsZero = IsZero(MultiplL);
                    MultiplDIsOne = IsOne(MultiplD);
                    MultiplUIsOne = IsOne(MultiplU);
                    MultiplLIsOne = IsOne(MultiplL);
                    MultiplVIsOne = IsOne(MultiplV);
                    MultiplDVIsOne = IsOne(MultiplDV);

                    /* multiplication of the vector V by the inverse matrix
                       of the diagonal of M */
                    if (!MultiplDIsZero && MultiplUIsZero && MultiplLIsZero) {
                        if (MultiplDVIsOne) {
			   for_AllCmp
			        VResCmp[Ind] = VCmp[Ind] * QInvDiagEl[Ind];
                        } else {
			   for_AllCmp
			        VResCmp[Ind] = MultiplDV * VCmp[Ind] * QInvDiagEl[Ind];
                        }
                    }

                    /* multiplication of the vector V by the inverse matrix
                       of the upper triangular part of M */
                    if (!MultiplUIsZero && MultiplLIsZero) {
                        if ((!Q->Symmetry && Q->ElOrder == Rowws)
                            || (Q->Symmetry && Q->ElOrder == Rowws)) {
                            for (Row = Dim; Row >= 1; Row--) {
                                Len = QLen[Row];
                                PtrEl = QEl[Row] + Len - 1;
                                Sum = 0.0;
        	  	        for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos > Row;
                                    ElCount--) {
                                    Sum -= (*PtrEl).Val * VResCmp[(*PtrEl).Pos];
                                    PtrEl--;
                                }
                                if (!MultiplUIsOne)
                                    Sum *= MultiplU;
                                if (MultiplVIsOne)
                                    Sum += VCmp[Row];
                                else
                                    Sum += MultiplV * VCmp[Row];
                                if (MultiplDIsOne) {
                                    VResCmp[Row] = Sum * QInvDiagEl[Row];
                                } else {
                                    VResCmp[Row] = Sum * MultiplD * QInvDiagEl[Row];
                                }
                            }
                        }
                        if ((!Q->Symmetry && Q->ElOrder == Clmws)
                            || (Q->Symmetry && Q->ElOrder == Clmws)) {
                            for (Clm = Dim; Clm >= 1; Clm--) {
                                Sum = VResCmp[Clm];
                                if (!MultiplUIsOne)
                                    Sum *= MultiplU;
                                if (MultiplVIsOne)
                                    Sum += VCmp[Clm];
                                else
                                    Sum += MultiplV * VCmp[Clm];
                                if (MultiplDIsOne) {
                                    Cmp = Sum * QInvDiagEl[Clm];
                                } else {
                                    Cmp = Sum * MultiplD * QInvDiagEl[Clm];
                                }
                                VResCmp[Clm] = Cmp;

                                Len = QLen[Clm];
                                PtrEl = QEl[Clm];
                                for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos < Clm;
                                    ElCount--) {
                                    VResCmp[(*PtrEl).Pos] -= (*PtrEl).Val * Cmp;
                                    PtrEl++;
                                }
                            }
                        }
                    }

                    /* multiplication of the vector V by the inverse matrix 
                       of the lower triangular part of M */
                    if (MultiplUIsZero && !MultiplLIsZero) {
                        if ((!Q->Symmetry && Q->ElOrder == Rowws)
                            || (Q->Symmetry && Q->ElOrder == Clmws)) {
                            for (Row = 1; Row <= Dim; Row++) {
                                Len = QLen[Row];
                                PtrEl = QEl[Row];
                                Sum = 0.0;
       			        for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos < Row;
                                    ElCount--) {
                                    Sum -= (*PtrEl).Val * VResCmp[(*PtrEl).Pos];
                                    PtrEl++;
                                }
                                if (!MultiplLIsOne)
                                    Sum *= MultiplL;
                                if (MultiplVIsOne)
                                    Sum += VCmp[Row];
                                else
                                    Sum += MultiplV * VCmp[Row];
                                if (MultiplDIsOne) {
                                    VResCmp[Row] = Sum * QInvDiagEl[Row];
                                } else {
                                    VResCmp[Row] = Sum * MultiplD * QInvDiagEl[Row];
                                }
                            }
                        }
                        if ((!Q->Symmetry && Q->ElOrder == Clmws)
                            || (Q->Symmetry && Q->ElOrder == Rowws)) {
                            for (Clm = 1; Clm <= Dim; Clm++) {
                                Sum = VResCmp[Clm];
                                if (!MultiplLIsOne)
                                    Sum *= MultiplL;
                                if (MultiplVIsOne)
                                    Sum += VCmp[Clm];
                                else
                                    Sum += MultiplV * VCmp[Clm];
                                if (MultiplDIsOne) {
                                    Cmp = Sum * QInvDiagEl[Clm];
                                } else {
                                    Cmp = Sum * MultiplD * QInvDiagEl[Clm];
                                }
                                VResCmp[Clm] = Cmp;

                                Len = QLen[Clm];
                                PtrEl = QEl[Clm] + Len - 1;
      			        for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos > Clm;
                                    ElCount--) {
                                    VResCmp[(*PtrEl).Pos] -= (*PtrEl).Val * Cmp;
                                    PtrEl--;
                                }
                            }
                        }
                    }
                    if (!MultiplUIsZero && !MultiplLIsZero) {
                        LASError(LASMulInvErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
                    }
                } else {
                    if (LASResult() == LASOK && !(*Q->ElSorted))
                        LASError(LASElNotSortedErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
                    if (LASResult() == LASOK && (*Q->ZeroInDiag || IsZero(Q->MultiplD)))
                        LASError(LASZeroInDiagErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
                }
	    } else {
                LASError(LASMemAllocErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
		if (VRes != NULL)
		    free(VRes);
            }
	    
	    if (VResName != NULL)
	        free(VResName);
	} else {
	    LASError(LASDimErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
	    VRes = NULL;
	}
    } else {
        VRes = NULL;
    }

    Q_Unlock(Q);
    V_Unlock(V);

    return(VRes);
}

Matrix *Transp_M(Matrix *M)
/* MRes = M^T, returns transposed matrix M */
{
    Matrix *MRes;

    char *MResName;

    M_Lock(M);
    
    if (LASResult() == LASOK) {
        MRes = (Matrix *)malloc(sizeof(Matrix));
	MResName = (char *)malloc((strlen(M_GetName(M)) + 10) * sizeof(char));
        if (MRes != NULL && MResName != NULL) {
	    sprintf(MResName, "(%s)^T", M_GetName(M));
            M_Constr(MRes, MResName,  M->ClmDim, M->RowDim, M->ElOrder, Tempor, False);
            if (LASResult() == LASOK) {
                if (M->ElOrder == Rowws)
                    MRes->ElOrder = Clmws;
                if (M->ElOrder == Clmws)
                    MRes->ElOrder = Rowws;
                MRes->Multipl = M->Multipl;
                if (M->Instance == Tempor && M->OwnData) {
                    M->OwnData = False;
                    MRes->OwnData = True;
                }
                MRes->Len = M->Len;
                MRes->El = M->El;
                MRes->ElSorted = M->ElSorted;
            }
        } else {
            LASError(LASMemAllocErr, "Transp_M", M_GetName(M), NULL, NULL);
	    if (MRes != NULL)
	        free(MRes);
        }
	    
	if (MResName != NULL)
	    free(MResName);
    } else {
        MRes = NULL;
    }

    M_Unlock(M);

    return(MRes);
}

QMatrix *Transp_Q(QMatrix *Q)
/* QRes = Q^T, returns transposed matrix Q */
{
    QMatrix *QRes;