qmatrix.c

OpenFVM-v1.1 open source cfd code

    }
}

Real Q_GetEl(QMatrix *Q, size_t Row, size_t Clm)
/* returns the value of a matrix element (all matrix elements are considered) */
{
    Real Val;
    
    size_t Len, ElCount;
    ElType *PtrEl;

    if (LASResult() == LASOK) {
        if (Row > 0 && Row <= Q->Dim && Clm > 0 && Clm <= Q->Dim) {
            Val = 0.0;
            if (Q->Symmetry && Q->ElOrder == Rowws) {
                if (Clm >= Row) {
                    Len = Q->Len[Row];
                    PtrEl = Q->El[Row];
                    for (ElCount = Len; ElCount > 0; ElCount--) {
                        if ((*PtrEl).Pos == Clm)
                            Val = (*PtrEl).Val;
                        PtrEl++;
                    }
                } else {
                    Len = Q->Len[Clm];
                    PtrEl = Q->El[Clm];
                    for (ElCount = Len; ElCount > 0; ElCount--) {
                        if ((*PtrEl).Pos == Row)
                            Val = (*PtrEl).Val;
                        PtrEl++;
                    }
                }
            } else if (Q->Symmetry && Q->ElOrder == Clmws) {
                if (Clm >= Row) {
                    Len = Q->Len[Clm];
                    PtrEl = Q->El[Clm];
                    for (ElCount = Len; ElCount > 0; ElCount--) {
                        if ((*PtrEl).Pos == Row)
                            Val = (*PtrEl).Val;
                        PtrEl++;
                    }
                } else {
                    Len = Q->Len[Row];
                    PtrEl = Q->El[Row];
                    for (ElCount = Len; ElCount > 0; ElCount--) {
                        if ((*PtrEl).Pos == Clm)
                            Val = (*PtrEl).Val;
                        PtrEl++;
                    }
                }
            } else if (!Q->Symmetry && Q->ElOrder == Rowws) {
                Len = Q->Len[Row];
                PtrEl = Q->El[Row];
                for (ElCount = Len; ElCount > 0; ElCount--) {
                    if ((*PtrEl).Pos == Clm)
                        Val = (*PtrEl).Val;
                    PtrEl++;
                }
            } else if (!Q->Symmetry && Q->ElOrder == Clmws) {
                Len = Q->Len[Clm];
                PtrEl = Q->El[Clm];
                for (ElCount = Len; ElCount > 0; ElCount--) {
                    if ((*PtrEl).Pos == Row)
                        Val = (*PtrEl).Val;
                    PtrEl++;
                }
            }

            if (Row == Clmws)
                Val *= Q->MultiplD;
            if (Row > Clm)
                Val *= Q->MultiplU;
            if (Row > Clm)
                Val *= Q->MultiplL;
        } else {
            LASError(LASRangeErr, "Q_GetEl", Q->Name, NULL, NULL);
            Val = 0.0;
        }
    } else {
        Val = 0.0;
    }
    return(Val);
}

void Q_SortEl(QMatrix *Q)
/* sorts elements of a row or column in ascended order */
{
    size_t Dim, RoC;
    Boolean UpperOnly;

    if (LASResult() == LASOK && !(*Q->ElSorted)) {
        Dim = Q->Dim;
        UpperOnly = True;
        for (RoC = 1; RoC <= Dim; RoC++) {
            /* sort of elements by the quick sort algorithms */
            qsort((void *)Q->El[RoC], Q->Len[RoC], sizeof(ElType), ElCompar);

            /* test whether elements contained in upper triangular part
               (incl. diagonal) of the matrix only */
            if (Q->ElOrder == Rowws) {
                if (Q->El[RoC][0].Pos < RoC)
                    UpperOnly = False;
            }
            if (Q->ElOrder == Clmws) {
                if (Q->El[RoC][Q->Len[RoC] - 1].Pos > RoC)
                    UpperOnly = False;
            }
        }
        
        *Q->ElSorted = True;
        *Q->DiagElAlloc = False;
        *Q->ZeroInDiag = True;
        
        if (Q->Symmetry) {
            if(!UpperOnly)
                LASError(LASSymStorErr, "Q_SortEl", Q->Name, NULL, NULL);
        }
    }
}

void Q_AllocInvDiagEl(QMatrix *Q)
/* allocate pointers and compute inverse for diagonal elements of the matrix Q */
{
    size_t Dim, RoC, Len, ElCount;
    Boolean Found;
    ElType *PtrEl;

    if (LASResult() == LASOK && !(*Q->DiagElAlloc)) {
        Dim = Q->Dim;
        *Q->ZeroInDiag = False;
        if (Q->Symmetry && Q->ElOrder == Rowws) {
            for (RoC = 1; RoC <= Dim; RoC++) {
                if (Q->El[RoC][0].Pos == RoC) {
                    Q->DiagEl[RoC] = Q->El[RoC];
                } else {
                    *Q->ZeroInDiag = True;
                    Q->DiagEl[RoC] = &ZeroEl;
                }
            }
        }
        if (Q->Symmetry && Q->ElOrder == Clmws) {
            for (RoC = 1; RoC <= Dim; RoC++) {
                Len = Q->Len[RoC];
                if (Q->El[RoC][Len - 1].Pos == RoC) {
                    Q->DiagEl[RoC] = Q->El[RoC] + Len - 1;
                } else {
                    *Q->ZeroInDiag = True;
                    Q->DiagEl[RoC] = &ZeroEl;
                }
            }
        }
        if (!Q->Symmetry) {
            for (RoC = 1; RoC <= Dim; RoC++) {
                Found = False;
                Len = Q->Len[RoC];
                PtrEl = Q->El[RoC] + Len - 1;
                for (ElCount = Len; ElCount > 0; ElCount--) {
                    if ((*PtrEl).Pos == RoC) {
                        Found = True;
                        Q->DiagEl[RoC] = PtrEl;
                    }
                    PtrEl--;
                }
                if (!Found) {
                    *Q->ZeroInDiag = True;
                    Q->DiagEl[RoC] = &ZeroEl;
                }
            }
        }
        *Q->DiagElAlloc = True;
        
        if (!(*Q->ZeroInDiag)) {
            for (RoC = 1; RoC <= Dim; RoC++)
                Q->InvDiagEl[RoC] = 1.0 / (*Q->DiagEl[RoC]).Val;
        }
    }
}

static int ElCompar(const void *El1, const void *El2)
/* compares positions of two matrix elements */
{
    int Compar;

    Compar = 0;
    if (((ElType *)El1)->Pos < ((ElType *)El2)->Pos)
        Compar = -1;
    if (((ElType *)El1)->Pos > ((ElType *)El2)->Pos)
        Compar = +1;

    return(Compar);
}

void Q_SetKer(QMatrix *Q, Vector *RightKer, Vector *LeftKer)
/* defines the null space in the case of a singular matrix */
{
    double Sum, Mean, Cmp, Norm;
    size_t Dim, Ind;
    Real *KerCmp;

    V_Lock(RightKer);
    V_Lock(LeftKer);
    
    if (LASResult() == LASOK) {
	if (Q->Dim == RightKer->Dim && (Q->Symmetry || Q->Dim == LeftKer->Dim)) {
  	    Dim = Q->Dim;
	    
	    /* release array for old null space components when it exists */ 
	    if (Q->RightKerCmp != NULL) {
	        free(Q->RightKerCmp);
	        Q->RightKerCmp = NULL;
	    }
	    if (Q->LeftKerCmp != NULL) {
	        free(Q->LeftKerCmp);
	        Q->LeftKerCmp = NULL;
	    }
	    Q->UnitRightKer = False;
	    Q->UnitLeftKer = False;
	    
	    /* right null space */
            KerCmp = RightKer->Cmp;
	    /* test whether the matrix Q has a unit right null space */
	    Sum = 0.0;
	    for(Ind = 1; Ind <= Dim; Ind++)
	        Sum += KerCmp[Ind];
	    Mean = Sum / (double)Dim;
	    Q->UnitRightKer = True;
	    if (!IsZero(Mean)) {
	        for(Ind = 1; Ind <= Dim; Ind++)
                    if (!IsOne(KerCmp[Ind] / Mean))
		        Q->UnitRightKer = False;
            } else {
	        Q->UnitRightKer = False;
	    }
	    if (!Q->UnitRightKer) {
	        Sum = 0.0;
	        for(Ind = 1; Ind <= Dim; Ind++) {
	            Cmp = KerCmp[Ind];
	            Sum += Cmp * Cmp;
	        }
		Norm = sqrt(Sum);
		if (!IsZero(Norm)) {
                    Q->RightKerCmp = (Real *)malloc((Dim + 1) * sizeof(Real));
		    if (Q->RightKerCmp != NULL) {
		        for(Ind = 1; Ind <= Dim; Ind++)
	                    Q->RightKerCmp[Ind] = KerCmp[Ind] / Norm;
		    } else {
	                LASError(LASMemAllocErr, "Q_SetKer", Q->Name, RightKer->Name,
		            LeftKer->Name);
		    }
		}
	    }
	    
	    if (!Q->Symmetry) {
	        /* left null space */
                KerCmp = LeftKer->Cmp;
	        /* test whether the matrix Q has a unit left null space */
	        Sum = 0.0;
	        for(Ind = 1; Ind <= Dim; Ind++)
	            Sum += KerCmp[Ind];
	        Mean = Sum / (double)Dim;
	        Q->UnitLeftKer = True;
	        if (!IsZero(Mean)) {
	            for(Ind = 1; Ind <= Dim; Ind++)
                        if (!IsOne(KerCmp[Ind] / Mean))
		            Q->UnitLeftKer = False;
                } else {
	            Q->UnitLeftKer = False;
	        }
	        if (!Q->UnitLeftKer) {
	            Sum = 0.0;
	            for(Ind = 1; Ind <= Dim; Ind++) {
	                Cmp = KerCmp[Ind];
	                Sum += Cmp * Cmp;
	            }
		    Norm = sqrt(Sum);
		    if (!IsZero(Norm)) {
                        Q->LeftKerCmp = (Real *)malloc((Dim + 1) * sizeof(Real));
		        if (Q->LeftKerCmp != NULL) {
		            for(Ind = 1; Ind <= Dim; Ind++)
	                        Q->LeftKerCmp[Ind] = KerCmp[Ind] / Norm;
		        } else {
	                    LASError(LASMemAllocErr, "Q_SetKer", Q->Name,
				RightKer->Name, LeftKer->Name);
		        }
	    	    }
	        }
	    } else {
	    }
	} else {
	    LASError(LASDimErr, "Q_SetKer", Q->Name, RightKer->Name, LeftKer->Name);
	}
    }

    V_Unlock(RightKer);
    V_Unlock(LeftKer);
}

Boolean Q_KerDefined(QMatrix *Q)
/* returns True if Q is singular and the null space has been defined
   otherwise False */ 
{
    Boolean KerDefined;

    if (LASResult() == LASOK) {
        if ((Q->UnitRightKer || Q->RightKerCmp != NULL) && !IsZero(Q->MultiplD)
	    && IsOne(Q->MultiplU / Q->MultiplD) && IsOne(Q->MultiplL / Q->MultiplD))
	    KerDefined = True;
	else
	    KerDefined = False;
    } else {
        KerDefined = (Boolean)0;
    }
    return(KerDefined);
}

void **Q_EigenvalInfo(QMatrix *Q)
/* return address of the infos for eigenvalues */
{
    return(&(Q->EigenvalInfo));
}

void Q_Lock(QMatrix *Q)
/* lock the matrix Q */
{
    if (Q != NULL) 
        Q->LockLevel++;
}

void Q_Unlock(QMatrix *Q)
/* unlock the matrix Q */
{
    if (Q != NULL) {
        Q->LockLevel--;
        if (Q->Instance == Tempor && Q->LockLevel <= 0) {
            Q_Destr(Q); 
	    free(Q);
	}
    }
}