qmatrix.c

OpenFVM-v1.1 open source cfd code

/****************************************************************************/
/*                                qmatrix.c                                 */
/****************************************************************************/
/*                                                                          */
/* type QMATRIX                                                             */
/*                                                                          */
/* Copyright (C) 1992-1996 Tomas Skalicky. All rights reserved.             */
/*                                                                          */
/****************************************************************************/
/*                                                                          */
/*        ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS          */
/*              OF THE COPYRIGHT NOTICE (SEE FILE COPYRGHT.H)               */
/*                                                                          */
/****************************************************************************/

#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "qmatrix.h"
#include "errhandl.h"
#include "copyrght.h"

static ElType ZeroEl = { 0, 0.0 };

static int ElCompar(const void *El1, const void *El2);

void Q_Constr(QMatrix *Q, char *Name, size_t Dim, Boolean Symmetry,
              ElOrderType ElOrder, InstanceType Instance, Boolean OwnData)
/* constructor of the type QMatrix */
{
    size_t RoC;

    Q->Name = (char *)malloc((strlen(Name) + 1) * sizeof(char));
    if (Q->Name != NULL)
        strcpy(Q->Name, Name);
    else
        LASError(LASMemAllocErr, "Q_Constr", Name, NULL, NULL);
    Q->Dim = Dim;
    Q->Symmetry = Symmetry;
    Q->ElOrder = ElOrder;
    Q->Instance = Instance;
    Q->LockLevel = 0;
    Q->MultiplD = 1.0;
    Q->MultiplU = 1.0;
    Q->MultiplL = 1.0;
    Q->OwnData = OwnData;
    if (OwnData) {
        if (LASResult() == LASOK) {
	    Q->Len = (size_t *)malloc((Dim + 1) * sizeof(size_t));
	    Q->El = (ElType **)malloc((Dim + 1) * sizeof(ElType *));
	    Q->ElSorted = (Boolean *)malloc(sizeof(Boolean));
	    Q->DiagElAlloc = (Boolean *)malloc(sizeof(Boolean));
	    Q->DiagEl = (ElType **)malloc((Dim + 1) * sizeof(ElType *));
	    Q->ZeroInDiag = (Boolean *)malloc(sizeof(Boolean));
            Q->InvDiagEl = (Real *)malloc((Dim + 1) * sizeof(Real));
	    Q->ILUExists = (Boolean *)malloc(sizeof(Boolean));
            Q->ILU = (QMatrix *)malloc(sizeof(QMatrix));
	    if (Q->Len != NULL && Q->El != NULL && Q->ElSorted != NULL
	        && Q->DiagElAlloc != NULL && Q->DiagEl != NULL && Q->ZeroInDiag != NULL
	        && Q->InvDiagEl != NULL && Q->ILUExists != NULL && Q->ILU != NULL) {
                for (RoC = 1; RoC <= Dim; RoC++) {
                    Q->Len[RoC] = 0;
                    Q->El[RoC] = NULL;
                    Q->DiagEl[RoC] = NULL;
                    Q->InvDiagEl[RoC] = 0.0;
                }
                *Q->ElSorted = False;
                *Q->DiagElAlloc = False;
                *Q->ZeroInDiag = True;
                *Q->ILUExists = False;
            } else {
	        LASError(LASMemAllocErr, "Q_Constr", Name, NULL, NULL);
            }
        } else {
	    Q->Len = NULL;
	    Q->El = NULL;
	    Q->ElSorted = NULL;
	    Q->DiagElAlloc = NULL;
	    Q->DiagEl = NULL;
	    Q->ZeroInDiag = NULL;
            Q->InvDiagEl = NULL;
	    Q->ILUExists = NULL;
	    Q->ILU = NULL;
        }
    }
    Q->UnitRightKer = False;
    Q->RightKerCmp = NULL;
    Q->UnitLeftKer = False;
    Q->LeftKerCmp = NULL;
    Q->EigenvalInfo = NULL;
}

void Q_Destr(QMatrix *Q)
/* destructor of the type QMatrix */
{
    size_t Dim, RoC;

    if (Q->Name != NULL)
        free(Q->Name);
    Dim = Q->Dim;
    if (Q->OwnData) {
	if (Q->Len != NULL && Q->El != NULL) {
            for (RoC = 1; RoC <= Dim; RoC++) {
                if (Q->Len[RoC] > 0) {
                    if (Q->El[RoC] != NULL)
                        free(Q->El[RoC]);
                }
            }
        }
        if (Q->Len != NULL) {
            free(Q->Len);
            Q->Len = NULL;
        }
        if (Q->El != NULL) {
            free(Q->El);
            Q->El = NULL;
        }
        if (Q->ElSorted != NULL) {
            free(Q->ElSorted);
            Q->ElSorted = NULL;
        }
        if (Q->DiagElAlloc != NULL) {
            free(Q->DiagElAlloc);
            Q->DiagElAlloc = NULL;
        }
        if (Q->DiagEl != NULL) {
            free(Q->DiagEl);
            Q->DiagEl = NULL;
        }
        if (Q->ZeroInDiag != NULL) {
            free(Q->ZeroInDiag);
            Q->ZeroInDiag = NULL;
        }
        if (Q->InvDiagEl != NULL) {
            free(Q->InvDiagEl);
            Q->InvDiagEl = NULL;
        }
        if (Q->ILUExists != NULL && Q->ILU != NULL) {
            if (*Q->ILUExists) 
                Q_Destr(Q->ILU);
        }
        if (Q->ILUExists != NULL) {
            free(Q->ILUExists);
            Q->ILUExists = NULL;
        }
        if (Q->ILU != NULL) {
            free(Q->ILU);
            Q->ILU = NULL;
        }
    }
    if (Q->RightKerCmp != NULL) {
        free(Q->RightKerCmp);
        Q->RightKerCmp = NULL;
    }
    if (Q->LeftKerCmp != NULL) {
        free(Q->LeftKerCmp);
        Q->LeftKerCmp = NULL;
    }
    if (Q->EigenvalInfo != NULL) {
        free(Q->EigenvalInfo);
        Q->EigenvalInfo = NULL;
    }
}

void Q_SetName(QMatrix *Q, char *Name)
/* (re)set name of the matrix Q */
{
    if (LASResult() == LASOK) {
        free(Q->Name);
        Q->Name = (char *)malloc((strlen(Name) + 1) * sizeof(char));
        if (Q->Name != NULL)
            strcpy(Q->Name, Name);
        else
            LASError(LASMemAllocErr, "Q_SetName", Name, NULL, NULL);
    }
}

char *Q_GetName(QMatrix *Q)
/* returns the name of the matrix Q */
{
    if (LASResult() == LASOK)
        return(Q->Name);
    else
        return("");
}

size_t Q_GetDim(QMatrix *Q)
/* returns the dimension of the matrix Q */
{
    size_t Dim;

    if (LASResult() == LASOK)
        Dim = Q->Dim;
    else
        Dim = 0;
    return(Dim);
}

Boolean Q_GetSymmetry(QMatrix *Q)
/* returns True if Q is symmetric otherwise False */
{
    Boolean Symmetry;

    if (LASResult() == LASOK) {
        Symmetry = Q->Symmetry;
    } else {
        Symmetry = (Boolean)0;
    }
    return(Symmetry);
}

ElOrderType Q_GetElOrder(QMatrix *Q)
/* returns element order of the matrix Q */
{
    ElOrderType ElOrder;

    if (LASResult() == LASOK) {
        ElOrder = Q->ElOrder;
    } else {
        ElOrder = (ElOrderType)0;
    }
    return(ElOrder);
}

void Q_SetLen(QMatrix *Q, size_t RoC, size_t Len)
/* set the lenght of a row or column of the matrix Q */
{
    size_t ElCount;
    ElType *PtrEl;

    if (LASResult() == LASOK) {
        if (Q->Instance == Normal && RoC > 0 && RoC <= Q->Dim) {
            Q->Len[RoC] = Len;

            PtrEl = Q->El[RoC];

            if (PtrEl != NULL) {
                free(PtrEl);
		PtrEl = NULL;
	    }

            if (Len > 0) {
                PtrEl = (ElType *)malloc(Len * sizeof(ElType));
                Q->El[RoC] = PtrEl;

                if (PtrEl != NULL) {
                    for (ElCount = Len; ElCount > 0; ElCount--) {
                        *PtrEl = ZeroEl;
                        PtrEl++;
                    }
                } else {
                    LASError(LASMemAllocErr, "Q_SetLen", Q->Name, NULL, NULL);
                }
            } else {
                Q->El[RoC] = NULL;
            }
        } else {
            if (Q->Instance != Normal)
                LASError(LASLValErr, "Q_SetLen", Q->Name, NULL, NULL);
            else
                LASError(LASRangeErr, "Q_SetLen", Q->Name, NULL, NULL);
        }
    }
}

size_t Q_GetLen(QMatrix *Q, size_t RoC)
/* returns the lenght of a row or column of the matrix Q */
{
    size_t Len;

    if (LASResult() == LASOK) {
        if (RoC > 0 && RoC <= Q->Dim) {
            Len = Q->Len[RoC];
        } else {
            LASError(LASRangeErr, "Q_GetLen", Q->Name, NULL, NULL);
            Len = 0;
        }
    } else {
        Len = 0;
    }
    return(Len);
}

void Q_SetEntry(QMatrix *Q, size_t RoC, size_t Entry, size_t Pos, Real Val)
/* set a new matrix entry */
{
    if (LASResult() == LASOK) {
        if ((RoC > 0 && RoC <= Q->Dim && Pos > 0 && Pos <= Q->Dim) &&
            (Entry < Q->Len[RoC])) {
            Q->El[RoC][Entry].Pos = Pos;
            Q->El[RoC][Entry].Val = Val;
        } else {
            LASError(LASRangeErr, "Q_SetEntry", Q->Name, NULL, NULL);
        }
    }
}

size_t Q_GetPos(QMatrix *Q, size_t RoC, size_t Entry)
/* returns the position of a matrix element */
{
    size_t Pos;

    if (LASResult() == LASOK)
        if (RoC > 0 && RoC <= Q->Dim && Entry < Q->Len[RoC]) {
            Pos = Q->El[RoC][Entry].Pos;
        } else {
            LASError(LASRangeErr, "Q_GetPos", Q->Name, NULL, NULL);
            Pos = 0;
        }
    else
        Pos = 0;
    return(Pos);
}

Real Q_GetVal(QMatrix *Q, size_t RoC, size_t Entry)
/* returns the value of a matrix element */
{
    Real Val;

    if (LASResult() == LASOK)
        if (RoC > 0 && RoC <= Q->Dim && Entry < Q->Len[RoC]) {
            Val = Q->El[RoC][Entry].Val;
        } else {
            LASError(LASRangeErr, "Q_GetVal", Q->Name, NULL, NULL);
            Val = 0.0;
        }
    else
        Val = 0.0;
    return(Val);
}

void Q_AddVal(QMatrix *Q, size_t RoC, size_t Entry, Real Val)
/* add a value to a matrix entry */
{
    if (LASResult() == LASOK) {
        if ((RoC > 0 && RoC <= Q->Dim) && (Entry < Q->Len[RoC]))
            Q->El[RoC][Entry].Val += Val;
        else
            LASError(LASRangeErr, "Q_AddVal", Q->Name, NULL, NULL);