operats.c

一个用来实现偏微分方程中网格的计算库

                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);
}

QVector *Mul_SV(_LPDouble S, QVector *V)
/* VRes = S * V */
{
    QVector *VRes;

    char *VResName;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
        VRes = (QVector *)malloc(sizeof(QVector));
	VResName = (char *)malloc((strlen(V_GetName(V)) + 20) * sizeof(char));
        if (VRes != NULL && VResName != NULL) {
	    sprintf(VResName, "%12.5e * (%s)", _LPPrintFormat(S), V_GetName(V));
            V_Constr(VRes, VResName, V->Dim, Tempor, _LPFalse);
            if (LASResult() == LASOK) {
                VRes->Multipl = S * V->Multipl;
                if (V->Instance == Tempor && V->OwnData) {
                    V->OwnData = _LPFalse;
                    VRes->OwnData = _LPTrue;
                }
                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(_LPDouble 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)", _LPPrintFormat(S), M_GetName(M));
            M_Constr(MRes, MResName,  M->RowDim, M->ClmDim, M->ElOrder, Tempor, _LPFalse);
            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(_LPDouble 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)", _LPPrintFormat(S), Q_GetName(Q));
            Q_Constr(QRes, QResName, Q->Dim, Q->Symmetry, Q->ElOrder, Tempor, _LPFalse);
            if (LASResult() == LASOK) {
                if (Q->Instance == Tempor && Q->OwnData) {
                    Q->OwnData = _LPFalse;
                    QRes->OwnData = _LPTrue;
                }
                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);
}

_LPDouble Mul_VV(QVector *V1, QVector *V2)
/* S = V1 * V2 */
{
    _LPDouble SRes;

    size_t Dim, Ind;
    _LPNumber *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);
}

_LPDouble InnerProd_VV(QVector *V1, QVector *V2)
/* S = conjugate V1 * V2 */
{

/* S = V1 * V2 for real numbers*/
#ifndef _LP_USE_COMPLEX_NUMBERS
    return (Mul_VV(V1, V2));
#else

    _LPDouble SRes;

    size_t Dim, Ind;
    _LPNumber *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 += conj(V1Cmp[Ind]) * V2Cmp[Ind];
	    SRes *= V1->Multipl * V2->Multipl;
        } else {
            LASError(LASDimErr, "InnerProd_VV", V_GetName(V1),  V_GetName(V2), NULL);
            SRes = 1.0;
        }
    } else {
        SRes = 1.0;
    }

    V_Unlock(V1);
    V_Unlock(V2);

    return(SRes);

#endif /* ifdef _LP_USE_COMPLEX_NUMBERS */
}

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

    char *VResName;
    _LPDouble MultiplMV;
    _LPDouble Sum, Cmp;
    size_t RowDim, ClmDim, Row, Clm, Len, ElCount;
    size_t *MLen;
    _LPBoolean MultiplMVIsOne;
    ElType **MEl, *PtrEl;
    _LPNumber *VCmp, *VResCmp;

    M_Lock(M);
    V_Lock(V);
    
    if (LASResult() == LASOK) {
	if (M->ClmDim == V->Dim) {
            RowDim = M->RowDim;
            ClmDim = M->ClmDim;
	    VRes = (QVector *)malloc(sizeof(QVector));
	    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, _LPTrue);
		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 (_LPIsOneNumber(MultiplMV)) {
                        MultiplMVIsOne = _LPTrue;
                    } else {
                        MultiplMVIsOne = _LPFalse;
                    }

                    /* 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);
}

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

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

    Q_Lock(Q);
    V_Lock(V);
    
    if (LASResult() == LASOK) {
	if (Q->Dim == V->Dim) {
	    Dim = V->Dim;
	    VRes = (QVector *)malloc(sizeof(QVector));
	    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, _LPTrue);
		
                /* 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;
		    /*
		      these casts are not necessary when _LPBoolean is a bool,
		      but when _LPBoolean is a struct it avoids a warning
		      on some compilers
		    */
                    MultiplDVIsZero = (_LPBoolean) _LPIsZeroNumber(MultiplDV);
                    MultiplUVIsZero = (_LPBoolean) _LPIsZeroNumber(MultiplUV);
                    MultiplLVIsZero = (_LPBoolean) _LPIsZeroNumber(MultiplLV);
                    MultiplDVIsOne  = (_LPBoolean) _LPIsOneNumber(MultiplDV);
                    MultiplUVIsOne  = (_LPBoolean) _LPIsOneNumber(MultiplUV);
                    MultiplLVIsOne  = (_LPBoolean) _LPIsOneNumber(MultiplLV);
                    if (!_LPIsZeroNumber(MultiplDV) && _LPIsOneNumber(MultiplUV / MultiplDV)
			&& _LPIsOneNumber(MultiplLV / MultiplDV)) {
                        MultiplDULVEquals = _LPTrue;
                    } else {
                        MultiplDULVEquals = _LPFalse;
                    }

                    /* 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;
                                }
                                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++;
	                            }
                                }
  			    }
                        }
                    }