diff --git a/src/sparse13/spalloc.cpp b/src/sparse13/spalloc.cpp index ab85ebefd9..233093c61f 100644 --- a/src/sparse13/spalloc.cpp +++ b/src/sparse13/spalloc.cpp @@ -15,7 +15,6 @@ * spWhereSingular * spGetSize * spSetReal - * spSetComplex * spFillinCount * spElementCount * @@ -45,11 +44,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -126,12 +120,10 @@ char* spCreate(int Size, BOOLEAN Complex, int* pError) } /* Test for valid type. */ -#if NOT spCOMPLEX if (Complex) { *pError = spPANIC; return NULL; } -#endif #if NOT REAL if (NOT Complex) { *pError = spPANIC; @@ -193,9 +185,6 @@ char* spCreate(int Size, BOOLEAN Complex, int* pError) /* Take out the trash. */ Matrix->TrashCan.Real = 0.0; -#if spCOMPLEX - Matrix->TrashCan.Imag = 0.0; -#endif Matrix->TrashCan.Row = 0; Matrix->TrashCan.Col = 0; Matrix->TrashCan.NextInRow = NULL; @@ -233,24 +222,6 @@ char* spCreate(int Size, BOOLEAN Complex, int* pError) Matrix->IntToExtColMap[I] = I; } -#if TRANSLATE - /* Allocate space in memory for ExtToIntColMap vector. */ - if ((Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL) - goto MemoryError; - - /* Allocate space in memory for ExtToIntRowMap vector. */ - if ((Matrix->ExtToIntRowMap = ALLOC(int, SizePlusOne)) == NULL) - goto MemoryError; - - /* Initialize MapExtToInt vectors. */ - for (I = 1; I <= AllocatedSize; I++) { - Matrix->ExtToIntColMap[I] = -1; - Matrix->ExtToIntRowMap[I] = -1; - } - Matrix->ExtToIntColMap[0] = 0; - Matrix->ExtToIntRowMap[0] = 0; -#endif - /* Allocate space for fill-ins and initial set of elements. */ InitializeElementBlocks(Matrix, SPACE_FOR_ELEMENTS * AllocatedSize, SPACE_FOR_FILL_INS * AllocatedSize); @@ -671,14 +642,7 @@ int spGetSize(char* eMatrix, BOOLEAN External) /* Begin `spGetSize'. */ ASSERT(IS_SPARSE(Matrix)); -#if TRANSLATE - if (External) - return Matrix->ExtSize; - else - return Matrix->Size; -#else return Matrix->Size; -#endif } /* @@ -700,15 +664,6 @@ void spSetReal(char* eMatrix) return; } -void spSetComplex(char* eMatrix) -{ - /* Begin `spSetComplex'. */ - - ASSERT(IS_SPARSE((MatrixPtr)eMatrix) AND spCOMPLEX); - ((MatrixPtr)eMatrix)->Complex = YES; - return; -} - /* * ELEMENT OR FILL-IN COUNT * diff --git a/src/sparse13/spbuild.cpp b/src/sparse13/spbuild.cpp index 3b7e9a1da2..f0a5867609 100644 --- a/src/sparse13/spbuild.cpp +++ b/src/sparse13/spbuild.cpp @@ -20,11 +20,9 @@ * * >>> Other functions contained in this file: * spcFindElementInCol - * Translate * spcCreateElement * spcLinkRows * EnlargeMatrix - * ExpandTranslationArrays */ /* @@ -43,11 +41,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -70,9 +63,7 @@ static char RCSid[] = "@(#)$Header$"; extern ElementPtr spcGetFillin(MatrixPtr Matrix); extern ElementPtr spcGetElement(MatrixPtr Matrix); -static void Translate(MatrixPtr Matrix, int* Row, int* Col); static void EnlargeMatrix(MatrixPtr Matrix, int NewSize); -static void ExpandTranslationArrays(MatrixPtr Matrix, int NewSize); /* * CLEAR MATRIX @@ -98,18 +89,6 @@ void spClear(char* eMatrix) ASSERT(IS_SPARSE(Matrix)); /* Clear matrix. */ -#if spCOMPLEX - 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 -#endif { for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInCol[I]; @@ -122,9 +101,6 @@ void spClear(char* eMatrix) /* Empty the trash. */ Matrix->TrashCan.Real = 0.0; -#if spCOMPLEX - Matrix->TrashCan.Imag = 0.0; -#endif Matrix->Error = spOKAY; Matrix->Factored = NO; @@ -183,17 +159,8 @@ RealNumber* spGetElement(char* eMatrix, int Row, int Col) if ((Row == 0) OR(Col == 0)) return &Matrix->TrashCan.Real; -#if NOT TRANSLATE ASSERT(Matrix->NeedsOrdering); -#endif - -#if TRANSLATE - Translate(Matrix, &Row, &Col); - if (Matrix->Error == spNO_MEMORY) - return NULL; -#endif -#if NOT TRANSLATE #if NOT EXPANDABLE ASSERT(Row <= Matrix->Size AND Col <= Matrix->Size); #endif @@ -204,7 +171,6 @@ RealNumber* spGetElement(char* eMatrix, int Row, int Col) EnlargeMatrix(Matrix, MAX(Row, Col)); if (Matrix->Error == spNO_MEMORY) return NULL; -#endif #endif /* @@ -290,109 +256,6 @@ ElementPtr spcFindElementInCol(MatrixPtr Matrix, ElementPtr* LastAddr, int Row, return NULL; } -#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) -{ - 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 - ASSERT(IntRow <= Matrix->Size); -#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 - ASSERT(IntCol <= Matrix->Size); -#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 /* * ADDITION OF ADMITTANCE TO MATRIX BY INDEX @@ -649,9 +512,6 @@ ElementPtr spcCreateElement(MatrixPtr Matrix, int Row, int Col, ElementPtr* Last pElement->Row = Row; pElement->Col = Col; pElement->Real = 0.0; -#if spCOMPLEX - pElement->Imag = 0.0; -#endif #if INITIALIZE pElement->pInitInfo = NULL; #endif @@ -708,9 +568,6 @@ ElementPtr spcCreateElement(MatrixPtr Matrix, int Row, int Col, ElementPtr* Last pElement->Col = Col; #endif pElement->Real = 0.0; -#if spCOMPLEX - pElement->Imag = 0.0; -#endif #if INITIALIZE pElement->pInitInfo = NULL; #endif @@ -847,58 +704,6 @@ static void EnlargeMatrix(MatrixPtr Matrix, int NewSize) 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, int NewSize) -{ - int I, OldAllocatedSize = Matrix->AllocatedExtSize; - - /* Begin `ExpandTranslationArrays'. */ - Matrix->ExtSize = NewSize; - - if (NewSize <= OldAllocatedSize) - return; - - /* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */ - NewSize = MAX(NewSize, 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 - #if INITIALIZE /* * INITIALIZE MATRIX @@ -960,18 +765,6 @@ int (*pInit)(); /* Begin `spInitialize'. */ ASSERT(IS_SPARSE(Matrix)); -#if spCOMPLEX - /* 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; - } - } - } -#endif /* spCOMPLEX */ /* Initialize the matrix. */ for (J = Matrix->Size; J > 0; J--) { @@ -980,9 +773,6 @@ int (*pInit)(); while (pElement != NULL) { if (pElement->pInitInfo == NULL) { pElement->Real = 0.0; -#if spCOMPLEX - pElement->Imag = 0.0; -#endif } else { Error = (*pInit)((RealNumber*)pElement, pElement->pInitInfo, Matrix->IntToExtRowMap[pElement->Row], Col); @@ -997,9 +787,6 @@ int (*pInit)(); /* Empty the trash. */ Matrix->TrashCan.Real = 0.0; -#if spCOMPLEX - Matrix->TrashCan.Imag = 0.0; -#endif Matrix->Error = spOKAY; Matrix->Factored = NO; diff --git a/src/sparse13/spconfig.h b/src/sparse13/spconfig.h index ae52e16216..22cd3f9dcf 100644 --- a/src/sparse13/spconfig.h +++ b/src/sparse13/spconfig.h @@ -67,9 +67,6 @@ * both real and complex systems at the same time, but there is a * slight speed and memory advantage if the routines are complied * to handle only real systems of equations. - * spCOMPLEX - * This specifies that the routines will be complied to handle - * complex systems of equations. * EXPANDABLE * Setting this compiler flag true (1) makes the matrix * expandable before it has been factored. If the matrix is @@ -140,16 +137,6 @@ * must have an allocated length of one plus the size of the * matrix. ARRAY_OFFSET must be either 0 or 1, no other offsets * are valid. - * spSEPARATED_COMPLEX_VECTORS - * This specifies the format for complex vectors. If this is set - * false then a complex vector is made up of one double sized - * array of RealNumber's in which the real and imaginary numbers - * are placed in the alternately array in the array. In other - * words, the first entry would be Complex[1].Real, then comes - * Complex[1].Imag, then Complex[1].Real, etc. If - * spSEPARATED_COMPLEX_VECTORS is set true, then each complex - * vector is represented by two arrays of RealNumbers, one with - * the real terms, the other with the imaginary. [NO] * MODIFIED_MARKOWITZ * This specifies that the modified Markowitz method of pivot * selection is to be used. The modified Markowitz method differs @@ -243,12 +230,6 @@ /* Begin options. */ #define REAL YES #define EXPANDABLE YES -#if defined(cmplx_spPrefix) -/* NEURON's nonlinz.cpp uses cmplx_spGetElement after previous use of matrix */ -#define TRANSLATE YES -#else -#define TRANSLATE NO /* instead of YES */ -#endif #define INITIALIZE NO /* instead of YES */ #define DIAGONAL_PIVOTING YES #define ARRAY_OFFSET NOT FORTRAN @@ -275,13 +256,6 @@ * with user code, so use 0 for NO and 1 for YES. */ #endif /* spINSIDE_SPARSE */ -#if defined(cmplx_spPrefix) -#define spCOMPLEX 1 -#define spSEPARATED_COMPLEX_VECTORS 1 -#else -#define spCOMPLEX 0 /* instead of 1 */ -#define spSEPARATED_COMPLEX_VECTORS 0 -#endif #ifdef spINSIDE_SPARSE /* diff --git a/src/sparse13/spdefs.h b/src/sparse13/spdefs.h index f7977c77ad..78ae57ea8b 100644 --- a/src/sparse13/spdefs.h +++ b/src/sparse13/spdefs.h @@ -41,7 +41,6 @@ #ifdef lint #undef REAL -#undef spCOMPLEX #undef EXPANDABLE #undef TRANSLATE #undef INITIALIZE @@ -60,7 +59,6 @@ #undef DEBUG #define REAL YES -#define spCOMPLEX YES #define EXPANDABLE YES #define TRANSLATE YES #define INITIALIZE YES @@ -133,225 +131,7 @@ } /* Macro function that returns the approx absolute value of a complex number. */ -#if spCOMPLEX -#define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag)) -#else #define ELEMENT_MAG(ptr) ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real) -#endif - -/* Complex assignment statements. */ -#define CMPLX_ASSIGN(to, from) \ - { \ - (to).Real = (from).Real; \ - (to).Imag = (from).Imag; \ - } -#define CMPLX_CONJ_ASSIGN(to, from) \ - { \ - (to).Real = (from).Real; \ - (to).Imag = -(from).Imag; \ - } -#define CMPLX_NEGATE_ASSIGN(to, from) \ - { \ - (to).Real = -(from).Real; \ - (to).Imag = -(from).Imag; \ - } -#define CMPLX_CONJ_NEGATE_ASSIGN(to, from) \ - { \ - (to).Real = -(from).Real; \ - (to).Imag = (from).Imag; \ - } -#define CMPLX_CONJ(a) (a).Imag = -(a).Imag -#define CMPLX_NEGATE(a) \ - { \ - (a).Real = -(a).Real; \ - (a).Imag = -(a).Imag; \ - } - -/* Macro that returns the approx magnitude (L-1 norm) of a complex number. */ -#define CMPLX_1_NORM(a) (ABS((a).Real) + ABS((a).Imag)) - -/* Macro that returns the approx magnitude (L-infinity norm) of a complex. */ -#define CMPLX_INF_NORM(a) (MAX(ABS((a).Real), ABS((a).Imag))) - -/* Macro function that returns the magnitude (L-2 norm) of a complex number. */ -#define CMPLX_2_NORM(a) (sqrt((a).Real * (a).Real + (a).Imag * (a).Imag)) - -/* Macro function that performs complex addition. */ -#define CMPLX_ADD(to, from_a, from_b) \ - { \ - (to).Real = (from_a).Real + (from_b).Real; \ - (to).Imag = (from_a).Imag + (from_b).Imag; \ - } - -/* Macro function that performs complex subtraction. */ -#define CMPLX_SUBT(to, from_a, from_b) \ - { \ - (to).Real = (from_a).Real - (from_b).Real; \ - (to).Imag = (from_a).Imag - (from_b).Imag; \ - } - -/* Macro function that is equivalent to += operator for complex numbers. */ -#define CMPLX_ADD_ASSIGN(to, from) \ - { \ - (to).Real += (from).Real; \ - (to).Imag += (from).Imag; \ - } - -/* Macro function that is equivalent to -= operator for complex numbers. */ -#define CMPLX_SUBT_ASSIGN(to, from) \ - { \ - (to).Real -= (from).Real; \ - (to).Imag -= (from).Imag; \ - } - -/* Macro function that multiplies a complex number by a scalar. */ -#define SCLR_MULT(to, sclr, cmplx) \ - { \ - (to).Real = (sclr) * (cmplx).Real; \ - (to).Imag = (sclr) * (cmplx).Imag; \ - } - -/* Macro function that multiply-assigns a complex number by a scalar. */ -#define SCLR_MULT_ASSIGN(to, sclr) \ - { \ - (to).Real *= (sclr); \ - (to).Imag *= (sclr); \ - } - -/* Macro function that multiplies two complex numbers. */ -#define CMPLX_MULT(to, from_a, from_b) \ - { \ - (to).Real = (from_a).Real * (from_b).Real - (from_a).Imag * (from_b).Imag; \ - (to).Imag = (from_a).Real * (from_b).Imag + (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that implements to *= from for complex numbers. */ -#define CMPLX_MULT_ASSIGN(to, from) \ - { \ - RealNumber to_real_ = (to).Real; \ - (to).Real = to_real_ * (from).Real - (to).Imag * (from).Imag; \ - (to).Imag = to_real_ * (from).Imag + (to).Imag * (from).Real; \ - } - -/* Macro function that multiplies two complex numbers, the first of which is - * conjugated. */ -#define CMPLX_CONJ_MULT(to, from_a, from_b) \ - { \ - (to).Real = (from_a).Real * (from_b).Real + (from_a).Imag * (from_b).Imag; \ - (to).Imag = (from_a).Real * (from_b).Imag - (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then adds them - * to another. to = add + mult_a * mult_b */ -#define CMPLX_MULT_ADD(to, mult_a, mult_b, add) \ - { \ - (to).Real = (mult_a).Real * (mult_b).Real - (mult_a).Imag * (mult_b).Imag + (add).Real; \ - (to).Imag = (mult_a).Real * (mult_b).Imag + (mult_a).Imag * (mult_b).Real + (add).Imag; \ - } - -/* Macro function that subtracts the product of two complex numbers from - * another. to = subt - mult_a * mult_b */ -#define CMPLX_MULT_SUBT(to, mult_a, mult_b, subt) \ - { \ - (to).Real = (subt).Real - (mult_a).Real * (mult_b).Real + (mult_a).Imag * (mult_b).Imag; \ - (to).Imag = (subt).Imag - (mult_a).Real * (mult_b).Imag - (mult_a).Imag * (mult_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then adds them - * to another. to = add + mult_a* * mult_b where mult_a* represents mult_a - * conjugate. */ -#define CMPLX_CONJ_MULT_ADD(to, mult_a, mult_b, add) \ - { \ - (to).Real = (mult_a).Real * (mult_b).Real + (mult_a).Imag * (mult_b).Imag + (add).Real; \ - (to).Imag = (mult_a).Real * (mult_b).Imag - (mult_a).Imag * (mult_b).Real + (add).Imag; \ - } - -/* Macro function that multiplies two complex numbers and then adds them - * to another. to += mult_a * mult_b */ -#define CMPLX_MULT_ADD_ASSIGN(to, from_a, from_b) \ - { \ - (to).Real += (from_a).Real * (from_b).Real - (from_a).Imag * (from_b).Imag; \ - (to).Imag += (from_a).Real * (from_b).Imag + (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then subtracts them - * from another. */ -#define CMPLX_MULT_SUBT_ASSIGN(to, from_a, from_b) \ - { \ - (to).Real -= (from_a).Real * (from_b).Real - (from_a).Imag * (from_b).Imag; \ - (to).Imag -= (from_a).Real * (from_b).Imag + (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then adds them - * to the destination. to += from_a* * from_b where from_a* represents from_a - * conjugate. */ -#define CMPLX_CONJ_MULT_ADD_ASSIGN(to, from_a, from_b) \ - { \ - (to).Real += (from_a).Real * (from_b).Real + (from_a).Imag * (from_b).Imag; \ - (to).Imag += (from_a).Real * (from_b).Imag - (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then subtracts them - * from the destination. to -= from_a* * from_b where from_a* represents from_a - * conjugate. */ -#define CMPLX_CONJ_MULT_SUBT_ASSIGN(to, from_a, from_b) \ - { \ - (to).Real -= (from_a).Real * (from_b).Real + (from_a).Imag * (from_b).Imag; \ - (to).Imag -= (from_a).Real * (from_b).Imag - (from_a).Imag * (from_b).Real; \ - } - -/* - * Macro functions that provide complex division. - */ - -/* Complex division: to = num / den */ -#define CMPLX_DIV(to, num, den) \ - { \ - RealNumber r_, s_; \ - if (((den).Real >= (den).Imag AND(den).Real > -(den).Imag) OR((den).Real < (den).Imag AND(den).Real <= -(den).Imag)) { \ - r_ = (den).Imag / (den).Real; \ - s_ = (den).Real + r_ * (den).Imag; \ - (to).Real = ((num).Real + r_ * (num).Imag) / s_; \ - (to).Imag = ((num).Imag - r_ * (num).Real) / s_; \ - } else { \ - r_ = (den).Real / (den).Imag; \ - s_ = (den).Imag + r_ * (den).Real; \ - (to).Real = (r_ * (num).Real + (num).Imag) / s_; \ - (to).Imag = (r_ * (num).Imag - (num).Real) / s_; \ - } \ - } - -/* Complex division and assignment: num /= den */ -#define CMPLX_DIV_ASSIGN(num, den) \ - { \ - RealNumber r_, s_, t_; \ - if (((den).Real >= (den).Imag AND(den).Real > -(den).Imag) OR((den).Real < (den).Imag AND(den).Real <= -(den).Imag)) { \ - r_ = (den).Imag / (den).Real; \ - s_ = (den).Real + r_ * (den).Imag; \ - t_ = ((num).Real + r_ * (num).Imag) / s_; \ - (num).Imag = ((num).Imag - r_ * (num).Real) / s_; \ - (num).Real = t_; \ - } else { \ - r_ = (den).Real / (den).Imag; \ - s_ = (den).Imag + r_ * (den).Real; \ - t_ = (r_ * (num).Real + (num).Imag) / s_; \ - (num).Imag = (r_ * (num).Imag - (num).Real) / s_; \ - (num).Real = t_; \ - } \ - } - -/* Complex reciprocation: to = 1.0 / den */ -#define CMPLX_RECIPROCAL(to, den) \ - { \ - RealNumber r_; \ - if (((den).Real >= (den).Imag AND(den).Real > -(den).Imag) OR((den).Real < (den).Imag AND(den).Real <= -(den).Imag)) { \ - r_ = (den).Imag / (den).Real; \ - (to).Imag = -r_ * ((to).Real = 1.0 / ((den).Real + r_ * (den).Imag)); \ - } else { \ - r_ = (den).Real / (den).Imag; \ - (to).Real = -r_ * ((to).Imag = -1.0 / ((den).Imag + r_ * (den).Real)); \ - } \ - } /* * ASSERT and ABORT @@ -421,26 +201,6 @@ typedef spREAL RealNumber, *RealVector; -/* - * COMPLEX NUMBER DATA STRUCTURE - * - * >>> Structure fields: - * Real (RealNumber) - * The real portion of the number. Real must be the first - * field in this structure. - * Imag (RealNumber) - * The imaginary portion of the number. This field must follow - * immediately after Real. - */ - -/* Begin `ComplexNumber'. */ - -typedef struct -{ - RealNumber Real; - RealNumber Imag; -} ComplexNumber, *ComplexVector; - /* * MATRIX ELEMENT DATA STRUCTURE * @@ -455,11 +215,6 @@ typedef struct * Real (RealNumber) * The real portion of the value of the element. Real must be the first * field in this structure. - * Imag (RealNumber) - * The imaginary portion of the value of the element. If the matrix - * routines are not compiled to handle complex matrices, then this - * field does not exist. If it exists, it must follow immediately after - * Real. * Row (int) * The row number of the element. * Col (int) @@ -488,9 +243,6 @@ typedef struct struct MatrixElement { RealNumber Real; -#if spCOMPLEX - RealNumber Imag; -#endif int Row; int Col; struct MatrixElement* NextInRow; diff --git a/src/sparse13/spfactor.cpp b/src/sparse13/spfactor.cpp index c3f960f6f4..f01c612deb 100644 --- a/src/sparse13/spfactor.cpp +++ b/src/sparse13/spfactor.cpp @@ -13,7 +13,6 @@ * spPartition * * >>> Other functions contained in this file: - * FactorComplexMatrix CreateInternalVectors * CountMarkowitz MarkowitzProducts * SearchForPivot SearchForSingleton * QuicklySearchDiagonal SearchDiagonal @@ -21,7 +20,7 @@ * FindBiggestInColExclude ExchangeRowsAndCols * spcRowExchange spcColExchange * ExchangeColElements ExchangeRowElements - * RealRowColElimination ComplexRowColElimination + * RealRowColElimination * UpdateMarkowitzNumbers CreateFillin * MatrixIsSingular ZeroPivot * WriteStatus @@ -43,11 +42,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -67,7 +61,6 @@ static char RCSid[] = "@(#)$Header$"; #include /* avoid "declared implicitly `extern' and later `static' " warnings. */ -static int FactorComplexMatrix(MatrixPtr Matrix); static void CreateInternalVectors(MatrixPtr Matrix); static void CountMarkowitz(MatrixPtr Matrix, RealVector RHS, int Step); static void MarkowitzProducts(MatrixPtr Matrix, int Step); @@ -82,7 +75,6 @@ static void ExchangeRowsAndCols(MatrixPtr Matrix, ElementPtr pPivot, int Step); static void ExchangeColElements(MatrixPtr Matrix, int Row1, ElementPtr Element1, int Row2, ElementPtr Element2, int Column); static void ExchangeRowElements(MatrixPtr Matrix, int Col1, ElementPtr Element1, int Col2, ElementPtr Element2, int Row); static void RealRowColElimination(MatrixPtr Matrix, ElementPtr pPivot); -static void ComplexRowColElimination(MatrixPtr Matrix, ElementPtr pPivot); static void UpdateMarkowitzNumbers(MatrixPtr Matrix, ElementPtr pPivot); static ElementPtr CreateFillin(MatrixPtr Matrix, int Row, int Col); static int MatrixIsSingular(MatrixPtr Matrix, int Step); @@ -215,10 +207,7 @@ int spOrderAndFactor(char* eMatrix, RealNumber* RHS, RealNumber RelThreshold, Re pPivot = Matrix->Diag[Step]; LargestInCol = FindLargestInCol(pPivot->NextInCol); if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot))) { - if (Matrix->Complex) - ComplexRowColElimination(Matrix, pPivot); - else - RealRowColElimination(Matrix, pPivot); + RealRowColElimination(Matrix, pPivot); } else { ReorderingRequired = YES; break; /* for loop */ @@ -265,10 +254,7 @@ int spOrderAndFactor(char* eMatrix, RealNumber* RHS, RealNumber RelThreshold, Re return MatrixIsSingular(Matrix, Step); ExchangeRowsAndCols(Matrix, pPivot, Step); - if (Matrix->Complex) - ComplexRowColElimination(Matrix, pPivot); - else - RealRowColElimination(Matrix, pPivot); + RealRowColElimination(Matrix, pPivot); if (Matrix->Error >= spFATAL) return Matrix->Error; @@ -335,10 +321,6 @@ int spFactor(char* eMatrix) } if (NOT Matrix->Partitioned) spPartition(eMatrix, spDEFAULT_PARTITION); -#if spCOMPLEX - if (Matrix->Complex) - return FactorComplexMatrix(Matrix); -#endif #if REAL Size = Matrix->Size; @@ -412,117 +394,6 @@ int spFactor(char* eMatrix) #endif /* REAL */ } -#if spCOMPLEX -/* - * FACTOR COMPLEX MATRIX - * - * This routine is the companion routine to spFactor(), it - * handles complex matrices. It is otherwise identical. - * - * >>> Returned: - * The error code is returned. Possible errors are listed below. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to matrix. - * - * >>> Possible errors: - * spSINGULAR - * Error is cleared in this function. - */ - -static int FactorComplexMatrix(MatrixPtr Matrix) -{ - ElementPtr pElement; - ElementPtr pColumn; - int Step, Size; - ComplexNumber Mult, Pivot; - - /* Begin `FactorComplexMatrix'. */ - ASSERT(Matrix->Complex); - - Size = Matrix->Size; - pElement = Matrix->Diag[1]; - if (ELEMENT_MAG(pElement) == 0.0) - return ZeroPivot(Matrix, 1); - /* Cmplx expr: *pPivot = 1.0 / *pPivot. */ - CMPLX_RECIPROCAL(*pElement, *pElement); - - /* Start factorization. */ - for (Step = 2; Step <= Size; Step++) { - if (Matrix->DoCmplxDirect[Step]) { /* Update column using direct addressing scatter-gather. */ - ComplexNumber* Dest; - Dest = (ComplexNumber*)Matrix->Intermediate; - - /* Scatter. */ - pElement = Matrix->FirstInCol[Step]; - while (pElement != NULL) { - Dest[pElement->Row] = *(ComplexNumber*)pElement; - pElement = pElement->NextInCol; - } - - /* Update column. */ - pColumn = Matrix->FirstInCol[Step]; - while (pColumn->Row < Step) { - pElement = Matrix->Diag[pColumn->Row]; - /* Cmplx expr: Mult = Dest[pColumn->Row] * (1.0 / *pPivot). */ - CMPLX_MULT(Mult, Dest[pColumn->Row], *pElement); - CMPLX_ASSIGN(*pColumn, Mult); - while ((pElement = pElement->NextInCol) != NULL) { /* Cmplx expr: Dest[pElement->Row] -= Mult * pElement */ - CMPLX_MULT_SUBT_ASSIGN(Dest[pElement->Row], Mult, *pElement); - } - pColumn = pColumn->NextInCol; - } - - /* Gather. */ - pElement = Matrix->Diag[Step]->NextInCol; - while (pElement != NULL) { - *(ComplexNumber*)pElement = Dest[pElement->Row]; - pElement = pElement->NextInCol; - } - - /* Check for singular matrix. */ - Pivot = Dest[Step]; - if (CMPLX_1_NORM(Pivot) == 0.0) - return ZeroPivot(Matrix, Step); - CMPLX_RECIPROCAL(*Matrix->Diag[Step], Pivot); - } else { /* Update column using direct addressing scatter-gather. */ - ComplexNumber** pDest; - pDest = (ComplexNumber**)Matrix->Intermediate; - - /* Scatter. */ - pElement = Matrix->FirstInCol[Step]; - while (pElement != NULL) { - pDest[pElement->Row] = (ComplexNumber*)pElement; - pElement = pElement->NextInCol; - } - - /* Update column. */ - pColumn = Matrix->FirstInCol[Step]; - while (pColumn->Row < Step) { - pElement = Matrix->Diag[pColumn->Row]; - /* Cmplx expr: Mult = *pDest[pColumn->Row] * (1.0 / *pPivot). */ - CMPLX_MULT(Mult, *pDest[pColumn->Row], *pElement); - CMPLX_ASSIGN(*pDest[pColumn->Row], Mult); - while ((pElement = pElement->NextInCol) != NULL) { /* Cmplx expr: *pDest[pElement->Row] -= Mult * pElement */ - CMPLX_MULT_SUBT_ASSIGN(*pDest[pElement->Row], Mult, *pElement); - } - pColumn = pColumn->NextInCol; - } - - /* Check for singular matrix. */ - pElement = Matrix->Diag[Step]; - if (ELEMENT_MAG(pElement) == 0.0) - return ZeroPivot(Matrix, Step); - CMPLX_RECIPROCAL(*pElement, *pElement); - } - } - - Matrix->Factored = YES; - return (Matrix->Error = spOKAY); -} -#endif /* spCOMPLEX */ - /* * PARTITION MATRIX * @@ -566,7 +437,7 @@ void spPartition(char* eMatrix, int Mode) ElementPtr pElement, pColumn; int Step, Size; int *Nc, *No, *Nm; - BOOLEAN *DoRealDirect, *DoCmplxDirect; + BOOLEAN *DoRealDirect; /* Begin `spPartition'. */ ASSERT(IS_SPARSE(Matrix)); @@ -574,7 +445,6 @@ void spPartition(char* eMatrix, int Mode) return; Size = Matrix->Size; DoRealDirect = Matrix->DoRealDirect; - DoCmplxDirect = Matrix->DoCmplxDirect; Matrix->Partitioned = YES; /* If partition is specified by the user, this is easy. */ @@ -584,18 +454,12 @@ void spPartition(char* eMatrix, int Mode) for (Step = 1; Step <= Size; Step++) #if REAL DoRealDirect[Step] = YES; -#endif -#if spCOMPLEX - DoCmplxDirect[Step] = YES; #endif return; } else if (Mode == spINDIRECT_PARTITION) { for (Step = 1; Step <= Size; Step++) #if REAL DoRealDirect[Step] = NO; -#endif -#if spCOMPLEX - DoCmplxDirect[Step] = NO; #endif return; } else @@ -641,9 +505,6 @@ void spPartition(char* eMatrix, int Mode) #if REAL DoRealDirect[Step] = (Nm[Step] + No[Step] > 3 * Nc[Step] - 2 * Nm[Step]); -#endif -#if spCOMPLEX - DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7 * Nc[Step] - 4 * Nm[Step]); #endif } @@ -703,25 +564,12 @@ static void CreateInternalVectors(MatrixPtr Matrix) Matrix->Error = spNO_MEMORY; } #endif -#if spCOMPLEX - if (Matrix->DoCmplxDirect == NULL) { - if ((Matrix->DoCmplxDirect = ALLOC(BOOLEAN, Size + 1)) == NULL) - Matrix->Error = spNO_MEMORY; - } -#endif /* Create Intermediate vectors for use in MatrixSolve. */ -#if spCOMPLEX - if (Matrix->Intermediate == NULL) { - if ((Matrix->Intermediate = ALLOC(RealNumber, 2 * (Size + 1))) == NULL) - Matrix->Error = spNO_MEMORY; - } -#else if (Matrix->Intermediate == NULL) { if ((Matrix->Intermediate = ALLOC(RealNumber, Size + 1)) == NULL) Matrix->Error = spNO_MEMORY; } -#endif if (Matrix->Error != spNO_MEMORY) Matrix->InternalVectorsAllocated = YES; @@ -765,17 +613,8 @@ static void CountMarkowitz(MatrixPtr Matrix, RealVector RHS, int Step) /* Correct array pointer for ARRAY_OFFSET. */ #if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX if (RHS != NULL) --RHS; -#else - if (RHS != NULL) { - if (Matrix->Complex) - RHS -= 2; - else - --RHS; - } -#endif #endif /* Generate MarkowitzRow Count for each row. */ @@ -793,19 +632,9 @@ static void CountMarkowitz(MatrixPtr Matrix, RealVector RHS, int Step) /* Include nonzero elements in the RHS vector. */ ExtRow = Matrix->IntToExtRowMap[I]; -#if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX if (RHS != NULL) if (RHS[ExtRow] != 0.0) Count++; -#else - if (RHS != NULL) { - if (Matrix->Complex) { - if ((RHS[2 * ExtRow] != 0.0) OR(RHS[2 * ExtRow + 1] != 0.0)) - Count++; - } else if (RHS[I] != 0.0) - Count++; - } -#endif Matrix->MarkowitzRow[I] = Count; } @@ -2015,10 +1844,6 @@ void spcRowExchange(MatrixPtr Matrix, int Row1, int Row2) SWAP(int, Matrix->MarkowitzRow[Row1], Matrix->MarkowitzRow[Row2]); SWAP(ElementPtr, Matrix->FirstInRow[Row1], Matrix->FirstInRow[Row2]); SWAP(int, Matrix->IntToExtRowMap[Row1], Matrix->IntToExtRowMap[Row2]); -#if TRANSLATE - Matrix->ExtToIntRowMap[Matrix->IntToExtRowMap[Row1]] = Row1; - Matrix->ExtToIntRowMap[Matrix->IntToExtRowMap[Row2]] = Row2; -#endif return; } @@ -2102,10 +1927,6 @@ void spcColExchange(MatrixPtr Matrix, int Col1, int Col2) SWAP(int, Matrix->MarkowitzCol[Col1], Matrix->MarkowitzCol[Col2]); SWAP(ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]); SWAP(int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]); -#if TRANSLATE - Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col1]] = Col1; - Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col2]] = Col2; -#endif return; } @@ -2454,57 +2275,6 @@ static void RealRowColElimination(MatrixPtr Matrix, ElementPtr pPivot) * spNO_MEMORY */ -static void ComplexRowColElimination(MatrixPtr Matrix, ElementPtr pPivot) -{ -#if spCOMPLEX - ElementPtr pSub; - int Row; - ElementPtr pLower, pUpper; - - /* Begin `ComplexRowColElimination'. */ - - /* Test for zero pivot. */ - if (ELEMENT_MAG(pPivot) == 0.0) { - (void)MatrixIsSingular(Matrix, pPivot->Row); - return; - } - CMPLX_RECIPROCAL(*pPivot, *pPivot); - - pUpper = pPivot->NextInRow; - while (pUpper != NULL) { - /* Calculate upper triangular element. */ - /* Cmplx expr: *pUpper = *pUpper * (1.0 / *pPivot). */ - CMPLX_MULT_ASSIGN(*pUpper, *pPivot); - - pSub = pUpper->NextInCol; - pLower = pPivot->NextInCol; - while (pLower != NULL) { - Row = pLower->Row; - - /* Find element in row that lines up with current lower triangular element. */ - while (pSub != NULL AND pSub->Row < Row) - pSub = pSub->NextInCol; - - /* Test to see if desired element was not found, if not, create fill-in. */ - if (pSub == NULL OR pSub->Row > Row) { - pSub = CreateFillin(Matrix, Row, pUpper->Col); - if (pSub == NULL) { - Matrix->Error = spNO_MEMORY; - return; - } - } - - /* Cmplx expr: pElement -= *pUpper * pLower. */ - CMPLX_MULT_SUBT_ASSIGN(*pSub, *pUpper, *pLower); - pSub = pSub->NextInCol; - pLower = pLower->NextInCol; - } - pUpper = pUpper->NextInRow; - } - return; -#endif /* spCOMPLEX */ -} - /* * UPDATE MARKOWITZ NUMBERS * diff --git a/src/sparse13/spmatrix.h b/src/sparse13/spmatrix.h index a923809ef8..ce5878c324 100644 --- a/src/sparse13/spmatrix.h +++ b/src/sparse13/spmatrix.h @@ -227,11 +227,6 @@ struct spTemplate { */ /* Begin function declarations. */ - -#if defined(__STDC__) || defined(__cplusplus) - -/* For compilers that understand function prototypes. */ - extern void spClear(char*); extern spREAL spCondition(char*, spREAL, int*); extern char* spCreate(int, int, int*); @@ -260,7 +255,6 @@ extern void spPrint(char*, int, int, int); extern spREAL spPseudoCondition(char*); extern spREAL spRoundoff(char*, spREAL); extern void spScale(char*, spREAL[], spREAL[]); -extern void spSetComplex(char*); extern void spSetReal(char*); extern void spStripFills(char*); extern void spWhereSingular(char*, int*, int*); @@ -273,48 +267,4 @@ extern void spMultTransposed(char* eMatrix, spREAL* RHS, spREAL* Solution, std:: extern void spSolve(char* eMatrix, spREAL* RHS, spREAL* Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); extern void spSolveTransposed(char*, spREAL*, spREAL*, std::optional = std::nullopt, std::optional = std::nullopt); -#else /* NOT defined(__STDC__) */ - -/* For compilers that do not understand function prototypes. */ - -extern void spClear(); -extern spREAL spCondition(); -extern char* spCreate(); -extern void spDeleteRowAndCol(); -extern void spDestroy(); -extern void spDeterminant(); -extern int spElementCount(); -extern int spError(); -extern int spFactor(); -extern int spFileMatrix(); -extern int spFileStats(); -extern int spFileVector(); -extern int spFillinCount(); -extern int spGetAdmittance(); -extern spREAL* spGetElement(char*, int, int); -extern char* spGetInitInfo(); -extern int spGetOnes(); -extern int spGetQuad(); -extern int spGetSize(); -extern int spInitialize(); -extern void spInstallInitInfo(); -extern spREAL spLargestElement(); -extern void spMNA_Preorder(); -extern void spMultiply(); -extern void spMultTransposed(); -extern spREAL spNorm(); -extern int spOrderAndFactor(); -extern void spPartition(); -extern void spPrint(); -extern spREAL spPseudoCondition(); -extern spREAL spRoundoff(); -extern void spScale(); -extern void spSetComplex(); -extern void spSetReal(); -extern void spSolve(); -extern void spSolveTransposed(); -extern void spStripFills(); -extern void spWhereSingular(); -#endif /* defined(__STDC__) */ - #endif /* spOKAY */ diff --git a/src/sparse13/spoutput.cpp b/src/sparse13/spoutput.cpp index 534025bcd2..5c33082f3e 100644 --- a/src/sparse13/spoutput.cpp +++ b/src/sparse13/spoutput.cpp @@ -33,11 +33,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "$Header$"; -#endif - /* * IMPORTS * @@ -143,11 +138,7 @@ void spPrint(char* eMatrix, int PrintReordered, int Data, int Header) Size = Matrix->Size; /* Create a packed external to internal row and column translation array. */ -#if TRANSLATE - Top = Matrix->AllocatedExtSize; -#else Top = Matrix->AllocatedSize; -#endif CALLOC(PrintOrdToIntRowMap, int, Top + 1); CALLOC(PrintOrdToIntColMap, int, Top + 1); if (PrintOrdToIntRowMap == NULL OR PrintOrdToIntColMap == NULL) { @@ -286,19 +277,6 @@ void spPrint(char* eMatrix, int PrintReordered, int Data, int Header) } printf("\n"); -#if spCOMPLEX - if (Matrix->Complex AND Data) { - printf(" "); - for (J = StartCol; J <= StopCol; J++) { - if (pImagElements[J - StartCol] != NULL) { - printf(" %8.2lgj", - (double)pImagElements[J - StartCol]->Imag); - } else - printf(" "); - } - printf("\n"); - } -#endif /* spCOMPLEX */ } /* Calculate index of first column in next group. */ @@ -440,31 +418,6 @@ int spFileMatrix(char* eMatrix, char* File, char* Label, int Reordered, int Data return 0; } -#if spCOMPLEX - if (Data AND Matrix->Complex) { - for (I = 1; I <= Size; I++) { - pElement = Matrix->FirstInCol[I]; - while (pElement != NULL) { - if (Reordered) { - Row = pElement->Row; - Col = I; - } else { - Row = Matrix->IntToExtRowMap[pElement->Row]; - Col = Matrix->IntToExtColMap[I]; - } - Err = fprintf(pMatrixFile, "%d\t%d\t%-.15lg\t%-.15lg\n", - Row, Col, (double)pElement->Real, (double)pElement->Imag); - if (Err < 0) - return 0; - pElement = pElement->NextInCol; - } - } - /* Output terminator, a line of zeros. */ - if (Header) - if (fprintf(pMatrixFile, "0\t0\t0.0\t0.0\n") < 0) - return 0; - } -#endif /* spCOMPLEX */ #if REAL if (Data AND NOT Matrix->Complex) { @@ -511,11 +464,9 @@ int spFileMatrix(char* eMatrix, char* File, char* Label, int Reordered, int Data * File (char *) * Name of file into which matrix is to be written. * RHS (RealNumber []) - * Right-hand side vector. This is only the real portion if - * spSEPARATED_COMPLEX_VECTORS is true. + * Right-hand side vector. * iRHS (RealNumber []) - * Right-hand side vector, imaginary portion. Not necessary if matrix - * is real or if spSEPARATED_COMPLEX_VECTORS is set false. + * Right-hand side vector, imaginary portion. * * >>> Local variables: * pMatrixFile (FILE *) @@ -523,17 +474,12 @@ int spFileMatrix(char* eMatrix, char* File, char* Label, int Reordered, int Data * Size (int) * The size of the matrix. * - * >>> Obscure Macros - * IMAG_RHS - * Replaces itself with `, iRHS' if the options spCOMPLEX and - * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears - * without a trace. */ int spFileVector(char* eMatrix, char* File, RealVector RHS, RealVector iRHS) { MatrixPtr Matrix = (MatrixPtr)eMatrix; - int I, Size, Err; + int I, Size; FILE* pMatrixFile; /* Begin `spFileVector'. */ @@ -545,44 +491,11 @@ int spFileVector(char* eMatrix, char* File, RealVector RHS, RealVector iRHS) /* Correct array pointers for ARRAY_OFFSET. */ #if NOT ARRAY_OFFSET -#if spCOMPLEX - if (Matrix->Complex) { -#if spSEPARATED_COMPLEX_VECTORS - ASSERT(iRHS != NULL) - --RHS; - --iRHS; -#else - RHS -= 2; -#endif - } else -#endif /* spCOMPLEX */ --RHS; #endif /* NOT ARRAY_OFFSET */ /* Output vector. */ Size = Matrix->Size; -#if spCOMPLEX - if (Matrix->Complex) { -#if spSEPARATED_COMPLEX_VECTORS - for (I = 1; I <= Size; I++) { - Err = fprintf(pMatrixFile, "%-.15lg\t%-.15lg\n", - (double)RHS[I], (double)iRHS[I]); - if (Err < 0) - return 0; - } -#else - for (I = 1; I <= Size; I++) { - Err = fprintf(pMatrixFile, "%-.15lg\t%-.15lg\n", - (double)RHS[2 * I], (double)RHS[2 * I + 1]); - if (Err < 0) - return 0; - } -#endif - } -#endif /* spCOMPLEX */ -#if REAL AND spCOMPLEX - else -#endif #if REAL { for (I = 1; I <= Size; I++) { @@ -656,10 +569,7 @@ int spFileStats(char* eMatrix, char* File, char* Label) fprintf(pStatsFile, "Matrix has not been factored.\n"); fprintf(pStatsFile, "||| Starting new matrix |||\n"); fprintf(pStatsFile, "%s\n", Label); - if (Matrix->Complex) - fprintf(pStatsFile, "Matrix is complex.\n"); - else - fprintf(pStatsFile, "Matrix is real.\n"); + fprintf(pStatsFile, "Matrix is real.\n"); fprintf(pStatsFile, " Size = %d\n", Size); /* Search matrix. */ diff --git a/src/sparse13/spsolve.cpp b/src/sparse13/spsolve.cpp index ee70d3b769..95bc4d436c 100644 --- a/src/sparse13/spsolve.cpp +++ b/src/sparse13/spsolve.cpp @@ -11,10 +11,6 @@ * >>> User accessible functions contained in this file: * spSolve * spSolveTransposed - * - * >>> Other functions contained in this file: - * SolveComplexMatrix - * SolveComplexTransposedMatrix */ /* @@ -33,11 +29,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -55,10 +46,6 @@ static char RCSid[] = "@(#)$Header$"; #include "spdefs.h" #include "spmatrix.h" -/* avoid "declared implicitly `extern' and later `static' " warnings. */ -static void SolveComplexMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); -static void SolveComplexTransposedMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); - /* * SOLVE MATRIX EQUATION * @@ -126,12 +113,6 @@ void spSolve(char* eMatrix, RealVector RHS, RealVector Solution, std::optionalComplex) { - SolveComplexMatrix(Matrix, RHS, Solution, iRHS, iSolution); - return; - } -#endif #if REAL Intermediate = Matrix->Intermediate; @@ -183,155 +164,6 @@ void spSolve(char* eMatrix, RealVector RHS, RealVector Solution, std::optional>> Arguments: - * Matrix (char *) - * Pointer to matrix. - * RHS (RealVector) - * RHS is the real portion of the input data array, the right hand - * side. This data is undisturbed and may be reused for other solves. - * Solution (RealVector) - * Solution is the real portion of the output data array. This routine - * is constructed such that RHS and Solution can be the same - * array. - * iRHS (RealVector) - * iRHS is the imaginary portion of the input data array, the right - * hand side. This data is undisturbed and may be reused for other solves. - * If spSEPARATED_COMPLEX_VECTOR is set false, there is no need to - * supply this array. - * iSolution (RealVector) - * iSolution is the imaginary portion of the output data array. This - * routine is constructed such that iRHS and iSolution can be - * the same array. If spSEPARATED_COMPLEX_VECTOR is set false, there is no - * need to supply this array. - * - * >>> Local variables: - * Intermediate (ComplexVector) - * Temporary storage for use in forward elimination and backward - * substitution. Commonly referred to as c, when the LU factorization - * equations are given as Ax = b, Lc = b, Ux = c. - * Local version of Matrix->Intermediate, which was created during - * the initial factorization in function CreateInternalVectors() in the - * matrix factorization module. - * pElement (ElementPtr) - * Pointer used to address elements in both the lower and upper triangle - * matrices. - * pExtOrder (int *) - * Pointer used to sequentially access each entry in IntToExtRowMap - * and IntToExtColMap arrays. Used to quickly scramble and unscramble - * RHS and Solution to account for row and column interchanges. - * pPivot (ElementPtr) - * Pointer that points to current pivot or diagonal element. - * Size (int) - * Size of matrix. Made local to reduce indirection. - * Temp (ComplexNumber) - * Temporary storage for entries in arrays. - */ - -static void SolveComplexMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS, std::optional iSolution) -{ - ElementPtr pElement; - ComplexVector Intermediate; - int I, *pExtOrder, Size; - ElementPtr pPivot; - ComplexVector ExtVector; - ComplexNumber Temp; - - /* Begin `SolveComplexMatrix'. */ - - Size = Matrix->Size; - Intermediate = (ComplexVector)Matrix->Intermediate; - -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS - --RHS; - --iRHS; - --Solution; - --iSolution; -#else - RHS -= 2; - Solution -= 2; -#endif -#endif - - /* Initialize Intermediate vector. */ - pExtOrder = &Matrix->IntToExtRowMap[Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Size; I > 0; I--) { - Intermediate[I].Real = RHS[*(pExtOrder)]; - ASSERT(iRHS.has_value()); - Intermediate[I].Imag = iRHS.value()[*(pExtOrder--)]; - } -#else - ExtVector = (ComplexVector)RHS; - for (I = Size; I > 0; I--) - Intermediate[I] = ExtVector[*(pExtOrder--)]; -#endif - - /* Forward substitution. Solves Lc = b.*/ - for (I = 1; I <= Size; I++) { - Temp = Intermediate[I]; - - /* This step of the substitution is skipped if Temp equals zero. */ - if ((Temp.Real != 0.0) OR(Temp.Imag != 0.0)) { - pPivot = Matrix->Diag[I]; - /* Cmplx expr: Temp *= (1.0 / Pivot). */ - CMPLX_MULT_ASSIGN(Temp, *pPivot); - Intermediate[I] = Temp; - pElement = pPivot->NextInCol; - while (pElement != NULL) { - /* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */ - CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row], - Temp, *pElement); - pElement = pElement->NextInCol; - } - } - } - - /* Backward Substitution. Solves Ux = c.*/ - for (I = Size; I > 0; I--) { - Temp = Intermediate[I]; - pElement = Matrix->Diag[I]->NextInRow; - - while (pElement != NULL) { - /* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */ - CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement, Intermediate[pElement->Col]); - pElement = pElement->NextInRow; - } - Intermediate[I] = Temp; - } - - /* Unscramble Intermediate vector while placing data in to Solution vector. */ - pExtOrder = &Matrix->IntToExtColMap[Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Size; I > 0; I--) { - Solution[*(pExtOrder)] = Intermediate[I].Real; - iSolution.value()[*(pExtOrder--)] = Intermediate[I].Imag; - } -#else - ExtVector = (ComplexVector)Solution; - for (I = Size; I > 0; I--) - ExtVector[*(pExtOrder--)] = Intermediate[I]; -#endif - - return; -} -#endif /* spCOMPLEX */ - #if TRANSPOSE /* * SOLVE TRANSPOSED MATRIX EQUATION @@ -402,13 +234,6 @@ void spSolveTransposed(char* eMatrix, RealVector RHS, RealVector Solution, std:: /* Begin `spSolveTransposed'. */ ASSERT(IS_VALID(Matrix) AND IS_FACTORED(Matrix)); -#if spCOMPLEX - if (Matrix->Complex) { - SolveComplexTransposedMatrix(Matrix, RHS, Solution, iRHS, iSolution); - return; - } -#endif - #if REAL Size = Matrix->Size; Intermediate = Matrix->Intermediate; @@ -457,156 +282,3 @@ void spSolveTransposed(char* eMatrix, RealVector RHS, RealVector Solution, std:: #endif /* REAL */ } #endif /* TRANSPOSE */ - -#if TRANSPOSE AND spCOMPLEX -/* - * SOLVE COMPLEX TRANSPOSED MATRIX EQUATION - * - * Performs forward elimination and back substitution to find the - * unknown vector from the RHS vector and transposed factored - * matrix. This routine is useful when performing sensitivity analysis - * on a circuit using the adjoint method. This routine assumes that - * the pivots are associated with the untransposed lower triangular - * (L) matrix and that the diagonal of the untransposed upper - * triangular (U) matrix consists of ones. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to matrix. - * RHS (RealVector) - * RHS is the input data array, the right hand - * side. This data is undisturbed and may be reused for other solves. - * This vector is only the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set true. - * Solution (RealVector) - * Solution is the real portion of the output data array. This routine - * is constructed such that RHS and Solution can be the same array. - * This vector is only the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set true. - * iRHS (RealVector) - * iRHS is the imaginary portion of the input data array, the right - * hand side. This data is undisturbed and may be reused for other solves. - * If either spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set false, there - * is no need to supply this array. - * iSolution (RealVector) - * iSolution is the imaginary portion of the output data array. This - * routine is constructed such that iRHS and iSolution can be - * the same array. If spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set - * false, there is no need to supply this array. - * - * >>> Local variables: - * Intermediate (ComplexVector) - * Temporary storage for use in forward elimination and backward - * substitution. Commonly referred to as c, when the LU factorization - * equations are given as Ax = b, Lc = b, Ux = c. Local version of - * Matrix->Intermediate, which was created during - * the initial factorization in function CreateInternalVectors() in the - * matrix factorization module. - * pElement (ElementPtr) - * Pointer used to address elements in both the lower and upper triangle - * matrices. - * pExtOrder (int *) - * Pointer used to sequentially access each entry in IntToExtRowMap - * and IntToExtColMap arrays. Used to quickly scramble and unscramble - * RHS and Solution to account for row and column interchanges. - * pPivot (ElementPtr) - * Pointer that points to current pivot or diagonal element. - * Size (int) - * Size of matrix. Made local to reduce indirection. - * Temp (ComplexNumber) - * Temporary storage for entries in arrays. - * - */ - -static void -SolveComplexTransposedMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS, std::optional iSolution) -{ - ElementPtr pElement; - ComplexVector Intermediate; - int I, *pExtOrder, Size; - ComplexVector ExtVector; - ElementPtr pPivot; - ComplexNumber Temp; - - /* Begin `SolveComplexTransposedMatrix'. */ - - Size = Matrix->Size; - Intermediate = (ComplexVector)Matrix->Intermediate; - -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS - --RHS; - --iRHS; - --Solution; - --iSolution; -#else - RHS -= 2; - Solution -= 2; -#endif -#endif - - /* Initialize Intermediate vector. */ - pExtOrder = &Matrix->IntToExtColMap[Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Size; I > 0; I--) { - Intermediate[I].Real = RHS[*(pExtOrder)]; - ASSERT(iRHS.has_value()); - Intermediate[I].Imag = iRHS.value()[*(pExtOrder--)]; - } -#else - ExtVector = (ComplexVector)RHS; - for (I = Size; I > 0; I--) - Intermediate[I] = ExtVector[*(pExtOrder--)]; -#endif - - /* Forward elimination. */ - for (I = 1; I <= Size; I++) { - Temp = Intermediate[I]; - - /* This step of the elimination is skipped if Temp equals zero. */ - if ((Temp.Real != 0.0) OR(Temp.Imag != 0.0)) { - pElement = Matrix->Diag[I]->NextInRow; - while (pElement != NULL) { - /* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */ - CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Col], - Temp, *pElement); - pElement = pElement->NextInRow; - } - } - } - - /* Backward Substitution. */ - for (I = Size; I > 0; I--) { - pPivot = Matrix->Diag[I]; - Temp = Intermediate[I]; - pElement = pPivot->NextInCol; - - while (pElement != NULL) { - /* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */ - CMPLX_MULT_SUBT_ASSIGN(Temp, Intermediate[pElement->Row], *pElement); - - pElement = pElement->NextInCol; - } - /* Cmplx expr: Intermediate = Temp * (1.0 / *pPivot). */ - CMPLX_MULT(Intermediate[I], Temp, *pPivot); - } - - /* Unscramble Intermediate vector while placing data in to Solution vector. */ - pExtOrder = &Matrix->IntToExtRowMap[Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Size; I > 0; I--) { - Solution[*(pExtOrder)] = Intermediate[I].Real; - iSolution.value()[*(pExtOrder--)] = Intermediate[I].Imag; - } -#else - ExtVector = (ComplexVector)Solution; - for (I = Size; I > 0; I--) - ExtVector[*(pExtOrder--)] = Intermediate[I]; -#endif - - return; -} -#endif /* TRANSPOSE AND spCOMPLEX */ diff --git a/src/sparse13/sputils.cpp b/src/sparse13/sputils.cpp index 66984fef77..613d10216f 100644 --- a/src/sparse13/sputils.cpp +++ b/src/sparse13/sputils.cpp @@ -24,9 +24,6 @@ * >>> Other functions contained in this file: * CountTwins * SwapCols - * ScaleComplexMatrix - * ComplexMatrixMultiply - * ComplexCondition */ /* @@ -45,11 +42,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -76,10 +68,6 @@ extern ElementPtr spcFindElementInCol(MatrixPtr Matrix, ElementPtr* LastAddr, in /* avoid "declared implicitly `extern' and later `static' " warnings. */ static int CountTwins(MatrixPtr Matrix, int Col, ElementPtr* ppTwin1, ElementPtr* ppTwin2); static void SwapCols(MatrixPtr Matrix, ElementPtr pTwin1, ElementPtr pTwin2); -static void ScaleComplexMatrix(MatrixPtr Matrix, RealVector RHS_ScaleFactors, RealVector SolutionScaleFactors); -static void ComplexMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); -static void ComplexTransposedMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); -static RealNumber ComplexCondition(MatrixPtr Matrix, RealNumber NormOfMatrix, int* pError); #if MODIFIED_NODAL /* @@ -267,10 +255,6 @@ static void SwapCols(MatrixPtr Matrix, ElementPtr pTwin1, ElementPtr pTwin2) SWAP(ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]); SWAP(int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]); -#if TRANSLATE - Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col2]] = Col2; - Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col1]] = Col1; -#endif Matrix->Diag[Col1] = pTwin2; Matrix->Diag[Col2] = pTwin1; @@ -348,13 +332,6 @@ void spScale(char* eMatrix, RealVector RHS_ScaleFactors, RealVector SolutionScal if (NOT Matrix->RowsLinked) spcLinkRows(Matrix); -#if spCOMPLEX - if (Matrix->Complex) { - ScaleComplexMatrix(Matrix, RHS_ScaleFactors, SolutionScaleFactors); - return; - } -#endif - #if REAL lSize = Matrix->Size; @@ -395,108 +372,6 @@ void spScale(char* eMatrix, RealVector RHS_ScaleFactors, RealVector SolutionScal } #endif /* SCALING */ -#if spCOMPLEX AND SCALING -/* - * SCALE COMPLEX MATRIX - * - * This function scales the matrix to enhance the possibility of - * finding a good pivoting order. Note that scaling enhances accuracy - * of the solution only if it affects the pivoting order, so it makes - * no sense to scale the matrix before spFactor(). If scaling is - * desired it should be done before spOrderAndFactor(). There - * are several things to take into account when choosing the scale - * factors. First, the scale factors are directly multiplied against - * the elements in the matrix. To prevent roundoff, each scale factor - * should be equal to an integer power of the number base of the - * machine. Since most machines operate in base two, scale factors - * should be a power of two. Second, the matrix should be scaled such - * that the matrix of element uncertainties is equilibrated. Third, - * this function multiplies the scale factors by the elements, so if - * one row tends to have uncertainties 1000 times smaller than the - * other rows, then its scale factor should be 1024, not 1/1024. - * Fourth, to save time, this function does not scale rows or columns - * if their scale factors are equal to one. Thus, the scale factors - * should be normalized to the most common scale factor. Rows and - * columns should be normalized separately. For example, if the size - * of the matrix is 100 and 10 rows tend to have uncertainties near - * 1e-6 and the remaining 90 have uncertainties near 1e-12, then the - * scale factor for the 10 should be 1/1,048,576 and the scale factors - * for the remaining 90 should be 1. Fifth, since this routine - * directly operates on the matrix, it is necessary to apply the scale - * factors to the RHS and Solution vectors. It may be easier to - * simply use spOrderAndFactor() on a scaled matrix to choose the - * pivoting order, and then throw away the matrix. Subsequent - * factorizations, performed with spFactor(), will not need to have - * the RHS and Solution vectors descaled. Lastly, this function - * should not be executed before the function spMNA_Preorder. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to the matrix to be scaled. - * SolutionScaleFactors (RealVector) - * The array of Solution scale factors. These factors scale the columns. - * All scale factors are real valued. - * RHS_ScaleFactors (RealVector) - * The array of RHS scale factors. These factors scale the rows. - * All scale factors are real valued. - * - * >>> Local variables: - * lSize (int) - * Local version of the size of the matrix. - * pElement (ElementPtr) - * Pointer to an element in the matrix. - * pExtOrder (int *) - * Pointer into either IntToExtRowMap or IntToExtColMap vector. Used to - * compensate for any row or column swaps that have been performed. - * ScaleFactor (RealNumber) - * The scale factor being used on the current row or column. - */ - -static void ScaleComplexMatrix(MatrixPtr Matrix, RealVector RHS_ScaleFactors, RealVector SolutionScaleFactors) -{ - ElementPtr pElement; - int I, lSize, *pExtOrder; - RealNumber ScaleFactor; - - /* Begin `ScaleComplexMatrix'. */ - lSize = Matrix->Size; - -/* Correct pointers to arrays for ARRAY_OFFSET */ -#if NOT ARRAY_OFFSET - --RHS_ScaleFactors; - --SolutionScaleFactors; -#endif - - /* Scale Rows */ - pExtOrder = &Matrix->IntToExtRowMap[1]; - for (I = 1; I <= lSize; I++) { - if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0) { - pElement = Matrix->FirstInRow[I]; - - while (pElement != NULL) { - pElement->Real *= ScaleFactor; - pElement->Imag *= ScaleFactor; - pElement = pElement->NextInRow; - } - } - } - - /* Scale Columns */ - pExtOrder = &Matrix->IntToExtColMap[1]; - for (I = 1; I <= lSize; I++) { - if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0) { - pElement = Matrix->FirstInCol[I]; - - while (pElement != NULL) { - pElement->Real *= ScaleFactor; - pElement->Imag *= ScaleFactor; - pElement = pElement->NextInCol; - } - } - } - return; -} -#endif /* SCALING AND spCOMPLEX */ #if MULTIPLICATION /* @@ -516,12 +391,10 @@ static void ScaleComplexMatrix(MatrixPtr Matrix, RealVector RHS_ScaleFactors, Re * Solution is the vector being multiplied by the matrix. * iRHS (RealVector) * iRHS is the imaginary portion of the right hand side. This is - * what is being solved for. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. + * what is being solved for. * iSolution (RealVector) * iSolution is the imaginary portion of the vector being multiplied - * by the matrix. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. + * by the matrix. * */ @@ -538,13 +411,6 @@ void spMultiply(char* eMatrix, RealVector RHS, RealVector Solution, std::optiona if (NOT Matrix->RowsLinked) spcLinkRows(Matrix); -#if spCOMPLEX - if (Matrix->Complex) { - ComplexMatrixMultiply(Matrix, RHS, Solution, iRHS, iSolution); - return; - } -#endif - #if REAL #if NOT ARRAY_OFFSET /* Correct array pointers for ARRAY_OFFSET. */ @@ -574,95 +440,6 @@ void spMultiply(char* eMatrix, RealVector RHS, RealVector Solution, std::optiona } #endif /* MULTIPLICATION */ -#if spCOMPLEX AND MULTIPLICATION -/* - * COMPLEX MATRIX MULTIPLICATION - * - * Multiplies matrix by solution vector to find source vector. - * Assumes matrix has not been factored. This routine can be used - * as a test to see if solutions are correct. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to the matrix. - * RHS (RealVector) - * RHS is the right hand side. This is what is being solved for. - * This is only the real portion of the right-hand side if the matrix - * is complex and spSEPARATED_COMPLEX_VECTORS is set true. - * Solution (RealVector) - * Solution is the vector being multiplied by the matrix. This is only - * the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set true. - * iRHS (RealVector) - * iRHS is the imaginary portion of the right hand side. This is - * what is being solved for. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. - * iSolution (RealVector) - * iSolution is the imaginary portion of the vector being multiplied - * by the matrix. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. - * - */ - -static void ComplexMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS, std::optional iSolution) -{ - ElementPtr pElement; - ComplexVector Vector; - ComplexNumber Sum; - int I, *pExtOrder; - -/* Begin `ComplexMatrixMultiply'. */ - -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS - --RHS; - --iRHS; - --Solution; - --iSolution; -#else - RHS -= 2; - Solution -= 2; -#endif -#endif - - /* Initialize Intermediate vector with reordered Solution vector. */ - Vector = (ComplexVector)Matrix->Intermediate; - pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Matrix->Size; I > 0; I--) { - Vector[I].Real = Solution[*pExtOrder]; - ASSERT(iSolution.has_value()); - Vector[I].Imag = iSolution.value()[*(pExtOrder--)]; - } -#else - for (I = Matrix->Size; I > 0; I--) - Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)]; -#endif - - pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInRow[I]; - Sum.Real = Sum.Imag = 0.0; - - while (pElement != NULL) { /* Cmplx expression : Sum += Element * Vector[Col] */ - CMPLX_MULT_ADD_ASSIGN(Sum, *pElement, Vector[pElement->Col]); - pElement = pElement->NextInRow; - } - -#if spSEPARATED_COMPLEX_VECTORS - RHS[*pExtOrder] = Sum.Real; - ASSERT(iRHS != std::nullopt); - iRHS.value()[*pExtOrder--] = Sum.Imag; -#else - ((ComplexVector)RHS)[*pExtOrder--] = Sum; -#endif - } - return; -} -#endif /* spCOMPLEX AND MULTIPLICATION */ - #if MULTIPLICATION AND TRANSPOSE /* * TRANSPOSED MATRIX MULTIPLICATION @@ -681,12 +458,10 @@ static void ComplexMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector S * Solution is the vector being multiplied by the matrix. * iRHS (RealVector) * iRHS is the imaginary portion of the right hand side. This is - * what is being solved for. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. + * what is being solved for. * iSolution (RealVector) * iSolution is the imaginary portion of the vector being multiplied - * by the matrix. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. + * by the matrix. * */ @@ -701,12 +476,6 @@ void spMultTransposed(char* eMatrix, RealVector RHS, RealVector Solution, std::o /* Begin `spMultTransposed'. */ ASSERT(IS_SPARSE(Matrix) AND NOT Matrix->Factored); -#if spCOMPLEX - if (Matrix->Complex) { - ComplexTransposedMatrixMultiply(Matrix, RHS, Solution, iRHS, iSolution); - return; - } -#endif #if REAL #if NOT ARRAY_OFFSET @@ -737,94 +506,6 @@ void spMultTransposed(char* eMatrix, RealVector RHS, RealVector Solution, std::o } #endif /* MULTIPLICATION AND TRANSPOSE */ -#if spCOMPLEX AND MULTIPLICATION AND TRANSPOSE -/* - * COMPLEX TRANSPOSED MATRIX MULTIPLICATION - * - * Multiplies transposed matrix by solution vector to find source vector. - * Assumes matrix has not been factored. This routine can be used - * as a test to see if solutions are correct. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to the matrix. - * RHS (RealVector) - * RHS is the right hand side. This is what is being solved for. - * This is only the real portion of the right-hand side if the matrix - * is complex and spSEPARATED_COMPLEX_VECTORS is set true. - * Solution (RealVector) - * Solution is the vector being multiplied by the matrix. This is only - * the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set true. - * iRHS (RealVector) - * iRHS is the imaginary portion of the right hand side. This is - * what is being solved for. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. - * iSolution (RealVector) - * iSolution is the imaginary portion of the vector being multiplied - * by the matrix. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. - * - */ - -static void ComplexTransposedMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS, std::optional iSolution) -{ - ElementPtr pElement; - ComplexVector Vector; - ComplexNumber Sum; - int I, *pExtOrder; - -/* Begin `ComplexMatrixMultiply'. */ - -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS - --RHS; - --iRHS; - --Solution; - --iSolution; -#else - RHS -= 2; - Solution -= 2; -#endif -#endif - - /* Initialize Intermediate vector with reordered Solution vector. */ - Vector = (ComplexVector)Matrix->Intermediate; - pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Matrix->Size; I > 0; I--) { - Vector[I].Real = Solution[*pExtOrder]; - ASSERT(iSolution.has_value()) - Vector[I].Imag = iSolution.value()[*(pExtOrder--)]; - } -#else - for (I = Matrix->Size; I > 0; I--) - Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)]; -#endif - - pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInCol[I]; - Sum.Real = Sum.Imag = 0.0; - - while (pElement != NULL) { /* Cmplx expression : Sum += Element * Vector[Row] */ - CMPLX_MULT_ADD_ASSIGN(Sum, *pElement, Vector[pElement->Row]); - pElement = pElement->NextInCol; - } - -#if spSEPARATED_COMPLEX_VECTORS - RHS[*pExtOrder] = Sum.Real; - ASSERT(iRHS != std::nullopt); - iRHS.value()[*pExtOrder--] = Sum.Imag; -#else - ((ComplexVector)RHS)[*pExtOrder--] = Sum; -#endif - } - return; -} -#endif /* spCOMPLEX AND MULTIPLICATION AND TRANSPOSE */ #if DETERMINANT /* @@ -859,8 +540,6 @@ static void ComplexTransposedMatrixMultiply(MatrixPtr Matrix, RealVector RHS, Re * number is scaled to be greater than or equal to 1.0 and less than 10.0. * * >>> Local variables: - * Norm (RealNumber) - * L-infinity norm of a complex number. * Size (int) * Local storage for Matrix->Size. * Temp (RealNumber) @@ -871,10 +550,6 @@ void spDeterminant(char* eMatrix, int* pExponent, RealNumber* pDeterminant, std: { MatrixPtr Matrix = (MatrixPtr)eMatrix; int I, Size; - RealNumber Norm, nr, ni; - ComplexNumber Pivot, cDeterminant; - -#define NORM(a) (nr = ABS((a).Real), ni = ABS((a).Imag), MAX(nr, ni)) /* Begin `spDeterminant'. */ ASSERT(IS_SPARSE(Matrix) AND IS_FACTORED(Matrix)); @@ -882,72 +557,12 @@ void spDeterminant(char* eMatrix, int* pExponent, RealNumber* pDeterminant, std: if (Matrix->Error == spSINGULAR) { *pDeterminant = 0.0; -#if spCOMPLEX - if (Matrix->Complex) { - ASSERT(piDeterminant != std::nullopt); - *piDeterminant.value() = 0.0; - } -#endif return; } Size = Matrix->Size; I = 0; -#if spCOMPLEX - if (Matrix->Complex) /* Complex Case. */ - { - cDeterminant.Real = 1.0; - cDeterminant.Imag = 0.0; - - while (++I <= Size) { - CMPLX_RECIPROCAL(Pivot, *Matrix->Diag[I]); - CMPLX_MULT_ASSIGN(cDeterminant, Pivot); - - /* Scale Determinant. */ - Norm = NORM(cDeterminant); - if (Norm != 0.0) { - while (Norm >= 1.0e12) { - cDeterminant.Real *= 1.0e-12; - cDeterminant.Imag *= 1.0e-12; - *pExponent += 12; - Norm = NORM(cDeterminant); - } - while (Norm < 1.0e-12) { - cDeterminant.Real *= 1.0e12; - cDeterminant.Imag *= 1.0e12; - *pExponent -= 12; - Norm = NORM(cDeterminant); - } - } - } - - /* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */ - Norm = NORM(cDeterminant); - if (Norm != 0.0) { - while (Norm >= 10.0) { - cDeterminant.Real *= 0.1; - cDeterminant.Imag *= 0.1; - (*pExponent)++; - Norm = NORM(cDeterminant); - } - while (Norm < 1.0) { - cDeterminant.Real *= 10.0; - cDeterminant.Imag *= 10.0; - (*pExponent)--; - Norm = NORM(cDeterminant); - } - } - if (Matrix->NumberOfInterchangesIsOdd) - CMPLX_NEGATE(cDeterminant); - - *pDeterminant = cDeterminant.Real; - *piDeterminant.value() = cDeterminant.Imag; - } -#endif /* spCOMPLEX */ -#if REAL AND spCOMPLEX - else -#endif #if REAL { /* Real Case. */ *pDeterminant = 1.0; @@ -1074,123 +689,6 @@ void spStripFills(char* eMatrix) } #endif -#if TRANSLATE AND DELETE -/* - * DELETE A ROW AND COLUMN FROM THE MATRIX - * - * Deletes a row and a column from a matrix. - * - * Sparse will abort if an attempt is made to delete a row or column that - * doesn't exist. - * - * >>> Arguments: - * eMatrix (char *) - * Pointer to the matrix in which the row and column are to be deleted. - * Row (int) - * Row to be deleted. - * Col (int) - * Column to be deleted. - * - * >>> Local variables: - * ExtCol (int) - * The external column that is being deleted. - * ExtRow (int) - * The external row that is being deleted. - * pElement (ElementPtr) - * Pointer to an element in the matrix. Used when scanning rows and - * columns in order to eliminate elements from the last row or column. - * ppElement (ElementPtr *) - * Pointer to the location of an ElementPtr. This location will be - * filled with a NULL pointer if it is the new last element in its row - * or column. - * pElement (ElementPtr) - * Pointer to an element in the last row or column of the matrix. - * Size (int) - * The local version Matrix->Size, the size of the matrix. - */ - -void spDeleteRowAndCol(char* eMatrix, int Row, int Col) -{ - MatrixPtr Matrix = (MatrixPtr)eMatrix; - ElementPtr pElement, *ppElement, pLastElement; - int Size, ExtRow, ExtCol; - - /* Begin `spDeleteRowAndCol'. */ - - ASSERT(IS_SPARSE(Matrix) AND Row > 0 AND Col > 0); - ASSERT(Row <= Matrix->ExtSize AND Col <= Matrix->ExtSize); - - Size = Matrix->Size; - ExtRow = Row; - ExtCol = Col; - if (NOT Matrix->RowsLinked) - spcLinkRows(Matrix); - - Row = Matrix->ExtToIntRowMap[Row]; - Col = Matrix->ExtToIntColMap[Col]; - ASSERT(Row > 0 AND Col > 0); - - /* Move Row so that it is the last row in the matrix. */ - if (Row != Size) - spcRowExchange(Matrix, Row, Size); - - /* Move Col so that it is the last column in the matrix. */ - if (Col != Size) - spcColExchange(Matrix, Col, Size); - - /* Correct Diag pointers. */ - if (Row == Col) - SWAP(ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size]) - else { - Matrix->Diag[Row] = spcFindElementInCol(Matrix, Matrix->FirstInCol + Row, - Row, Row, NO); - Matrix->Diag[Col] = spcFindElementInCol(Matrix, Matrix->FirstInCol + Col, - Col, Col, NO); - } - - /* - * Delete last row and column of the matrix. - */ - /* Break the column links to every element in the last row. */ - pLastElement = Matrix->FirstInRow[Size]; - while (pLastElement != NULL) { - ppElement = &(Matrix->FirstInCol[pLastElement->Col]); - while ((pElement = *ppElement) != NULL) { - if (pElement == pLastElement) - *ppElement = NULL; /* Unlink last element in column. */ - else - ppElement = &pElement->NextInCol; /* Skip element. */ - } - pLastElement = pLastElement->NextInRow; - } - - /* Break the row links to every element in the last column. */ - pLastElement = Matrix->FirstInCol[Size]; - while (pLastElement != NULL) { - ppElement = &(Matrix->FirstInRow[pLastElement->Row]); - while ((pElement = *ppElement) != NULL) { - if (pElement == pLastElement) - *ppElement = NULL; /* Unlink last element in row. */ - else - ppElement = &pElement->NextInRow; /* Skip element. */ - } - pLastElement = pLastElement->NextInCol; - } - - /* Clean up some details. */ - Matrix->Size = Size - 1; - Matrix->Diag[Size] = NULL; - Matrix->FirstInRow[Size] = NULL; - Matrix->FirstInCol[Size] = NULL; - Matrix->CurrentSize--; - Matrix->ExtToIntRowMap[ExtRow] = -1; - Matrix->ExtToIntColMap[ExtCol] = -1; - Matrix->NeedsOrdering = YES; - - return; -} -#endif - #if PSEUDOCONDITION /* * CALCULATE PSEUDOCONDITION @@ -1321,23 +819,14 @@ RealNumber spCondition(char* eMatrix, RealNumber NormOfMatrix, int* pError) return 0.0; } -#if spCOMPLEX - if (Matrix->Complex) - return ComplexCondition(Matrix, NormOfMatrix, pError); -#endif - #if REAL Size = Matrix->Size; T = Matrix->Intermediate; -#if spCOMPLEX - Tm = Matrix->Intermediate + Size; -#else Tm = ALLOC(RealNumber, Size + 1); if (Tm == NULL) { *pError = spNO_MEMORY; return 0.0; } -#endif for (I = Size; I > 0; I--) T[I] = 0.0; @@ -1484,9 +973,7 @@ RealNumber spCondition(char* eMatrix, RealNumber NormOfMatrix, int* pError) for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]); -#if NOT spCOMPLEX FREE(Tm); -#endif Linpack = ASy / ASz; OLeary = E / MaxY; @@ -1495,210 +982,6 @@ RealNumber spCondition(char* eMatrix, RealNumber NormOfMatrix, int* pError) #endif /* REAL */ } -#if spCOMPLEX -/* - * ESTIMATE CONDITION NUMBER - * - * Complex version of spCondition(). - * - * >>> Returns: - * The reciprocal of the condition number. - * - * >>> Arguments: - * Matrix (MatrixPtr) - * Pointer to the matrix. - * NormOfMatrix (RealNumber) - * The L-infinity norm of the unfactored matrix as computed by - * spNorm(). - * pError (int *) - * Used to return error code. - * - * >>> Possible errors: - * spNO_MEMORY - */ - -static RealNumber ComplexCondition(MatrixPtr Matrix, RealNumber NormOfMatrix, int* pError) -{ - ElementPtr pElement; - ComplexVector T, Tm; - int I, K, Row; - ElementPtr pPivot; - int Size; - RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor; - RealNumber Linpack, OLeary, InvNormOfInverse; - ComplexNumber Wp, Wm; - - /* Begin `ComplexCondition'. */ - - Size = Matrix->Size; - T = (ComplexVector)Matrix->Intermediate; - Tm = ALLOC(ComplexNumber, Size + 1); - if (Tm == NULL) { - *pError = spNO_MEMORY; - return 0.0; - } - for (I = Size; I > 0; I--) - T[I].Real = T[I].Imag = 0.0; - - /* - * Part 1. Ay = e. - * Solve Ay = LUy = e where e consists of +1 and -1 terms with the sign - * chosen to maximize the size of w in Lw = e. Since the terms in w can - * get very large, scaling is used to avoid overflow. - */ - - /* Forward elimination. Solves Lw = e while choosing e. */ - E = 1.0; - for (I = 1; I <= Size; I++) { - pPivot = Matrix->Diag[I]; - if (T[I].Real < 0.0) - Em = -E; - else - Em = E; - Wm = T[I]; - Wm.Real += Em; - ASm = CMPLX_1_NORM(Wm); - CMPLX_MULT_ASSIGN(Wm, *pPivot); - if (CMPLX_1_NORM(Wm) > SLACK) { - ScaleFactor = 1.0 / MAX(SQR(SLACK), CMPLX_1_NORM(Wm)); - for (K = Size; K > 0; K--) - SCLR_MULT_ASSIGN(T[K], ScaleFactor); - E *= ScaleFactor; - Em *= ScaleFactor; - ASm *= ScaleFactor; - SCLR_MULT_ASSIGN(Wm, ScaleFactor); - } - Wp = T[I]; - Wp.Real -= Em; - ASp = CMPLX_1_NORM(Wp); - CMPLX_MULT_ASSIGN(Wp, *pPivot); - - /* Update T for both values of W, minus value is placed in Tm. */ - pElement = pPivot->NextInCol; - while (pElement != NULL) { - Row = pElement->Row; - /* Cmplx expr: Tm[Row] = T[Row] - (Wp * *pElement). */ - CMPLX_MULT_SUBT(Tm[Row], Wm, *pElement, T[Row]); - /* Cmplx expr: T[Row] -= Wp * *pElement. */ - CMPLX_MULT_SUBT_ASSIGN(T[Row], Wm, *pElement); - ASp += CMPLX_1_NORM(T[Row]); - ASm += CMPLX_1_NORM(Tm[Row]); - pElement = pElement->NextInCol; - } - - /* If minus value causes more growth, overwrite T with its values. */ - if (ASm > ASp) { - T[I] = Wm; - pElement = pPivot->NextInCol; - while (pElement != NULL) { - T[pElement->Row] = Tm[pElement->Row]; - pElement = pElement->NextInCol; - } - } else - T[I] = Wp; - } - - /* Compute 1-norm of T, which now contains w, and scale ||T|| to 1/SLACK. */ - for (ASw = 0.0, I = Size; I > 0; I--) - ASw += CMPLX_1_NORM(T[I]); - ScaleFactor = 1.0 / (SLACK * ASw); - if (ScaleFactor < 0.5) { - for (I = Size; I > 0; I--) - SCLR_MULT_ASSIGN(T[I], ScaleFactor); - E *= ScaleFactor; - } - - /* Backward Substitution. Solves Uy = w.*/ - for (I = Size; I >= 1; I--) { - pElement = Matrix->Diag[I]->NextInRow; - while (pElement != NULL) { /* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */ - CMPLX_MULT_SUBT_ASSIGN(T[I], T[pElement->Col], *pElement); - pElement = pElement->NextInRow; - } - if (CMPLX_1_NORM(T[I]) > SLACK) { - ScaleFactor = 1.0 / MAX(SQR(SLACK), CMPLX_1_NORM(T[I])); - for (K = Size; K > 0; K--) - SCLR_MULT_ASSIGN(T[K], ScaleFactor); - E *= ScaleFactor; - } - } - - /* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */ - for (ASy = 0.0, I = Size; I > 0; I--) - ASy += CMPLX_1_NORM(T[I]); - ScaleFactor = 1.0 / (SLACK * ASy); - if (ScaleFactor < 0.5) { - for (I = Size; I > 0; I--) - SCLR_MULT_ASSIGN(T[I], ScaleFactor); - ASy = 1.0 / SLACK; - E *= ScaleFactor; - } - - /* Compute infinity-norm of T for O'Leary's estimate. */ - for (MaxY = 0.0, I = Size; I > 0; I--) - if (MaxY < CMPLX_1_NORM(T[I])) - MaxY = CMPLX_1_NORM(T[I]); - - /* - * Part 2. A* z = y where the * represents the transpose. - * Recall that A = LU implies A* = U* L*. - */ - - /* Forward elimination, U* v = y. */ - for (I = 1; I <= Size; I++) { - pElement = Matrix->Diag[I]->NextInRow; - while (pElement != NULL) { /* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */ - CMPLX_MULT_SUBT_ASSIGN(T[pElement->Col], T[I], *pElement); - pElement = pElement->NextInRow; - } - if (CMPLX_1_NORM(T[I]) > SLACK) { - ScaleFactor = 1.0 / MAX(SQR(SLACK), CMPLX_1_NORM(T[I])); - for (K = Size; K > 0; K--) - SCLR_MULT_ASSIGN(T[K], ScaleFactor); - ASy *= ScaleFactor; - } - } - - /* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */ - for (ASv = 0.0, I = Size; I > 0; I--) - ASv += CMPLX_1_NORM(T[I]); - ScaleFactor = 1.0 / (SLACK * ASv); - if (ScaleFactor < 0.5) { - for (I = Size; I > 0; I--) - SCLR_MULT_ASSIGN(T[I], ScaleFactor); - ASy *= ScaleFactor; - } - - /* Backward Substitution, L* z = v. */ - for (I = Size; I >= 1; I--) { - pPivot = Matrix->Diag[I]; - pElement = pPivot->NextInCol; - while (pElement != NULL) { /* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */ - CMPLX_MULT_SUBT_ASSIGN(T[I], T[pElement->Row], *pElement); - pElement = pElement->NextInCol; - } - CMPLX_MULT_ASSIGN(T[I], *pPivot); - if (CMPLX_1_NORM(T[I]) > SLACK) { - ScaleFactor = 1.0 / MAX(SQR(SLACK), CMPLX_1_NORM(T[I])); - for (K = Size; K > 0; K--) - SCLR_MULT_ASSIGN(T[K], ScaleFactor); - ASy *= ScaleFactor; - } - } - - /* Compute 1-norm of T, which now contains z. */ - for (ASz = 0.0, I = Size; I > 0; I--) - ASz += CMPLX_1_NORM(T[I]); - - FREE(Tm); - - Linpack = ASy / ASz; - OLeary = E / MaxY; - InvNormOfInverse = MIN(Linpack, OLeary); - return InvNormOfInverse / NormOfMatrix; -} -#endif /* spCOMPLEX */ - /* * L-INFINITY MATRIX NORM * @@ -1727,33 +1010,17 @@ RealNumber spNorm(char* eMatrix) if (NOT Matrix->RowsLinked) spcLinkRows(Matrix); -/* Compute row sums. */ + /* Compute row sums. */ #if REAL - if (NOT Matrix->Complex) { - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInRow[I]; - AbsRowSum = 0.0; - while (pElement != NULL) { - AbsRowSum += ABS(pElement->Real); - pElement = pElement->NextInRow; - } - if (Max < AbsRowSum) - Max = AbsRowSum; - } - } -#endif -#if spCOMPLEX - if (Matrix->Complex) { - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInRow[I]; - AbsRowSum = 0.0; - while (pElement != NULL) { - AbsRowSum += CMPLX_1_NORM(*pElement); - pElement = pElement->NextInRow; - } - if (Max < AbsRowSum) - Max = AbsRowSum; + for (I = Matrix->Size; I > 0; I--) { + pElement = Matrix->FirstInRow[I]; + AbsRowSum = 0.0; + while (pElement != NULL) { + AbsRowSum += ABS(pElement->Real); + pElement = pElement->NextInRow; } + if (Max < AbsRowSum) + Max = AbsRowSum; } #endif return Max; @@ -1833,14 +1100,13 @@ RealNumber spLargestElement(char* eMatrix) int I; RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0; RealNumber Pivot; - ComplexNumber cPivot; ElementPtr pElement, pDiag; /* Begin `spLargestElement'. */ ASSERT(IS_SPARSE(Matrix)); #if REAL - if (Matrix->Factored AND NOT Matrix->Complex) { + if (Matrix->Factored) { if (Matrix->Error == spSINGULAR) return 0.0; @@ -1871,7 +1137,7 @@ RealNumber spLargestElement(char* eMatrix) if (AbsColSum > MaxCol) MaxCol = AbsColSum; } - } else if (NOT Matrix->Complex) { + } else { for (I = 1; I <= Matrix->Size; I++) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { @@ -1883,51 +1149,6 @@ RealNumber spLargestElement(char* eMatrix) } return Max; } -#endif -#if spCOMPLEX - if (Matrix->Factored AND Matrix->Complex) { - if (Matrix->Error == spSINGULAR) - return 0.0; - - /* Find the bound on the size of the largest element over all factorization. */ - for (I = 1; I <= Matrix->Size; I++) { - pDiag = Matrix->Diag[I]; - - /* Lower triangular matrix. */ - CMPLX_RECIPROCAL(cPivot, *pDiag); - Mag = CMPLX_1_NORM(cPivot); - if (Mag > MaxRow) - MaxRow = Mag; - pElement = Matrix->FirstInRow[I]; - while (pElement != pDiag) { - Mag = CMPLX_1_NORM(*pElement); - if (Mag > MaxRow) - MaxRow = Mag; - pElement = pElement->NextInRow; - } - - /* Upper triangular matrix. */ - pElement = Matrix->FirstInCol[I]; - AbsColSum = 1.0; /* Diagonal of U is unity. */ - while (pElement != pDiag) { - AbsColSum += CMPLX_1_NORM(*pElement); - pElement = pElement->NextInCol; - } - if (AbsColSum > MaxCol) - MaxCol = AbsColSum; - } - } else if (Matrix->Complex) { - for (I = 1; I <= Matrix->Size; I++) { - pElement = Matrix->FirstInCol[I]; - while (pElement != NULL) { - Mag = CMPLX_1_NORM(*pElement); - if (Mag > Max) - Max = Mag; - pElement = pElement->NextInCol; - } - } - return Max; - } #endif return MaxRow * MaxCol; }