operats.c

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

                    if (!Q->Symmetry && Q->ElOrder == Rowws) {
                        if (MultiplDULVEquals) {
                            for (Row = 1; Row <= Dim; Row++) {
                                Len = QLen[Row];
                                PtrEl = QEl[Row];
		    	        Sum = 0.0;
			        for (ElCount = Len; ElCount > 0; ElCount--) {
			            Sum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                    PtrEl++;
			        }
                                if (MultiplDVIsOne)
			            VResCmp[Row] = Sum;
                                else
			            VResCmp[Row] = MultiplDV * Sum;
  			    }
  			} else {
                            for (Row = 1; Row <= Dim; Row++) {
                                Len = QLen[Row];
                                Sum = 0.0;
                                if (!MultiplLVIsZero) {
                                    PtrEl = QEl[Row];
		    	            PartSum = 0.0;
  			            for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos < Row;
                                        ElCount--) {
			                PartSum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                        PtrEl++;
			            }
                                    if (MultiplLVIsOne)
			                Sum += PartSum;
                                    else
			                Sum += MultiplLV * PartSum;
                                }
                                if (!MultiplDVIsZero) {
                                    if (MultiplDVIsOne)
			                Sum += (*QDiagEl[Row]).Val * VCmp[Row];
                                    else
			                Sum += MultiplDV * (*QDiagEl[Row]).Val * VCmp[Row];
                                }
                                if (!MultiplUVIsZero) {
                                    PtrEl = QEl[Row] + Len - 1;
		    	            PartSum = 0.0;
  			            for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos > Row;
                                        ElCount--) {
			                PartSum += (*PtrEl).Val * VCmp[(*PtrEl).Pos];
                                        PtrEl--;
			            }
                                    if (MultiplUVIsOne)
			                Sum += PartSum;
                                    else
			                Sum += MultiplUV * PartSum;
                                }
			        VResCmp[Row] = Sum;
  			    }
                        }
                    }
                    if (!Q->Symmetry && Q->ElOrder == Clmws) {
                        if (MultiplDULVEquals) {
 			    for (Clm = 1; Clm <= Dim; Clm++) {
                                Len = QLen[Clm];
                                PtrEl = QEl[Clm];
                                if (MultiplDVIsOne)
   			            Cmp = VCmp[Clm];
                                else
        		            Cmp = MultiplDV * VCmp[Clm];
			        for (ElCount = Len; ElCount > 0; ElCount--) {
				    VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * Cmp;
                                    PtrEl++;
			        }
			    }
  			} else {
                            for (Clm = 1; Clm <= Dim; Clm++) {
                                Len = QLen[Clm];
		                Cmp = VCmp[Clm];
                                if (!MultiplUVIsZero) {
                                    PtrEl = QEl[Clm];
                                    if (MultiplUVIsOne)
                                        PartCmp = Cmp;
                                    else
                                        PartCmp = MultiplUV * Cmp;
  			            for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos < Clm;
                                        ElCount--) {
 				        VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * PartCmp;
                                        PtrEl++;
			            }
                                }
                                if (!MultiplDVIsZero) {
                                    if (MultiplDVIsOne)
 				        VResCmp[Clm] += (*QDiagEl[Clm]).Val * Cmp;
                                    else
 				        VResCmp[Clm] += MultiplDV * (*QDiagEl[Clm]).Val * Cmp;
                                }
                                if (!MultiplLVIsZero) {
                                    PtrEl = QEl[Clm] + Len - 1;
                                    if (MultiplLVIsOne)
                                        PartCmp = Cmp;
                                    else
                                        PartCmp = MultiplLV * Cmp;
  			            for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos > Clm;
                                        ElCount--) {
 				        VResCmp[(*PtrEl).Pos] += (*PtrEl).Val * Cmp;
                                        PtrEl--;
			            }
                                }
  			    }
                        }
                    }
                } else {
                    if (LASResult() == LASOK && !(*Q->ElSorted))
                        LASError(LASElNotSortedErr, "Mul_QV", Q_GetName(Q), V_GetName(V), NULL);
                }
	    } else {
		LASError(LASMemAllocErr, "Mul_QV", Q_GetName(Q), V_GetName(V), NULL);
		if (VRes != NULL)
		    free(VRes);
            }
	    
	    if (VResName != NULL)
	        free(VResName);
	} else {
	    LASError(LASDimErr, "Mul_QV", Q_GetName(Q), V_GetName(V), NULL);
	    VRes = NULL;
	}
    } else {
        VRes = NULL;
    }

    Q_Unlock(Q);
    V_Unlock(V);

    return(VRes);
}

QVector *MulInv_QV(QMatrix *Q, QVector *V)
/* VRes = Q^(-1) * V, this operation is limited to diagonal or triangular
   matrices */
{
    QVector *VRes;

    char *VResName;
    _LPDouble MultiplD, MultiplU, MultiplL, MultiplV, MultiplDV;
    _LPDouble Sum, Cmp;
    size_t Dim, Row, Clm, Ind, Len, ElCount;
    size_t *QLen;
    _LPBoolean MultiplDIsZero, MultiplUIsZero, MultiplLIsZero;
    _LPBoolean MultiplDIsOne, MultiplUIsOne, MultiplLIsOne, MultiplVIsOne,
        MultiplDVIsOne;
    ElType **QEl, *PtrEl;
    _LPNumber *QInvDiagEl;
    _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)) + 20)
		           * sizeof(char));
            if (VRes != NULL && VResName != NULL) {
	        sprintf(VResName, "(%s)^(-1) * (%s)", Q_GetName(Q), V_GetName(V));
                V_Constr(VRes, VResName, Dim, Tempor, _LPTrue);

                /* sort of elements and allocation of the inverse of diagonal elements
                   of the matrix Q */
                Q_SortEl(Q);
                Q_AllocInvDiagEl(Q);

		if (LASResult() == LASOK && *Q->ElSorted && !(*Q->ZeroInDiag)
		    && !_LPIsZeroNumber(Q->MultiplD)) {
                    /* assignment of auxiliary lokal variables */
                    QLen = Q->Len;
                    QEl = Q->El;
                    QInvDiagEl = Q->InvDiagEl;
                    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 */
                    MultiplD = 1.0 / Q->MultiplD;  /* attention here !!! */
                    MultiplU = Q->MultiplU;
                    MultiplL = Q->MultiplL;
                    MultiplV = V->Multipl;
                    MultiplDV = V->Multipl / Q->MultiplD;  /* attention here !!! */
		    /*
		      these casts are not necessary when _LPBoolean is a bool,
		      but when _LPBoolean is a struct it avoids a warning
		      on some compilers
		    */
                    MultiplDIsZero = (_LPBoolean) _LPIsZeroNumber(MultiplD);
                    MultiplUIsZero = (_LPBoolean) _LPIsZeroNumber(MultiplU);
                    MultiplLIsZero = (_LPBoolean) _LPIsZeroNumber(MultiplL);
                    MultiplDIsOne  = (_LPBoolean) _LPIsOneNumber(MultiplD);
                    MultiplUIsOne  = (_LPBoolean) _LPIsOneNumber(MultiplU);
                    MultiplLIsOne  = (_LPBoolean) _LPIsOneNumber(MultiplL);
                    MultiplVIsOne  = (_LPBoolean) _LPIsOneNumber(MultiplV);
                    MultiplDVIsOne = (_LPBoolean) _LPIsOneNumber(MultiplDV);

                    /* multiplication of the vector V by the inverse matrix
                       of the diagonal of M */
                    if (!MultiplDIsZero && MultiplUIsZero && MultiplLIsZero) {
                        if (MultiplDVIsOne) {
			   for_AllCmp
			        VResCmp[Ind] = VCmp[Ind] * QInvDiagEl[Ind];
                        } else {
			   for_AllCmp
			        VResCmp[Ind] = MultiplDV * VCmp[Ind] * QInvDiagEl[Ind];
                        }
                    }

                    /* multiplication of the vector V by the inverse matrix
                       of the upper triangular part of M */
                    if (!MultiplUIsZero && MultiplLIsZero) {
                        if ((!Q->Symmetry && Q->ElOrder == Rowws)
                            || (Q->Symmetry && Q->ElOrder == Rowws)) {
                            for (Row = Dim; Row >= 1; Row--) {
                                Len = QLen[Row];
                                PtrEl = QEl[Row] + Len - 1;
                                Sum = 0.0;
        	  	        for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos > Row;
                                    ElCount--) {
                                    Sum -= (*PtrEl).Val * VResCmp[(*PtrEl).Pos];
                                    PtrEl--;
                                }
                                if (!MultiplUIsOne)
                                    Sum *= MultiplU;
                                if (MultiplVIsOne)
                                    Sum += VCmp[Row];
                                else
                                    Sum += MultiplV * VCmp[Row];
                                if (MultiplDIsOne) {
                                    VResCmp[Row] = Sum * QInvDiagEl[Row];
                                } else {
                                    VResCmp[Row] = Sum * MultiplD * QInvDiagEl[Row];
                                }
                            }
                        }
                        if ((!Q->Symmetry && Q->ElOrder == Clmws)
                            || (Q->Symmetry && Q->ElOrder == Clmws)) {
                            for (Clm = Dim; Clm >= 1; Clm--) {
                                Sum = VResCmp[Clm];
                                if (!MultiplUIsOne)
                                    Sum *= MultiplU;
                                if (MultiplVIsOne)
                                    Sum += VCmp[Clm];
                                else
                                    Sum += MultiplV * VCmp[Clm];
                                if (MultiplDIsOne) {
                                    Cmp = Sum * QInvDiagEl[Clm];
                                } else {
                                    Cmp = Sum * MultiplD * QInvDiagEl[Clm];
                                }
                                VResCmp[Clm] = Cmp;

                                Len = QLen[Clm];
                                PtrEl = QEl[Clm];
                                for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos < Clm;
                                    ElCount--) {
                                    VResCmp[(*PtrEl).Pos] -= (*PtrEl).Val * Cmp;
                                    PtrEl++;
                                }
                            }
                        }
                    }

                    /* multiplication of the vector V by the inverse matrix 
                       of the lower triangular part of M */
                    if (MultiplUIsZero && !MultiplLIsZero) {
                        if ((!Q->Symmetry && Q->ElOrder == Rowws)
                            || (Q->Symmetry && Q->ElOrder == Clmws)) {
                            for (Row = 1; Row <= Dim; Row++) {
                                Len = QLen[Row];
                                PtrEl = QEl[Row];
                                Sum = 0.0;
       			        for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos < Row;
                                    ElCount--) {
                                    Sum -= (*PtrEl).Val * VResCmp[(*PtrEl).Pos];
                                    PtrEl++;
                                }
                                if (!MultiplLIsOne)
                                    Sum *= MultiplL;
                                if (MultiplVIsOne)
                                    Sum += VCmp[Row];
                                else
                                    Sum += MultiplV * VCmp[Row];
                                if (MultiplDIsOne) {
                                    VResCmp[Row] = Sum * QInvDiagEl[Row];
                                } else {
                                    VResCmp[Row] = Sum * MultiplD * QInvDiagEl[Row];
                                }
                            }
                        }
                        if ((!Q->Symmetry && Q->ElOrder == Clmws)
                            || (Q->Symmetry && Q->ElOrder == Rowws)) {
                            for (Clm = 1; Clm <= Dim; Clm++) {
                                Sum = VResCmp[Clm];
                                if (!MultiplLIsOne)
                                    Sum *= MultiplL;
                                if (MultiplVIsOne)
                                    Sum += VCmp[Clm];
                                else
                                    Sum += MultiplV * VCmp[Clm];
                                if (MultiplDIsOne) {
                                    Cmp = Sum * QInvDiagEl[Clm];
                                } else {
                                    Cmp = Sum * MultiplD * QInvDiagEl[Clm];
                                }
                                VResCmp[Clm] = Cmp;

                                Len = QLen[Clm];
                                PtrEl = QEl[Clm] + Len - 1;
      			        for (ElCount = Len; ElCount > 0 && (*PtrEl).Pos > Clm;
                                    ElCount--) {
                                    VResCmp[(*PtrEl).Pos] -= (*PtrEl).Val * Cmp;
                                    PtrEl--;
                                }
                            }
                        }
                    }
                    if (!MultiplUIsZero && !MultiplLIsZero) {
                        LASError(LASMulInvErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
                    }
                } else {
                    if (LASResult() == LASOK && !(*Q->ElSorted))
                        LASError(LASElNotSortedErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
                    if (LASResult() == LASOK && (*Q->ZeroInDiag || _LPIsZeroNumber(Q->MultiplD)))
                        LASError(LASZeroInDiagErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
                }
	    } else {
                LASError(LASMemAllocErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
		if (VRes != NULL)
		    free(VRes);
            }
	    
	    if (VResName != NULL)
	        free(VResName);
	} else {
	    LASError(LASDimErr, "MulInv_QV", Q_GetName(Q), V_GetName(V), NULL);
	    VRes = NULL;
	}
    } else {
        VRes = NULL;
    }

    Q_Unlock(Q);
    V_Unlock(V);

    return(VRes);
}

Matrix *Transp_M(Matrix *M)
/* MRes = M^T, returns transposed matrix M */
{
#if defined(_LP_USE_COMPLEX_NUMBERS)
    printf("ERROR: Transpose of a complex matrix not implemented.\n");
    abort();

    Matrix *MRes = NULL;
    return(MRes);
#else
    Matrix *MRes;

    char *MResName;

    M_Lock(M);
    
    if (LASResult() == LASOK) {
        MRes = (Matrix *)malloc(sizeof(Matrix));
	MResName = (char *)malloc((strlen(M_GetName(M)) + 10) * sizeof(char));
        if (MRes != NULL && MResName != NULL) {
	    sprintf(MResName, "(%s)^T", M_GetName(M));
            M_Constr(MRes, MResName,  M->ClmDim, M->RowDim, M->ElOrder, Tempor, _LPFalse);
            if (LASResult() == LASOK) {
                if (M->ElOrder == Rowws)
                    MRes->ElOrder = Clmws;
                if (M->ElOrder == Clmws)
                    MRes->ElOrder = Rowws;
                MRes->Multipl = M->Multipl;
                if (M->Instance == Tempor && M->OwnData) {
                    M->OwnData = _LPFalse;
                    MRes->OwnData = _LPTrue;
                }
                MRes->Len = M->Len;
                MRes->El = M->El;
                MRes->ElSorted = M->ElSorted;
            }
        } else {
            LASError(LASMemAllocErr, "Transp_M", M_GetName(M), NULL, NULL);
	    if (MRes != NULL)
	        free(MRes);
        }
	    
	if (MResName != NULL)
	    free(MResName);
    } else {
        MRes = NULL;
    }

    M_Unlock(M);

    return(MRes);
#endif
}

QMatrix *Transp_Q(QMatrix *Q)
/* QRes = Q^T, returns transposed 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, "(%s)^T", 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;
                }
                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;
            }