/* * Sparse1.4 is distributed as open-source software under the Berkeley * license model. Redistribution and use in source and binary forms, with * or without modification, are permitted provided that the following * conditions are met: * * Redistributions of source code must retain the original copyright notice, * this list of conditions and the following disclaimer. Redistributions * in binary form must reproduce the above copyright notice, this list * of conditions and the following disclaimer in the documentation and/or * other materials provided with the distribution. Neither the name of * the copyright holder nor the names of the authors may be used to endorse * or promote products derived from this software without specific prior * written permission. * * This software is provided by the copyright holders and contributors * ``as is'' and any express or implied warranties, including, but not * limited to, the implied warranties of merchantability and fitness for * a particular purpose are disclaimed. In no event shall the copyright * owner or contributors be liable for any direct, indirect, incidental, * special, exemplary, or consequential damages (including, but not * limited to, procurement of substitute goods or services; loss of use, * data, or profits; or business interruption) however caused and on any * theory of liability, whether in contract, strict liability, or tort * (including negligence or otherwise) arising in any way out of the use * of this software, even if advised of the possibility of such damage. */ /* * MATRIX BUILD MODULE * * Author: Advising professor: * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli * UC Berkeley * * This file contains the routines associated with clearing, loading and * preprocessing the matrix. */ #include #include #include static void Translate( MatrixPtr, int*, int* ); static void EnlargeMatrix( MatrixPtr, int ); static void ExpandTranslationArrays( MatrixPtr, int ); /*! * Sets every element of the matrix to zero and clears the error flag. * * \param eMatrix * Pointer to matrix that is to be cleared. */ /* >>> Local variables: * pElement (ElementPtr) * A pointer to the element being cleared. */ void spClear( spMatrix eMatrix ) { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement; register int I; /* Begin `spClear'. */ ASSERT_IS_SPARSE( Matrix ); /* Clear matrix. */ if (Matrix->PreviousMatrixWasComplex OR Matrix->Complex) { for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { pElement->Real = 0.0; pElement->Imag = 0.0; pElement = pElement->NextInCol; } } } else { for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { pElement->Real = 0.0; pElement = pElement->NextInCol; } } } /* Empty the trash. */ Matrix->TrashCan.Real = 0.0; Matrix->TrashCan.Imag = 0.0; Matrix->Error = spOKAY; Matrix->Factored = NO; Matrix->SingularCol = 0; Matrix->SingularRow = 0; Matrix->PreviousMatrixWasComplex = Matrix->Complex; return; } /*! * This routine is used to find an element given its indices. It will not * create it if it does not exist. * * \return * A pointer to the desired element, or \a NULL if it does not exist. * * \param eMatrix * Pointer to matrix. * \param Row * Row index for element. * \param Col * Column index for element. * * \see spGetElement() */ /* >>> Local variables: * pElement (ElementPtr) * Pointer to an element in the matrix. */ spElement * spFindElement( spMatrix eMatrix, int Row, int Col ) { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement; int StartAt = 0, Min = AG_INT_MAX; #define BorderRight 0 /* Start at left border, move right. */ #define BorderDown 1 /* Start at top border, move down. */ #define DiagRight 2 /* Start at diagonal, move right. */ #define DiagDown 3 /* Start at diagonal, move down. */ ASSERT_IS_SPARSE( Matrix ); vASSERT( Row >= 0 AND Col >= 0, "Negative row or column number" ); #if TRANSLATE Translate( Matrix, &Row, &Col ); if (Matrix->Error == spNO_MEMORY) return NULL; #endif #if NOT TRANSLATE #if NOT EXPANDABLE vASSERT( (Row <= Matrix->Size) AND (Col <= Matrix->Size), "Row or column number too large" ); #endif #endif /* Begin `spFindElement'. */ if (Row == Col) return &Matrix->Diag[Row]->Real; /* Determine where to start the search. */ if (Matrix->RowsLinked) { if ((Col >= Row) AND Matrix->Diag[Row]) { Min = Col - Row; StartAt = DiagRight; } else { Min = Col; StartAt = BorderRight; } } if ((Row >= Col) AND Matrix->Diag[Col]) { if (Row - Col < Min) StartAt = DiagDown; } else if (Row < Min) StartAt = BorderDown; /* Search column for element. */ if ((StartAt == BorderDown) OR (StartAt == DiagDown)) { if (StartAt == BorderDown) pElement = Matrix->FirstInCol[Col]; else pElement = Matrix->Diag[Col]; while ((pElement != NULL) AND (pElement->Row < Row)) pElement = pElement->NextInCol; if (pElement AND (pElement->Row == Row)) return &pElement->Real; else return NULL; } /* Search row for element. */ if (StartAt == BorderRight) pElement = Matrix->FirstInRow[Row]; else pElement = Matrix->Diag[Row]; while ((pElement != NULL) AND (pElement->Col < Col)) pElement = pElement->NextInRow; if (pElement AND (pElement->Col == Col)) return &pElement->Real; else return NULL; } /*! * Finds element [Row,Col] and returns a pointer to it. If element is * not found then it is created and spliced into matrix. This routine * is only to be used after spCreate() and before spMNA_Preorder(), * spFactor() or spOrderAndFactor(). Returns a pointer to the * real portion of an \a spElement. This pointer is later used by * \a spADD_xxx_ELEMENT to directly access element. * * \return * Returns a pointer to the element. This pointer is then used to directly * access the element during successive builds. * * \param eMatrix * Pointer to the matrix that the element is to be added to. * \param Row * Row index for element. Must be in the range of [0..Size] unless * the options \a EXPANDABLE or \a TRANSLATE are used. Elements placed in * row zero are discarded. In no case may \a Row be less than zero. * \param Col * Column index for element. Must be in the range of [0..Size] unless * the options \a EXPANDABLE or \a TRANSLATE are used. Elements placed in * column zero are discarded. In no case may \a Col be less than zero. * \see spFindElement() */ /* >>> Local variables: * pElement (RealNumber *) * Pointer to the element. * * >>> Possible errors: * spNO_MEMORY * Error is not cleared in this routine. */ spElement * spGetElement( spMatrix eMatrix, int Row, int Col ) { MatrixPtr Matrix = (MatrixPtr)eMatrix; ElementPtr pElement; /* Begin `spGetElement'. */ ASSERT_IS_SPARSE( Matrix ); vASSERT( Row >= 0 AND Col >= 0, "Negative row or column number" ); if ((Row == 0) OR (Col == 0)) return &Matrix->TrashCan.Real; #if NOT TRANSLATE vASSERT( NOT Matrix->Reordered, "Set TRANSLATE to add elements to a reordered matrix" ); #endif #if TRANSLATE Translate( Matrix, &Row, &Col ); if (Matrix->Error == spNO_MEMORY) return NULL; #endif #if NOT TRANSLATE #if NOT EXPANDABLE vASSERT( (Row <= Matrix->Size) AND (Col <= Matrix->Size), "Row or column number too large" ); #endif #if EXPANDABLE /* Re-size Matrix if necessary. */ if ((Row > Matrix->Size) OR (Col > Matrix->Size)) EnlargeMatrix( Matrix, MAX(Row, Col) ); if (Matrix->Error == spNO_MEMORY) return NULL; #endif #endif if ((Row != Col) OR ((pElement = Matrix->Diag[Row]) == NULL)) { /* * Element does not exist or does not reside along diagonal. Search * for element and if it does not exist, create it. */ pElement = spcCreateElement( Matrix, Row, Col, &(Matrix->FirstInRow[Row]), &(Matrix->FirstInCol[Col]), NO ); } /* * Cast pointer into a pointer to a RealNumber. This requires that Real * be the first record in the MatrixElement structure. */ return &pElement->Real; } #if TRANSLATE /* * TRANSLATE EXTERNAL INDICES TO INTERNAL * * Convert internal row and column numbers to internal row and column numbers. * Also updates Ext/Int maps. * * * >>> Arguments: * Matrix (MatrixPtr) * Pointer to the matrix. * Row (int *) * Upon entry Row is either a external row number of an external node * number. Upon entry, the internal equivalent is supplied. * Col (int *) * Upon entry Column is either a external column number of an external node * number. Upon entry, the internal equivalent is supplied. * * >>> Local variables: * ExtCol (int) * Temporary variable used to hold the external column or node number * during the external to internal column number translation. * ExtRow (int) * Temporary variable used to hold the external row or node number during * the external to internal row number translation. * IntCol (int) * Temporary variable used to hold the internal column or node number * during the external to internal column number translation. * IntRow (int) * Temporary variable used to hold the internal row or node number during * the external to internal row number translation. */ static void Translate( MatrixPtr Matrix, int *Row, int *Col ) { register int IntRow, IntCol, ExtRow, ExtCol; /* Begin `Translate'. */ ExtRow = *Row; ExtCol = *Col; /* Expand translation arrays if necessary. */ if ((ExtRow > Matrix->AllocatedExtSize) OR (ExtCol > Matrix->AllocatedExtSize)) { ExpandTranslationArrays( Matrix, MAX(ExtRow, ExtCol) ); if (Matrix->Error == spNO_MEMORY) return; } /* Set ExtSize if necessary. */ if ((ExtRow > Matrix->ExtSize) OR (ExtCol > Matrix->ExtSize)) Matrix->ExtSize = MAX(ExtRow, ExtCol); /* Translate external row or node number to internal row or node number. */ if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1) { Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize; Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize; IntRow = Matrix->CurrentSize; #if NOT EXPANDABLE vASSERT( IntRow <= Matrix->Size, "Matrix size fixed" ); #endif #if EXPANDABLE /* Re-size Matrix if necessary. */ if (IntRow > Matrix->Size) EnlargeMatrix( Matrix, IntRow ); if (Matrix->Error == spNO_MEMORY) return; #endif Matrix->IntToExtRowMap[IntRow] = ExtRow; Matrix->IntToExtColMap[IntRow] = ExtRow; } /* Translate external column or node number to internal column or node number.*/ if ((IntCol = Matrix->ExtToIntColMap[ExtCol]) == -1) { Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize; Matrix->ExtToIntColMap[ExtCol] = Matrix->CurrentSize; IntCol = Matrix->CurrentSize; #if NOT EXPANDABLE vASSERT( IntCol <= Matrix->Size, "Matrix size fixed" ); #endif #if EXPANDABLE /* Re-size Matrix if necessary. */ if (IntCol > Matrix->Size) EnlargeMatrix( Matrix, IntCol ); if (Matrix->Error == spNO_MEMORY) return; #endif Matrix->IntToExtRowMap[IntCol] = ExtCol; Matrix->IntToExtColMap[IntCol] = ExtCol; } *Row = IntRow; *Col = IntCol; return; } #endif #if QUAD_ELEMENT /*! * Performs same function as spGetElement() except rather than one * element, all four matrix elements for a floating two terminal * admittance component are added. This routine also works if component * is grounded. Positive elements are placed at [Node1,Node2] and * [Node2,Node1]. This routine is only to be used after spCreate() * and before spMNA_Preorder(), spFactor() or spOrderAndFactor(). * * \return * Error code. Possible errors include \a spNO_MEMORY. * Error is not cleared in this routine. * * \param Matrix * Pointer to the matrix that component is to be entered in. * \param Node1 * Row and column indices for elements. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Node zero is the * ground node. In no case may \a Node1 be less than zero. * \param Node2 * Row and column indices for elements. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Node zero is the * ground node. In no case may \a Node2 be less than zero. * \param Template * Collection of pointers to four elements that are later used to directly * address elements. User must supply the template, this routine will * fill it. */ spError spGetAdmittance( spMatrix Matrix, int Node1, int Node2, struct spTemplate *Template ) { /* Begin `spGetAdmittance'. */ Template->Element1 = spGetElement(Matrix, Node1, Node1 ); Template->Element2 = spGetElement(Matrix, Node2, Node2 ); Template->Element3Negated = spGetElement( Matrix, Node2, Node1 ); Template->Element4Negated = spGetElement( Matrix, Node1, Node2 ); if ( (Template->Element1 == NULL) OR (Template->Element2 == NULL) OR (Template->Element3Negated == NULL) OR (Template->Element4Negated == NULL) ) return spNO_MEMORY; if (Node1 == 0) SWAP( RealNumber*, Template->Element1, Template->Element2 ); return spOKAY; } #endif /* QUAD_ELEMENT */ #if QUAD_ELEMENT /*! * Similar to spGetAdmittance(), except that spGetAdmittance() only * handles 2-terminal components, whereas spGetQuad() handles simple * 4-terminals as well. These 4-terminals are simply generalized * 2-terminals with the option of having the sense terminals different * from the source and sink terminals. spGetQuad() adds four * elements to the matrix. Positive elements occur at [Row1,Col1] * [Row2,Col2] while negative elements occur at [Row1,Col2] and [Row2,Col1]. * The routine works fine if any of the rows and columns are zero. * This routine is only to be used after spCreate() and before * spMNA_Preorder(), spFactor() or spOrderAndFactor() * unless \a TRANSLATE is set true. * * \return * Error code. Possible errors include \a spNO_MEMORY. * Error is not cleared in this routine. * * \param Matrix * Pointer to the matrix that component is to be entered in. * \param Row1 * First row index for elements. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Zero is the * ground row. In no case may Row1 be less than zero. * \param Row2 * Second row index for elements. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Zero is the * ground row. In no case may Row2 be less than zero. * \param Col1 * First column index for elements. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Zero is the * ground column. In no case may Col1 be less than zero. * \param Col2 * Second column index for elements. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Zero is the * ground column. In no case may Col2 be less than zero. * \param Template * Collection of pointers to four elements that are later used to directly * address elements. User must supply the template, this routine will * fill it. */ spError spGetQuad( spMatrix Matrix, int Row1, int Row2, int Col1, int Col2, struct spTemplate *Template ) { /* Begin `spGetQuad'. */ Template->Element1 = spGetElement( Matrix, Row1, Col1); Template->Element2 = spGetElement( Matrix, Row2, Col2 ); Template->Element3Negated = spGetElement( Matrix, Row2, Col1 ); Template->Element4Negated = spGetElement( Matrix, Row1, Col2 ); if ( (Template->Element1 == NULL) OR (Template->Element2 == NULL) OR (Template->Element3Negated == NULL) OR (Template->Element4Negated == NULL) ) return spNO_MEMORY; if (Template->Element1 == &((MatrixPtr)Matrix)->TrashCan.Real) SWAP( RealNumber *, Template->Element1, Template->Element2 ); return spOKAY; } #endif /* QUAD_ELEMENT */ #if QUAD_ELEMENT /*! * Addition of four structural ones to matrix by index. * Performs similar function to spGetQuad() except this routine is * meant for components that do not have an admittance representation. * * The following stamp is used: \code * Pos Neg Eqn * Pos [ . . 1 ] * Neg [ . . -1 ] * Eqn [ 1 -1 . ] * \endcode * * \return * Error code. Possible errors include \a spNO_MEMORY. * Error is not cleared in this routine. * * \param Matrix * Pointer to the matrix that component is to be entered in. * \param Pos * See stamp above. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Zero is the * ground row. In no case may \a Pos be less than zero. * \param Neg * See stamp above. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Zero is the * ground row. In no case may \a Neg be less than zero. * \param Eqn * See stamp above. Must be in the range of [0..Size] * unless the options \a EXPANDABLE or \a TRANSLATE are used. Zero is the * ground row. In no case may \a Eqn be less than zero. * \param Template * Collection of pointers to four elements that are later used to directly * address elements. User must supply the template, this routine will * fill it. */ spError spGetOnes( spMatrix Matrix, int Pos, int Neg, int Eqn, struct spTemplate *Template ) { /* Begin `spGetOnes'. */ Template->Element4Negated = spGetElement( Matrix, Neg, Eqn ); Template->Element3Negated = spGetElement( Matrix, Eqn, Neg ); Template->Element2 = spGetElement( Matrix, Pos, Eqn ); Template->Element1 = spGetElement( Matrix, Eqn, Pos ); if ( (Template->Element1 == NULL) OR (Template->Element2 == NULL) OR (Template->Element3Negated == NULL) OR (Template->Element4Negated == NULL) ) return spNO_MEMORY; spADD_REAL_QUAD( *Template, 1.0 ); return spOKAY; } #endif /* QUAD_ELEMENT */ /* * FIND DIAGONAL * * This routine is used to find a diagonal element. It will not * create it if it does not exist. * * >>> Returned: * A pointer to the desired element, or NULL if it does not exist. * * >>> Arguments: * Matrix (MatrixPtr) * Pointer to matrix. * Index (int) * Row, Col index for diagonal element. * * >>> Local variables: * pElement (ElementPtr) * Pointer to an element in the matrix. */ ElementPtr spcFindDiag( MatrixPtr Matrix, register int Index ) { register ElementPtr pElement; /* Begin `spcFindDiag'. */ pElement = Matrix->FirstInCol[Index]; /* Search column for element. */ while ((pElement != NULL) AND (pElement->Row < Index)) pElement = pElement->NextInCol; if (pElement AND (pElement->Row == Index)) return pElement; else return NULL; } /* * CREATE AND SPLICE ELEMENT INTO MATRIX * * This routine is used to create new matrix elements and splice them into the * matrix. * * >>> Returned: * A pointer to the element that was created is returned. * * >>> Arguments: * Matrix (MatrixPtr) * Pointer to matrix. * Row (int) * Row index for element. * Col (int) * Column index for element. * ppToLeft (ElementPtr *) * This contains the address of the pointer to an element to the left * of the one being created. It is used to speed the search and if it * is immediately to the left, it is updated with address of the * created element. * ppAbove (ElementPtr *) * This contains the address of the pointer to an element above the * one being created. It is used to speed the search and it if it * is immediatley above, it is updated with address of the created * element. * Fillin (BOOLEAN) * Flag that indicates if created element is to be a fill-in. * * >>> Local variables: * pElement (ElementPtr) * Pointer to an element in the matrix. * pCreatedElement (ElementPtr) * Pointer to the desired element, the one that was just created. * * >>> Possible errors: * spNO_MEMORY */ ElementPtr spcCreateElement( MatrixPtr Matrix, int Row, register int Col, register ElementPtr *ppToLeft, register ElementPtr *ppAbove, BOOLEAN Fillin ) { register ElementPtr pElement, pCreatedElement; /* Begin `spcCreateElement'. */ /* Find element immediately above the desired element. */ pElement = *ppAbove; while ((pElement != NULL) AND (pElement->Row < Row)) { ppAbove = &pElement->NextInCol; pElement = *ppAbove; } if ((pElement != NULL) AND (pElement->Row == Row)) return pElement; /* The desired element does not exist, create it. */ if (Fillin) { pCreatedElement = spcGetFillin( Matrix ); Matrix->Fillins++; /* Update Markowitz counts and products. */ ++Matrix->MarkowitzRow[Row]; spcMarkoProd( Matrix->MarkowitzProd[Row], Matrix->MarkowitzRow[Row], Matrix->MarkowitzCol[Row] ); if ((Matrix->MarkowitzRow[Row] == 1) AND (Matrix->MarkowitzCol[Row] != 0)) { Matrix->Singletons--; } ++Matrix->MarkowitzCol[Col]; spcMarkoProd( Matrix->MarkowitzProd[Col], Matrix->MarkowitzCol[Col], Matrix->MarkowitzRow[Col] ); if ((Matrix->MarkowitzRow[Col] != 0) AND (Matrix->MarkowitzCol[Col] == 1)) { Matrix->Singletons--; } } else { pCreatedElement = spcGetElement( Matrix ); Matrix->NeedsOrdering = YES; } if (pCreatedElement == NULL) return NULL; Matrix->Elements++; /* Initialize Element. */ pCreatedElement->Row = Row; pCreatedElement->Col = Col; pCreatedElement->Real = 0.0; pCreatedElement->Imag = 0.0; pCreatedElement->pInitInfo = NULL; /* If element is on diagonal, store pointer in Diag. */ if (Row == Col) Matrix->Diag[Row] = pCreatedElement; /* Splice element into column. */ pCreatedElement->NextInCol = *ppAbove; *ppAbove = pCreatedElement; /* Find Element immediately to the left of the fill-in. */ if (Matrix->RowsLinked) { pElement = *ppToLeft; while (pElement != NULL) { if (pElement->Col < Col) { ppToLeft = &pElement->NextInRow; pElement = *ppToLeft; } else break; /* while loop */ } /* Splice element into row. */ pCreatedElement->NextInRow = *ppToLeft; *ppToLeft = pCreatedElement; } return pCreatedElement; } /* * * LINK ROWS * * This routine is used to generate the row links. The spGetElement() * routines do not create row links, which are needed by the spFactor() * routines. * * >>> Arguments: * Matrix (MatrixPtr) * Pointer to the matrix. * * >>> Local variables: * pElement (ElementPtr) * Pointer to an element in the matrix. * FirstInRowEntry (ElementPtr *) * A pointer into the FirstInRow array. Points to the FirstInRow entry * currently being operated upon. * FirstInRowArray (ArrayOfElementPtrs) * A pointer to the FirstInRow array. Same as Matrix->FirstInRow but * resides in a register and requires less indirection so is faster to * use. * Col (int) * Column currently being operated upon. */ void spcLinkRows( MatrixPtr Matrix ) { register ElementPtr pElement, *FirstInRowEntry; register ArrayOfElementPtrs FirstInRowArray; register int Col; /* Begin `spcLinkRows'. */ FirstInRowArray = Matrix->FirstInRow; for (Col = Matrix->Size; Col >= 1; Col--) FirstInRowArray[Col] = NULL; for (Col = Matrix->Size; Col >= 1; Col--) { /* Generate row links for the elements in the Col'th column. */ pElement = Matrix->FirstInCol[Col]; while (pElement != NULL) { pElement->Col = Col; FirstInRowEntry = &FirstInRowArray[pElement->Row]; pElement->NextInRow = *FirstInRowEntry; *FirstInRowEntry = pElement; pElement = pElement->NextInCol; } } Matrix->RowsLinked = YES; return; } /* * ENLARGE MATRIX * * Increases the size of the matrix. * * >>> Arguments: * Matrix (MatrixPtr) * Pointer to the matrix. * NewSize (int) * The new size of the matrix. * * >>> Local variables: * OldAllocatedSize (int) * The allocated size of the matrix before it is expanded. */ static void EnlargeMatrix( MatrixPtr Matrix, register int NewSize ) { register int I, OldAllocatedSize = Matrix->AllocatedSize; /* Begin `EnlargeMatrix'. */ Matrix->Size = NewSize; if (NewSize <= OldAllocatedSize) return; /* Expand the matrix frame. */ NewSize = MAX( NewSize, (int)(EXPANSION_FACTOR * OldAllocatedSize) ); Matrix->AllocatedSize = NewSize; if (( REALLOC(Matrix->IntToExtColMap, int, NewSize+1)) == NULL) { Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->IntToExtRowMap, int, NewSize+1)) == NULL) { Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->Diag, ElementPtr, NewSize+1)) == NULL) { Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->FirstInCol, ElementPtr, NewSize+1)) == NULL) { Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->FirstInRow, ElementPtr, NewSize+1)) == NULL) { Matrix->Error = spNO_MEMORY; return; } /* * Destroy the Markowitz and Intermediate vectors, they will be recreated * in spOrderAndFactor(). */ FREE( Matrix->MarkowitzRow ); FREE( Matrix->MarkowitzCol ); FREE( Matrix->MarkowitzProd ); FREE( Matrix->DoRealDirect ); FREE( Matrix->DoCmplxDirect ); FREE( Matrix->Intermediate ); Matrix->InternalVectorsAllocated = NO; /* Initialize the new portion of the vectors. */ for (I = OldAllocatedSize+1; I <= NewSize; I++) { Matrix->IntToExtColMap[I] = I; Matrix->IntToExtRowMap[I] = I; Matrix->Diag[I] = NULL; Matrix->FirstInRow[I] = NULL; Matrix->FirstInCol[I] = NULL; } return; } #if TRANSLATE /* * EXPAND TRANSLATION ARRAYS * * Increases the size arrays that are used to translate external to internal * row and column numbers. * * >>> Arguments: * Matrix (MatrixPtr) * Pointer to the matrix. * NewSize (int) * The new size of the translation arrays. * * >>> Local variables: * OldAllocatedSize (int) * The allocated size of the translation arrays before being expanded. */ static void ExpandTranslationArrays( MatrixPtr Matrix, register int NewSize ) { register int I, OldAllocatedSize = Matrix->AllocatedExtSize; /* Begin `ExpandTranslationArrays'. */ Matrix->ExtSize = NewSize; if (NewSize <= OldAllocatedSize) return; /* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */ NewSize = MAX( NewSize, (int)(EXPANSION_FACTOR * OldAllocatedSize) ); Matrix->AllocatedExtSize = NewSize; if (( REALLOC(Matrix->ExtToIntRowMap, int, NewSize+1)) == NULL) { Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->ExtToIntColMap, int, NewSize+1)) == NULL) { Matrix->Error = spNO_MEMORY; return; } /* Initialize the new portion of the vectors. */ for (I = OldAllocatedSize+1; I <= NewSize; I++) { Matrix->ExtToIntRowMap[I] = -1; Matrix->ExtToIntColMap[I] = -1; } return; } #endif /*! * Initialize the matrix. * * Sparse allows the user to keep initialization information with each * structurally nonzero matrix element. Each element has a pointer * that is set and used by the user. The user can set this pointer * using spInstallInitInfo() and may be read using spGetInitInfo(). Both * may be used only after the element exists. The function * spInitialize() is a user customizable way to initialize the matrix. * Passed to this routine is a function pointer. spInitialize() sweeps * through every element in the matrix and checks the \a pInitInfo * pointer (the user supplied pointer). If the \a pInitInfo is \a NULL, * which is true unless the user changes it (almost always true for * fill-ins), then the element is zeroed. Otherwise, the function * pointer is called and passed the \a pInitInfo pointer as well as the * element pointer and the external row and column numbers. If the * user sets the value of each element, then spInitialize() replaces * spClear(). * * The user function is expected to return a nonzero integer if there * is a fatal error and zero otherwise. Upon encountering a nonzero * return code, spInitialize() terminates, sets the error state of * the matrix to be \a spMANGLED, and returns the error code. * * \return * Returns the return value of the \a pInit() function. * \param eMatrix * Pointer to matrix. * \param pInit * Pointer to a function that initializes an element. * \see spClear() */ int spInitialize( spMatrix eMatrix, int (*pInit)( spElement *pElement, void * pInitInfo, int Row, int Col ) ) { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement; int J, Error, Col; /* Begin `spInitialize'. */ ASSERT_IS_SPARSE( Matrix ); /* Clear imaginary part of matrix if matrix is real but was complex. */ if (Matrix->PreviousMatrixWasComplex AND NOT Matrix->Complex) { for (J = Matrix->Size; J > 0; J--) { pElement = Matrix->FirstInCol[J]; while (pElement != NULL) { pElement->Imag = 0.0; pElement = pElement->NextInCol; } } } /* Initialize the matrix. */ for (J = Matrix->Size; J > 0; J--) { pElement = Matrix->FirstInCol[J]; Col = Matrix->IntToExtColMap[J]; while (pElement != NULL) { if (pElement->pInitInfo == NULL) { pElement->Real = 0.0; pElement->Imag = 0.0; } else { Error = (*pInit)((RealNumber *)pElement, pElement->pInitInfo, Matrix->IntToExtRowMap[pElement->Row], Col); if (Error) { Matrix->Error = spMANGLED; return Error; } } pElement = pElement->NextInCol; } } /* Empty the trash. */ Matrix->TrashCan.Real = 0.0; Matrix->TrashCan.Imag = 0.0; Matrix->Error = spOKAY; Matrix->Factored = NO; Matrix->SingularCol = 0; Matrix->SingularRow = 0; Matrix->PreviousMatrixWasComplex = Matrix->Complex; return 0; } /*! * This function installs a pointer to a data structure that is used * to contain initialization information to a matrix element. It is * is then used by spInitialize() to initialize the matrix. * * \param pElement * Pointer to matrix element. * \param pInitInfo * Pointer to the data structure that will contain initialiation * information. * \see spInitialize() */ void spInstallInitInfo( spElement *pElement, void * pInitInfo ) { /* Begin `spInstallInitInfo'. */ vASSERT( pElement != NULL, "Invalid element pointer" ); ((ElementPtr)pElement)->pInitInfo = pInitInfo; } /*! * This function returns a pointer to a data structure that is used * to contain initialization information to a matrix element. * * \return * The pointer to the initialiation information data structure * that is associated with a particular matrix element. * * \param pElement * Pointer to the matrix element. * * \see spInitialize() */ void * spGetInitInfo( spElement *pElement ) { /* Begin `spGetInitInfo'. */ vASSERT( pElement != NULL, "Invalid element pointer" ); return (void *)((ElementPtr)pElement)->pInitInfo; }