diff --git a/src/maths/sparse/spBuild.c b/src/maths/sparse/spBuild.c index ff767d0f2..c0baca585 100644 --- a/src/maths/sparse/spBuild.c +++ b/src/maths/sparse/spBuild.c @@ -402,7 +402,7 @@ spcFindElementInCol(MatrixPtr Matrix, ElementPtr *LastAddr, /* Element does not exist and must be created. */ if (CreateIfMissing) - return spcCreateElement( Matrix, Row, Col, LastAddr, NO, 0); + return spcCreateElementOld( Matrix, Row, Col, LastAddr, NO); else return NULL; } @@ -911,7 +911,114 @@ ElementPtr pElement, pCreatedElement; } +ElementPtr +spcCreateElementOld(MatrixPtr Matrix, int Row, int Col, + ElementPtr *LastAddr, int Fillin) +{ + ElementPtr pElement, pLastElement; + ElementPtr pCreatedElement; + /* Begin `spcCreateElement'. */ + + if (Matrix->RowsLinked) + { + /* Row pointers cannot be ignored. */ + + if (Fillin) + { + pElement = spcGetFillin( Matrix ); + Matrix->Fillins++; + } + else + { + pElement = spcGetElement( Matrix ); + Matrix->Elements++; + + Matrix->NeedsOrdering = YES; + } + if (pElement == NULL) return NULL; + + /* If element is on diagonal, store pointer in Diag. */ + if (Row == Col) Matrix->Diag[Row] = pElement; + + /* Initialize Element. */ + pCreatedElement = pElement; + pElement->Row = Row; + pElement->Col = Col; + pElement->Real = 0.0; + pElement->Imag = 0.0; +#if INITIALIZE + pElement->pInitInfo = NULL; +#endif + + /* Splice element into column. */ + pElement->NextInCol = *LastAddr; + *LastAddr = pElement; + + /* Search row for proper element position. */ + pElement = Matrix->FirstInRow[Row]; + pLastElement = NULL; + while (pElement != NULL) + { + /* Search for element row position. */ + if (pElement->Col < Col) + { + /* Have not reached desired element. */ + pLastElement = pElement; + pElement = pElement->NextInRow; + } + else pElement = NULL; + } + + /* Splice element into row. */ + pElement = pCreatedElement; + if (pLastElement == NULL) + { + /* Element is first in row. */ + pElement->NextInRow = Matrix->FirstInRow[Row]; + Matrix->FirstInRow[Row] = pElement; + } + else + { + /* Element is not first in row. */ + pElement->NextInRow = pLastElement->NextInRow; + pLastElement->NextInRow = pElement; + } + + } + else + { + /* Matrix has not been factored yet. Thus get element rather + * than fill-in. Also, row pointers can be ignored. */ + + /* Allocate memory for Element. */ + pElement = spcGetElement( Matrix ); + Matrix->Elements++; + if (pElement == NULL) return NULL; + + /* If element is on diagonal, store pointer in Diag. */ + if (Row == Col) Matrix->Diag[Row] = pElement; + + /* Initialize Element. */ + pCreatedElement = pElement; + pElement->Row = Row; +#if DEBUG + pElement->Col = Col; +#endif + pElement->Real = 0.0; + pElement->Imag = 0.0; +#if INITIALIZE + pElement->pInitInfo = NULL; +#endif + + /* Splice element into column. */ + pElement->NextInCol = *LastAddr; + *LastAddr = pElement; + } + + Matrix->Elements++; + return pCreatedElement; +} diff --git a/src/maths/sparse/spDefs.h b/src/maths/sparse/spDefs.h index 9bde6bba4..7aa8ee2b3 100644 --- a/src/maths/sparse/spDefs.h +++ b/src/maths/sparse/spDefs.h @@ -944,6 +944,8 @@ spcEXTERN ElementPtr spcFindElementInCol( MatrixPtr, ElementPtr*, int, int, int spcEXTERN ElementPtr spcFindDiag( MatrixPtr, int ); spcEXTERN ElementPtr spcCreateElement( MatrixPtr, int, int, ElementPtr*, ElementPtr*, int ); +spcEXTERN ElementPtr spcCreateElementOld( MatrixPtr, int, int, + ElementPtr*, int ); spcEXTERN void spcCreateInternalVectors( MatrixPtr ); spcEXTERN void spcLinkRows( MatrixPtr ); spcEXTERN void spcColExchange( MatrixPtr, int, int ); diff --git a/src/maths/sparse/spSMP.c b/src/maths/sparse/spSMP.c index d16f76966..1f308a36c 100644 --- a/src/maths/sparse/spSMP.c +++ b/src/maths/sparse/spSMP.c @@ -196,7 +196,7 @@ SMPmatrix *Matrix, double PivTol, double PivRel, int *NumSwaps) { - *NumSwaps = 0; + *NumSwaps = 1; spSetComplex( Matrix ); return spOrderAndFactor( Matrix, NULL, PivRel, PivTol, YES );