operats.c

OpenFVM-v1.1 open source cfd code

    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, "(%s)^T", 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;
                }
                if (!Q->Symmetry) {
                    if (Q->ElOrder == Rowws)
                        QRes->ElOrder = Clmws;
                    if (Q->ElOrder == Clmws)
                        QRes->ElOrder = Rowws;
                }
                QRes->MultiplD = Q->MultiplD;
                QRes->MultiplU = Q->MultiplL;
                QRes->MultiplL = Q->MultiplU;
                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;
                if (!Q->Symmetry) {
                    QRes->UnitRightKer = Q->UnitLeftKer;
                    QRes->RightKerCmp = Q->LeftKerCmp;
                    QRes->UnitLeftKer = Q->UnitRightKer;
                    QRes->LeftKerCmp = Q->RightKerCmp;
		} else {
                    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, "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, False);
            if (LASResult() == LASOK) {
                if (Q->Instance == Tempor && Q->OwnData) {
                    Q->OwnData = False;
                    QRes->OwnData = True;
                }
                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, False);
            if (LASResult() == LASOK) {
                if (Q->Instance == Tempor && Q->OwnData) {
                    Q->OwnData = False;
                    QRes->OwnData = True;
                }
                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, False);
            if (LASResult() == LASOK) {
                if (Q->Instance == Tempor && Q->OwnData) {
                    Q->OwnData = False;
                    QRes->OwnData = True;
                }
                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);
}

double l1Norm_V(Vector *V)
/* SRes = l1-Norm of the vector V */
{
    double SRes;

    double Sum;
    size_t Dim, Ind;
    Real *VCmp;

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

    V_Unlock(V);

    return(SRes);
}

double l2Norm_V(Vector *V)
/* SRes = l2-Norm of the vector V */
{
    double SRes;

    double Sum, Cmp;
    size_t Dim, Ind;
    Real *VCmp;

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

    V_Unlock(V);

    return(SRes);
}

double MaxNorm_V(Vector *V)
/* SRes = max-Norm of the vector V */
{
    double SRes;

    double MaxCmp, Cmp;
    size_t Dim, Ind;
    Real *VCmp;

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

    V_Unlock(V);

    return(SRes);
}

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

    double Sum, Mean;
    size_t Dim, Ind;
    Real *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 / (double)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);
}

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

    double Sum, Mean;
    size_t Dim, Ind;
    Real *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 / (double)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);
}