factor.c

OpenFVM-v1.1 open source cfd code

                         *      (Q->ILU ~ (D + U)^T D^(-1) (D + U))  
                         *  and incomplete LU factorization
                         *      (Q->ILU ~ (D + L) D^(-1) (D + U))
                         *  respectively 
                         */
                        for (RoC = 1; RoC <= Dim && LASResult() == LASOK; RoC++) {
                            Len = Q->ILU->Len[RoC];
                            /* set index mapping */
                            PtrEl = Q->ILU->El[RoC];
                            LDim = 0;
                            for (ElCount = 0; ElCount < Len && (*PtrEl).Pos <= RoC;
                                ElCount++) {
                                IndexMapp[(*PtrEl).Pos] = ElCount + 1;
                                PtrEl++;
                                LDim++;
                            }

                            /* initialization of L */
                            for (j = 1; j <= LDim; j++)
                                for (k = 1; k <= LDim; k++)
                                    L[j][k] = 0.0;			 
                            
                            /* fill matrix L with elements which have influence
                               on the factorization */
                            PtrEl = Q->ILU->El[RoC];
                            for (ElCount = 0; ElCount < Len && (*PtrEl).Pos <= RoC;
                                ElCount++) {
                                /* for row or column RoC_ */
                                RoC_ = (*PtrEl).Pos;
                                Len_ = Q->ILU->Len[RoC_];
                                PtrEl_ = Q->ILU->El[RoC_];
                                for (ElCount_ = 0; ElCount_ < Len_ && (*PtrEl_).Pos <= RoC;
                                    ElCount_++) {
                                    L[ElCount + 1][IndexMapp[(*PtrEl_).Pos]] = (*PtrEl_).Val;
                                    PtrEl_++;
                                }
                                PtrEl++;
                            }
                            
                            /* factorize L */
                            if (Q->ILU->Symmetry) {
                                for (j = 1; j < LDim; j++)
                                    for (k = j + 1; k < LDim; k++)
                                        L[LDim][k] -= L[LDim][j] * L[k][j] / L[j][j];
                                for (j = 1; j < LDim; j++)
                                    L[LDim][LDim] -= L[j][LDim] * L[j][LDim] / L[j][j];
                            } else {
                                for (j = 1; j < LDim; j++)
                                    for (k = j + 1; k < LDim; k++)
                                        L[LDim][k] -= L[LDim][j] * L[j][k] / L[j][j];
                                for (j = 1; j < LDim; j++)
                                    for (k = j + 1; k < LDim; k++)
                                        L[k][LDim] -= L[j][LDim] * L[k][j] / L[j][j];
                                for (j = 1; j < LDim; j++)
                                    L[LDim][LDim] -= L[j][LDim] * L[LDim][j] / L[j][j];
                            }
                            if (IsZero(L[LDim][LDim]))
                                LASError(LASZeroPivotErr, "ILUFactor", Q_GetName(Q), NULL, NULL);
                            
                            /* set back factorized elements */
                            PtrEl = Q->ILU->El[RoC];
                            for (ElCount = 0; ElCount < Len && (*PtrEl).Pos <= RoC;
                                ElCount++) {
                                (*PtrEl).Val = L[LDim][IndexMapp[(*PtrEl).Pos]];
                                PtrEl++;
                            }
                            if (!Q->ILU->Symmetry) {
                                PtrEl = Q->ILU->El[RoC];
                                for (ElCount = 0; ElCount < Len && (*PtrEl).Pos < RoC;
                                    ElCount++) {
                                    /* for row or column RoC_ */
                                    RoC_ = (*PtrEl).Pos;
                                    Len_ = Q->ILU->Len[RoC_];
                                    PtrEl_ = Q->ILU->El[RoC_];
                                    ElFound = False;
                                    for (ElCount_ = 0; ElCount_ < Len_ && (*PtrEl_).Pos <= RoC;
                                        ElCount_++) {
                                        if ((*PtrEl_).Pos == RoC) {
                                            (*PtrEl_).Val = L[ElCount + 1][LDim];
                                            ElFound = True;
                                        }
                                        PtrEl_++;
                                    }
                                    if (!ElFound)
                                        LASError(LASILUStructErr, "ILUFactor", Q_GetName(Q), NULL, NULL);
                                    PtrEl++;
                                }
                            }
                            
                            /* reset index mapping */
                            PtrEl = Q->ILU->El[RoC];
                            for (ElCount = 0; ElCount < Len && (*PtrEl).Pos <= RoC;
                                ElCount++) {
                                IndexMapp[(*PtrEl).Pos] = 0;
                                PtrEl++;
                            }
                        }
                    }
                    
                    /* invert diagonal elements */
                    *Q->ILU->DiagElAlloc = False;
                    Q_AllocInvDiagEl(Q->ILU);
                } else {
                    LASError(LASMemAllocErr, "ILUFactor", Q_GetName(Q), NULL, NULL);
                }

                if (IndexMapp != NULL)
                    free(IndexMapp);
                if (L != NULL) {
                    for (j = 0; j <= MaxLen; j++) {
	                if (L[j] != NULL)
                            free(L[j]);
	            }
                    free(L);
	        }
            } else {
                if (LASResult() == LASOK && !(*Q->ILU->ElSorted))
                    LASError(LASElNotSortedErr, "ILUFactor", Q_GetName(Q), NULL, NULL);
                if (LASResult() == LASOK && *Q->ILU->ZeroInDiag)
                    LASError(LASZeroInDiagErr, "ILUFactor", Q_GetName(Q), NULL, NULL);
            }
        }
        
        if (LASResult() == LASOK) {
	    QILUName = (char *)malloc((strlen(Q_GetName(Q)) + 10) * sizeof(char));
	    if (QILUName != NULL) {
	        sprintf(QILUName, "ILU(%s)", Q_GetName(Q));
		Q_SetName(Q->ILU, QILUName);
	    
                /* element ordering of matrix Q which can be modified by Transp_Q
                   is valid for Q->ILU */
                Q->ILU->ElOrder = Q->ElOrder;
		
                /* check for multipliers of the matrix Q */
                if (Q->MultiplU == Q->MultiplD && Q->MultiplL == Q->MultiplD) {
                    /* multipliers of matrix Q are valid for Q->ILU too */
                    Q->ILU->MultiplD = Q->MultiplD;
                    Q->ILU->MultiplU = Q->MultiplU;
                    Q->ILU->MultiplL = Q->MultiplL;
                    QRes = Q->ILU;
                } else {
                    LASError(LASILUStructErr, "ILUFactor", Q_GetName(Q), NULL, NULL);
                    QRes = NULL;
                }
		
                free(QILUName);
	    } else {
                LASError(LASMemAllocErr, "ILUFactor", Q_GetName(Q), NULL, NULL);
                QRes = NULL;
	    }
        } else {
            QRes = NULL;
        }
    } else {
        QRes = NULL;
    }

    Q_Unlock(Q);

    return(QRes);
}