operats.c

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

        } else {
            LASError(LASMemAllocErr, "Transp_Q", Q_GetName(Q), NULL, NULL);
	    if (QRes != NULL)
	        free(QRes);
        }
	    
        if (QResName != NULL)
            free(QResName);
    } else {
        QRes = NULL;
    }

    Q_Unlock(Q);

    return(QRes);
}

QMatrix *Diag_Q(QMatrix *Q)
/* QRes = Diag(Q), returns the diagonal of the matrix Q */
{
    QMatrix *QRes;

    char *QResName;

    Q_Lock(Q);
    
    if (LASResult() == LASOK) {
        QRes = (QMatrix *)malloc(sizeof(QMatrix));
	QResName = (char *)malloc((strlen(Q_GetName(Q)) + 10) * sizeof(char));
        if (QRes != NULL && QResName != NULL) {
	    sprintf(QResName, "Diag(%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 = Q->MultiplD;
                QRes->MultiplU = 0.0;
                QRes->MultiplL = 0.0;
                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, "Diag_Q", Q_GetName(Q), NULL, NULL);
	    if (QRes != NULL)
	        free(QRes);
        }
	    
        if (QResName != NULL)
            free(QResName);
    } else {
        QRes = NULL;
    }

    Q_Unlock(Q);
 
    return(QRes);
}

QMatrix *Upper_Q(QMatrix *Q)
/* QRes = Upper(Q), returns the upper triagonal part of the matrix Q */
{
    QMatrix *QRes;

    char *QResName;

    Q_Lock(Q);
    
    if (LASResult() == LASOK) {
        QRes = (QMatrix *)malloc(sizeof(QMatrix));
	QResName = (char *)malloc((strlen(Q_GetName(Q)) + 10) * sizeof(char));
        if (QRes != NULL && QResName != NULL) {
	    sprintf(QResName, "Upper(%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 = 0.0;
                QRes->MultiplU = Q->MultiplU;
                QRes->MultiplL = 0.0;
                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, "Upper_Q", Q_GetName(Q), NULL, NULL);
	    if (QRes != NULL)
	        free(QRes);
        }
	    
        if (QResName != NULL)
            free(QResName);
    } else {
        QRes = NULL;
    }

    Q_Unlock(Q);

    return(QRes);
}

QMatrix *Lower_Q(QMatrix *Q)
/* QRes = Lower(Q), returns the lower triagonal part of the matrix Q */
{
    QMatrix *QRes;

    char *QResName;

    Q_Lock(Q);
    
    if (LASResult() == LASOK) {
        QRes = (QMatrix *)malloc(sizeof(QMatrix));
	QResName = (char *)malloc((strlen(Q_GetName(Q)) + 10) * sizeof(char));
        if (QRes != NULL && QResName != NULL) {
	    sprintf(QResName, "Lower(%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 = 0.0;
                QRes->MultiplU = 0.0;
                QRes->MultiplL = 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, "Lower_Q", Q_GetName(Q), NULL, NULL);
	    if (QRes != NULL)
	        free(QRes);
        }
	    
        if (QResName != NULL)
            free(QResName);
    } else {
        QRes = NULL;
    }

    Q_Unlock(Q);

    return(QRes);
}

_LPReal l1Norm_V(QVector *V)
/* SRes = l1-Norm of the vector V */
{
#if defined(_LP_USE_COMPLEX_NUMBERS)
    _LPReal SRes, Sum;

    _LPNumber Multipl;
    size_t Dim, Ind;
    _LPNumber *VCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
	Dim = V->Dim;
        VCmp = V->Cmp;
	Sum = 0.0;
	Multipl = V->Multipl;
	for_AllCmp
	   Sum += _LPfabs(VCmp[Ind]*Multipl);
	SRes = Sum;
    } else {
	SRes = 1.0;
    }

    V_Unlock(V);

    return(SRes);
#else
    _LPReal SRes;

    _LPReal Sum;
    size_t Dim, Ind;
    _LPNumber *VCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
	Dim = V->Dim;
        VCmp = V->Cmp;
	Sum = 0.0;
	for_AllCmp
	   Sum += _LPfabs(VCmp[Ind]);
        Sum *= V->Multipl;
	SRes = Sum;
    } else {
	SRes = 1.0;
    }

    V_Unlock(V);

    return(SRes);
#endif
}

_LPReal l2Norm_V(QVector *V)
/* SRes = l2-Norm of the vector V */
{
#if defined(_LP_USE_COMPLEX_NUMBERS)
    _LPReal SRes, Sum;
    _LPNumber Cmp, Multipl;
    size_t Dim, Ind;
    _LPNumber *VCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
	Dim = V->Dim;
        VCmp = V->Cmp;
	Multipl = V->Multipl;
	Sum = 0.0;
	for_AllCmp {
	    Cmp = VCmp[Ind]*Multipl;
	    Sum += _LPNormxNorm(Cmp);
	}
	SRes = sqrt(Sum);
    } else {
	SRes = 1.0;
    }

    V_Unlock(V);

    return(SRes);
#else
    _LPReal SRes, Sum;
    _LPNumber Cmp;
    size_t Dim, Ind;
    _LPNumber *VCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
	Dim = V->Dim;
        VCmp = V->Cmp;
	Sum = 0.0;
	for_AllCmp {
	    Cmp = VCmp[Ind];
	    Sum += _LPNormxNorm(Cmp);
	}
        Sum *= V->Multipl * V->Multipl;
	SRes = sqrt(Sum);
    } else {
	SRes = 1.0;
    }

    V_Unlock(V);

    return(SRes);
#endif
}

_LPReal MaxNorm_V(QVector *V)
/* SRes = max-Norm of the vector V */
{
#if defined(_LP_USE_COMPLEX_NUMBERS)
    _LPReal SRes;

    _LPNumber MaxCmp, Cmp, Multipl;
    size_t Dim, Ind;
    _LPNumber *VCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
	Dim = V->Dim;
        VCmp = V->Cmp;
	Multipl = V->Multipl;
	MaxCmp = 0.0;
	for_AllCmp {
	   Cmp = _LPAbsRealPart(VCmp[Ind] * Multipl);
	   if (_LPIsGreater(Cmp, MaxCmp)) 
	       MaxCmp = Cmp;
	}
	SRes = _LPRealPart(MaxCmp);
    } else {
	SRes = 1.0;
    }

    V_Unlock(V);

    return(SRes);
#else
    _LPReal SRes;

    _LPNumber MaxCmp, Cmp;
    size_t Dim, Ind;
    _LPNumber *VCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
	Dim = V->Dim;
        VCmp = V->Cmp;
	MaxCmp = 0.0;
	for_AllCmp {
	   Cmp = _LPAbsRealPart(VCmp[Ind]);
	   if (_LPIsGreater(Cmp, MaxCmp)) 
	       MaxCmp = Cmp;
	}
        MaxCmp *= V->Multipl;
	SRes = _LPRealPart(MaxCmp);
    } else {
	SRes = 1.0;
    }

    V_Unlock(V);

    return(SRes);
#endif
}

QVector *OrthoRightKer_VQ(QVector *V, QMatrix *Q)
/* orthogonalize vector V to the "right" null space of matrix Q */
{
    QVector *VRes;

    _LPDouble Sum, Mean;
    size_t Dim, Ind;
    _LPNumber *VCmp, *KerCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
        if (Q->UnitRightKer || Q->RightKerCmp != NULL) {
            if (V->Instance == Normal && V->Dim == Q->Dim) {
                Dim = V->Dim;
                VCmp = V->Cmp;
                KerCmp = Q->RightKerCmp;
	    
	        if (Q->UnitRightKer) {
	            Sum = 0.0;
   	            for_AllCmp 
	                Sum += VCmp[Ind];
	    	    Mean = Sum / (_LPDouble)Dim;
   	            for_AllCmp 
	                VCmp[Ind] -= Mean;
                } else {
	            Sum = 0.0;
   	            for_AllCmp
	                Sum += VCmp[Ind] * KerCmp[Ind];
   	            for_AllCmp 
	                VCmp[Ind] -= Sum * KerCmp[Ind];
    	        }	   
	   	    
                VRes = V;
            } else {
                if (V->Instance != Normal) 
                    LASError(LASLValErr, "OrthoRightKer_VQ", V_GetName(V), Q_GetName(Q), NULL);
                if (V->Dim != Q->Dim) 
                    LASError(LASDimErr, "OrthoRightKer_VQ", V_GetName(V), Q_GetName(Q), NULL);
                VRes = NULL;
            }
	} else {
            VRes = V;
	}
    } else {
        VRes = NULL;
    }

    V_Unlock(V);

    return(VRes);
}

QVector *OrthoLeftKer_VQ(QVector *V, QMatrix *Q)
/* orthogonalize vector V to the "left" null space of matrix Q */
{
    QVector *VRes;

    _LPDouble Sum, Mean;
    size_t Dim, Ind;
    _LPNumber *VCmp, *KerCmp;

    V_Lock(V);
    
    if (LASResult() == LASOK) {
        if (Q->UnitRightKer || Q->RightKerCmp != NULL) {
            if (V->Instance == Normal && V->Dim == Q->Dim) {
                Dim = V->Dim;
                VCmp = V->Cmp;
		if (Q->Symmetry)
                    KerCmp = Q->RightKerCmp;
		else
                    KerCmp = Q->LeftKerCmp;
	    
	        if ((Q->Symmetry && Q->UnitRightKer) || Q->UnitLeftKer) {
	            Sum = 0.0;
   	            for_AllCmp 
	                Sum += VCmp[Ind];
	    	    Mean = Sum / (_LPDouble)Dim;
   	            for_AllCmp 
	                VCmp[Ind] -= Mean;
                } else {
	            Sum = 0.0;
   	            for_AllCmp
	                Sum += VCmp[Ind] * KerCmp[Ind];
   	            for_AllCmp 
	                VCmp[Ind] -= Sum * KerCmp[Ind];
    	        }	   
	   	    
                VRes = V;
            } else {
                if (V->Instance != Normal) 
                    LASError(LASLValErr, "OrthoLeftKer_VQ", V_GetName(V), Q_GetName(Q), NULL);
                if (V->Dim != Q->Dim) 
                    LASError(LASDimErr, "OrthoLeftKer_VQ", V_GetName(V), Q_GetName(Q), NULL);
                VRes = NULL;
            }
	} else {
            VRes = V;
	}
    } else {
        VRes = NULL;
    }

    V_Unlock(V);

    return(VRes);
}