/* * ====================================================================== * Functions for Matrices having Indirect storage pattern. * ====================================================================== * Copyright (C) 1993-96 by Mark Austin and David Mazzoni. * * This software is provided "as is" without express or implied warranty. * Permission is granted to use this software on any computer system, * and to redistribute it freely, subject to the following restrictions: * * 1. The authors are not responsible for the consequences of use of * this software, even if they arise from defects in the software. * 2. The origin of this software must not be misrepresented, either * by explicit claim or by omission. * 3. Altered versions must be plainly marked as such, and must not * be misrepresented as being the original software. * 4. This notice is to remain intact. * * Written by: Mark Austin July 1992 - October 1993 * ====================================================================== */ #include #include #ifdef __STDC__ #include #include #include #endif #include "defs.h" #include "matrix.h" #include "vector.h" #include "miscellaneous.h" /* #define DEBUG */ /* #define LU_DEBUG */ /* * ======================================================================= * matAllocIndirect() : Allocate memory for Matrix data structure with * eIndirect storage pattern * * Input : char *cpName -- Pointer to name of matrix. * : DATA_TYPE Type -- Data type to be stored in Matrix. * : int iNoRows -- No of Rows in Matrix. * : int iNoColumns -- No of Columns in Matrix. * Output : MATRIX *Matrix -- Pointer to matrix data structure. * ======================================================================= */ #ifdef __STDC__ MATRIX *matAllocIndirect( char *cpName, DATA_TYPE Type, int iNoRows, int iNoColumns) #else /* Start case not STDC */ MATRIX *matAllocIndirect( cpName, Type, iNoRows, iNoColumns) char *cpName; DATA_TYPE Type; int iNoRows; int iNoColumns; #endif /* End case not STDC */ { MATRIX *spA; /* [a] : Allocate matrix parent data structure, and cpName */ spA = (MATRIX *) safeMalloc( sizeof(MATRIX), __FILE__, __LINE__ ); if(cpName != NULL) spA->cpName = saveString(cpName, __FILE__, __LINE__ ); else spA->cpName = (char *) NULL; /* [b] : Set parameters and allocate memory for matrix uMatrix */ spA->Rep = Indirect; spA->Type = Type; spA->iNoRows = iNoRows; spA->iNoColumns = iNoColumns; switch((int) spA->Type) { case DoubleArray: spA->uMatrix.dpp = matAllocIndirectDouble( iNoRows, iNoColumns); break; case IntegerArray: spA->uMatrix.ipp = matAllocIndirectInteger( iNoRows, iNoColumns); break; case ComplexArray: fatalError("In matAllocIndirect() : spA->Type not implemented", (char *) NULL); break; default: fatalError("In matAllocIndirect() : Undefined spA->Type", (char *) NULL); break; } return (spA); } /* * ======================================================================= * matAllocIndirectDouble() : Allocate eIndirect storage pattern for * DOUBLE data types * * Input : int iNoRows -- No of Rows in Matrix. * : int iNoColumns -- No of Columns in Matrix. * Output : double **Matrix -- Pointer to matrix data structure. * ======================================================================= */ #ifdef __STDC__ double **matAllocIndirectDouble( int iNoRows, int iNoColumns) #else /* Start case not STDC */ double **matAllocIndirectDouble( iNoRows, iNoColumns) int iNoRows; int iNoColumns; #endif /* End case not STDC */ { double **Matrix; int ii; Matrix = (double **) safeCalloc( iNoRows, sizeof(double *), __FILE__, __LINE__ ); for(ii = 1; ii <= iNoRows; ii++) Matrix[ii-1] = (double *) safeCalloc( iNoColumns, sizeof(double), __FILE__, __LINE__ ); return (Matrix); } /* * ======================================================================= * matAllocIndirectInteger() : Allocate eIndirect storage pattern for * INTEGER data types * * Input : int iNoRows -- No of Rows in Matrix. * : int iNoColumns -- No of Columns in Matrix. * Output : double **Matrix -- Pointer to matrix data structure. * ======================================================================= */ #ifdef __STDC__ int **matAllocIndirectInteger( int iNoRows, int iNoColumns) #else /* Start case not STDC */ int **matAllocIndirectInteger( iNoRows, iNoColumns) int iNoRows; int iNoColumns; #endif /* End case not STDC */ { int **Matrix; int ii; Matrix = (int **) safeCalloc( iNoRows, sizeof(int *), __FILE__, __LINE__ ); for(ii = 1; ii <= iNoRows; ii++) Matrix[ii-1] = (int *) safeCalloc( iNoColumns, sizeof(int), __FILE__, __LINE__ ); return (Matrix); } /* * ======================================================================= * matPrintIndirectDouble() : Print a Matrix [iNoRows][iNoColumns] of * data type DOUBLE * * matPrintIndirectInteger() : Print a Matrix [iNoRows][iNoColumns] of * data type INTEGER * * Where -- ColumnsAcrossPage = Number of columns printed across page. * spA->iNoRows = Number of rows in matrix. * spA->iNoColumns = Number of columns in matrix. * iFirstColumn = Number of first column in block * iLastColumn = Number of last column in block * ib = Current No of Matrix Block. * * Input : MATRIX *spA -- Pointer to matrix data structure. * Output : void * ======================================================================= */ enum { ColumnsAcrossPage = 6 }; /* Item [a] */ #ifdef __STDC__ void matPrintIndirectDouble( MATRIX *spA ) #else /* Start case not STDC */ void matPrintIndirectDouble( spA ) MATRIX *spA; #endif /* End case not STDC */ { int ii, ij, ib; int iNoBlocks; int iFirstColumn, iLastColumn; /* [a] : Compute no of blocks of rows to be printed */ /* Item [c] */ if( spA->iNoColumns % ((int) ColumnsAcrossPage) == 0 ) iNoBlocks = (spA->iNoColumns/((int) ColumnsAcrossPage)); else iNoBlocks = (spA->iNoColumns/((int) ColumnsAcrossPage)) + 1; /* [b] : Loop over blocks of rows */ for( ib = 1; ib <= iNoBlocks; ib++ ) { /* Item [d] */ iFirstColumn = (ib-1)*((int) ColumnsAcrossPage) + 1; iLastColumn = MIN( ib*((int) ColumnsAcrossPage) , spA->iNoColumns ); /* [c] : Print title of matrix at top of each block */ /* Item [e] */ if( spA->cpName != NULL ) printf("\nMATRIX : \"%s\"\n\n", spA->cpName); else printf("\nMATRIX : \"UNTITLED\"\n\n"); /* [d] : Label row and column nos */ printf ("row/col "); for( ii = iFirstColumn; ii <= iLastColumn; ii++ ) printf("%3d ", ii); printf("\n"); /* [e] : Print Contents of Matrix */ /* Item [f] */ for( ii = 1; ii <= spA->iNoRows; ii++ ) { printf(" %3d ", ii); for( ij = iFirstColumn; ij <= iLastColumn; ij++) printf(" %12.5e", spA->uMatrix.dpp[ ii-1 ][ ij-1 ]); printf("\n"); } } } #ifdef __STDC__ void matPrintIndirectInteger( MATRIX *spA ) #else /* Start case not STDC */ void matPrintIndirectInteger( spA ) MATRIX *spA; #endif /* End case not STDC */ { int ii, ij, ib; int iNoBlocks; int iFirstColumn, iLastColumn; /* [a] : Compute no of blocks of rows to be printed */ /* Item [c] */ if( spA->iNoColumns % ((int) ColumnsAcrossPage) == 0 ) iNoBlocks = (spA->iNoColumns/((int) ColumnsAcrossPage)); else iNoBlocks = (spA->iNoColumns/((int) ColumnsAcrossPage)) + 1; /* [b] : Loop over blocks of rows */ for( ib = 1; ib <= iNoBlocks; ib++ ) { /* Item [d] */ iFirstColumn = (ib-1)*((int) ColumnsAcrossPage) + 1; iLastColumn = MIN( ib*((int) ColumnsAcrossPage) , spA->iNoColumns ); /* [c] : Print title of matrix at top of each block */ /* Item [e] */ if( spA->cpName != NULL ) printf("\nMATRIX : \"%s\"\n\n", spA->cpName); else printf("\nMATRIX : \"UNTITLED\"\n\n"); /* [d] : Label row and column nos */ printf ("row/col "); for( ii = iFirstColumn; ii <= iLastColumn; ii++ ) printf("%3d ", ii); printf("\n"); /* [e] : Print Contents of Matrix */ /* Item [f] */ for( ii = 1; ii <= spA->iNoRows; ii++ ) { printf(" %3d ", ii); for( ij = iFirstColumn; ij <= iLastColumn; ij++) printf(" %12d", spA->uMatrix.ipp[ ii-1 ][ ij-1 ]); printf("\n"); } } } /* * ======================================================================= * matFreeIndirect() : Free memory for eIndirect storage. * * Input : MATRIX *Matrix -- Pointer to matrix data structure. * Output : void * ======================================================================= */ #ifdef __STDC__ void matFreeIndirect( MATRIX *spA ) #else /* Start case not STDC */ void matFreeIndirect( spA ) MATRIX *spA; #endif /* End case not STDC */ { /* [a] : Free body of matrix[iNoRows][iNoColumns] */ switch(spA->Type) { case DoubleArray: matFreeIndirectDouble( spA->uMatrix.dpp, spA->iNoRows); break; case IntegerArray: matFreeIndirectInteger( spA->uMatrix.ipp, spA->iNoRows); break; case ComplexArray: fatalError("In matFreeIndirect() : spA->Type not implemented", (char *) NULL); break; default: break; } /* (b) Free cpName and parent data structure for matrix */ free ((char *) spA->cpName); free ((char *) spA); } /* * ======================================================================= * matFreeIndirectDouble() : Free memory for eIndirect storage of data * type DOUBLE. * * Input : double **d -- Pointer to matrix of data type double. * : int iNoRows -- No of Rows in Matrix. * Output : void * ======================================================================= */ #ifdef __STDC__ void matFreeIndirectDouble( double **d, int iNoRows) #else /* Start case not STDC */ void matFreeIndirectDouble( d, iNoRows) double **d; int iNoRows; #endif /* End case not STDC */ { int ii; for(ii = 1; ii <= iNoRows; ii++) free ((char *) d[ii-1]); free ((char *) d); } /* * ======================================================================= * matFreeIndirectInteger() : Free memory for eIndirect storage of data * type INTEGER. * * Input : int **i -- Pointer to matrix of data type double. * : int iNoRows -- No of Rows in Matrix. * Output : void * ======================================================================= */ #ifdef __STDC__ void matFreeIndirectInteger( int **i, int iNoRows) #else /* Start case not STDC */ void matFreeIndirectInteger( i, iNoRows) int **i; int iNoRows; #endif /* End case not STDC */ { int ii; for(ii = 1; ii <= iNoRows; ii++) free ((char *) i[ii-1]); free ((char *) i); } /* * ================================================================= * matScaleIndirectDouble() : Multiply matrix by scalar [c] = d.[A]. * * Input : MATRIX *spA -- Pointer to (mxn) matrix A. * : double dScaleFactor -- scale factor. * Output : MATRIX *spC -- Pointer to scaled matrix C. * ================================================================= */ #ifdef __STDC__ MATRIX *matScaleIndirectDouble( MATRIX *spA, double dScale) #else /* Start case not STDC */ MATRIX *matScaleIndirectDouble( spA, dScale) MATRIX *spA; double dScale; #endif /* End case not STDC */ { MATRIX *spB; int ii, ij; /* [a] : Check Input matrix [A] */ if(spA == (MATRIX *) NULL) fatalError("In matScaleIndirectDouble() : [A] == NULL", (char *) NULL); /* [b] : Scale Indirect Double Matrix */ spB = matCopyIndirectDouble( spA ); for( ii = 1; ii <= spA->iNoRows; ii++) for( ij = 1; ij <= spA->iNoColumns; ij++) spB->uMatrix.dpp[ii-1][ij-1] = dScale * spA->uMatrix.dpp[ii-1][ij-1]; return ( spB); } /* * ======================================================================= * matAddIndirectDouble() : Add Matrices [C] = [A] + [B]. * * Input : MATRIX *spA -- Pointer to (nxn) matrix A. * : MATRIX *spB -- Pointer to (nxn) matrix B. * Output : MATRIX *spC -- Pointer to (nxn) matrix C. * ======================================================================= */ #ifdef __STDC__ MATRIX *matAddIndirectDouble( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matAddIndirectDouble( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { MATRIX *spC; int ii, ij; /* [a] : Check compatibility of matrix dimensions */ if(spA->iNoRows != spB->iNoRows || spA->iNoColumns != spB->iNoColumns) { printf(" ERROR >> spA->iNoRows = %4d spA->iNoColumns = %4d\n", spA->iNoRows, spA->iNoColumns ); printf(" ERROR >> spB->iNoRows = %4d spB->iNoColumns = %4d\n", spB->iNoRows, spB->iNoColumns ); fatalError("Execution halted in matAddIndirectDouble()", "Inconsistent Matrix Dimensions", (char *) NULL); } /* [b] : Add Matrices : [C] = [A] + [B] */ spC = matAllocIndirect((char *) NULL, DoubleArray, spA->iNoRows, spB->iNoColumns); for(ij = 1; ij <= spA->iNoColumns; ij++) for(ii = 1; ii <= spA->iNoRows; ii++) spC->uMatrix.dpp[ii-1][ij-1] = spA->uMatrix.dpp[ ii-1 ][ ij-1 ] + spB->uMatrix.dpp[ ii-1 ][ ij-1 ]; return ( spC ); } /* * ======================================================================= * matAddReplaceIndirectDouble() : Matrix Replacement [A] = [A] + [B]. * * Input : MATRIX *spA -- Pointer to (nxn) matrix A. * : MATRIX *spB -- Pointer to (nxn) matrix B. * Output : MATRIX *spC -- Pointer to (nxn) matrix A. * ======================================================================= */ #ifdef __STDC__ MATRIX *matAddReplaceIndirectDouble( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matAddReplaceIndirectDouble( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { int ii, ij; /* [a] : Compare Dimensions of Matrices */ if(spA->iNoRows != spB->iNoRows || spA->iNoColumns != spB->iNoColumns) { printf(" ERROR >> spA->iNoRows = %4d spA->iNoColumns = %4d\n", spA->iNoRows, spA->iNoColumns ); printf(" ERROR >> spB->iNoRows = %4d spB->iNoColumns = %4d\n", spB->iNoRows, spB->iNoColumns ); fatalError("Execution halted in matAddReplaceIndirectDouble()", "Inconsistent Matrix Dimensions", (char *) NULL); } /* [b] : Matrix Replacement : [A] = [A] + [B] */ for(ij = 1; ij <= spA->iNoColumns; ij++) for(ii = 1; ii <= spA->iNoRows; ii++) spA->uMatrix.dpp[ii-1][ij-1] += spB->uMatrix.dpp[ ii-1 ][ ij-1 ]; return (spA); } /* * ======================================================================= * matSubIndirectDouble() : Add Matrices [C] = [A] - [B]. * * Input : MATRIX *spA -- Pointer to (nxn) matrix A. * : MATRIX *spB -- Pointer to (nxn) matrix B. * Output : MATRIX *spC -- Pointer to (nxn) matrix C. * ======================================================================= */ #ifdef __STDC__ MATRIX *matSubIndirectDouble( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matSubIndirectDouble( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { MATRIX *spC; int ii, ij; /* [a] : Compare Dimensions of Matrices */ if(spA->iNoRows != spB->iNoRows || spA->iNoColumns != spB->iNoColumns) { printf(" ERROR >> spA->iNoRows = %4d spA->iNoColumns = %4d\n", spA->iNoRows, spA->iNoColumns ); printf(" ERROR >> spB->iNoRows = %4d spB->iNoColumns = %4d\n", spB->iNoRows, spB->iNoColumns ); fatalError("Execution halted in matSubnndirectDouble()", "Inconsistent Matrix Dimensions", (char *) NULL); } /* [b] : Subtract Matrices : [C] = [A] - [B] */ spC = matAllocIndirect((char *) NULL, DoubleArray, spA->iNoRows, spB->iNoColumns); for(ij = 1; ij <= spA->iNoColumns; ij++) for(ii = 1; ii <= spA->iNoRows; ii++) spC->uMatrix.dpp[ ii-1 ][ ij-1 ] = spA->uMatrix.dpp[ ii-1 ][ ij-1 ] - spB->uMatrix.dpp[ ii-1 ][ ij-1 ]; return (spC); } /* * ======================================================================= * matSubReplaceIndirectDouble() : Matrix Replacement [A] = [A] - [B]. * * Input : MATRIX *spA -- Pointer to (nxn) matrix A. * : MATRIX *spB -- Pointer to (nxn) matrix B. * Output : MATRIX *spC -- Pointer to (nxn) matrix A. * ======================================================================= */ #ifdef __STDC__ MATRIX *matSubReplaceIndirectDouble( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matSubReplaceIndirectDouble( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { int ii, ij; /* [a] : Compare Dimensions of Matrices */ if(spA->iNoRows != spB->iNoRows || spA->iNoColumns != spB->iNoColumns) { printf(" ERROR >> spA->iNoRows = %4d spA->iNoColumns = %4d\n", spA->iNoRows, spA->iNoColumns ); printf(" ERROR >> spB->iNoRows = %4d spB->iNoColumns = %4d\n", spB->iNoRows, spB->iNoColumns ); fatalError("Execution halted in matSubReplaceIndirectDouble()", "Inconsistent Matrix Dimensions", (char *) NULL); } /* [b] : Replace Matrices : [A] = [A] - [B] */ for(ij = 1; ij <= spA->iNoColumns; ij++) for(ii = 1; ii <= spA->iNoRows; ii++) spA->uMatrix.dpp[ ii-1 ][ ij-1 ] -= spB->uMatrix.dpp[ ii-1 ][ ij-1 ]; return (spA); } /* * ======================================================================= * matMultIndirectDouble() : Matrix Multiplication [C] = [A].[B]. * * Input : MATRIX *spA -- Pointer to matrix A. * : MATRIX *spB -- Pointer to matrix B. * Output : MATRIX *spC -- Pointer to matrix C. * ======================================================================= */ #ifdef __STDC__ MATRIX *matMultIndirectDouble( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matMultIndirectDouble( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { MATRIX *spC; int ii, ij, ik; /* [a] : Compare Inside Dimensions of Matrices */ if(spA->iNoColumns != spB->iNoRows ) { printf(" ERROR >> spA->iNoRows = %4d spA->iNoColumns = %4d\n", spA->iNoRows, spA->iNoColumns ); printf(" ERROR >> spB->iNoRows = %4d spB->iNoColumns = %4d\n", spB->iNoRows, spB->iNoColumns ); fatalError("Execution halted in matMultIndirectDouble()", "Inconsistent Inside Matrix Dimensions", (char *) NULL); } /* [b] : Multiply Matrices */ spC = matAllocIndirect((char *) NULL, DoubleArray, spA->iNoRows, spB->iNoColumns); for(ij = 1; ij <= spB->iNoColumns; ij++) for(ii = 1; ii <= spA->iNoRows; ii++) for(ik = 1; ik <= spA->iNoColumns; ik++) spC->uMatrix.dpp[ ii-1 ][ ij-1 ] += spA->uMatrix.dpp[ ii-1 ][ ik-1 ]*spB->uMatrix.dpp[ ik-1 ][ ij-1 ]; return ( spC ); } /* * ======================================================================= * matMultReplaceIndirectDouble() : Matrix Replacement [A] = [A] * [B]. * * Input : MATRIX *spA -- Pointer to (nxn) matrix A. * : MATRIX *spB -- Pointer to (nxn) matrix B. * Output : MATRIX *spC -- Pointer to (nxn) matrix A. * * Note : Matrices [A] and [B] must be square. * ======================================================================= */ #ifdef __STDC__ MATRIX *matMultReplaceIndirectDouble( MATRIX *spA, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matMultReplaceIndirectDouble( spA, spB ) MATRIX *spA; MATRIX *spB; #endif /* End case not STDC */ { MATRIX *spC; int ii, ij; /* [a] : Check that [A] and [B] are square */ if( spA->iNoRows != spA->iNoColumns || spB->iNoRows != spB->iNoColumns || spA->iNoColumns != spB->iNoRows ) { printf(" ERROR >> spA->iNoRows = %4d spA->iNoColumns = %4d\n", spA->iNoRows, spA->iNoColumns ); printf(" ERROR >> spB->iNoRows = %4d spB->iNoColumns = %4d\n", spB->iNoRows, spB->iNoColumns ); fatalError("Execution halted in matMultReplaceIndirectDouble()", "Inconsistent Matrix Dimensions", (char *) NULL); } /* [b] : Compute [A] = [A] * [B] */ spC = matMultIndirectDouble( spA, spB ); for(ii = 1; ii <= spA->iNoRows; ii++) for(ij = 1; ij <= spA->iNoColumns; ij++) spA->uMatrix.dpp[ ii-1 ][ ij-1 ] = spC->uMatrix.dpp[ ii-1 ][ ij-1 ]; matFreeIndirect( spC ); return (spA); } /* * ======================================================================= * matNegateIndirectDouble() : Matrix Negation [B] = -[A]. * * Input : MATRIX *spA -- Pointer to matrix A. * Output : MATRIX *spB -- Pointer to matrix B. * ======================================================================= */ #ifdef __STDC__ MATRIX *matNegateIndirectDouble( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matNegateIndirectDouble( spA ) MATRIX *spA; #endif /* End case not STDC */ { MATRIX *spB; int ii, ij; spB = matAllocIndirect((char *) NULL, DoubleArray, spA->iNoRows, spA->iNoColumns); for(ii = 1; ii <= spA->iNoRows; ii = ii + 1) for(ij = 1; ij <= spA->iNoColumns; ij = ij + 1) spB->uMatrix.dpp[ii-1][ij-1] = -spA->uMatrix.dpp[ii-1][ij-1]; return (spB); } /* * ======================================================================= * matNegateReplaceIndirectDouble() : Matrix Negation [A] = -[A]. * * Input : MATRIX *spA -- Pointer to matrix A. * Output : MATRIX *spB -- Pointer to matrix A. * ======================================================================= */ #ifdef __STDC__ MATRIX *matNegateReplaceIndirectDouble( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matNegateReplaceIndirectDouble( spA ) MATRIX *spA; #endif /* End case not STDC */ { int ii, ij; for(ii = 1; ii <= spA->iNoRows; ii = ii + 1) for(ij = 1; ij <= spA->iNoColumns; ij = ij + 1) spA->uMatrix.dpp[ii-1][ij-1] *= -1.0; return (spA); } /* * ======================================================================= * matTransposeIndirectDouble() : Matrix Transpose [B] = [A]^T. * * Input : MATRIX *spA -- Pointer to matrix A. * Output : MATRIX *spB -- Pointer to matrix B. * ======================================================================= */ #ifdef __STDC__ MATRIX *matTransposeIndirectDouble( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matTransposeIndirectDouble( spA ) MATRIX *spA; #endif /* End case not STDC */ { MATRIX *spB; int ii, ij; /* [a] : Transpose Matrix */ spB = matAllocIndirect((char *) NULL, DoubleArray, spA->iNoColumns, spA->iNoRows); for(ii = 1; ii <= spA->iNoRows; ii++) for(ij = 1; ij <= spA->iNoColumns; ij++) spB->uMatrix.dpp[ ij-1 ][ ii-1 ] = spA->uMatrix.dpp[ ii-1 ][ ij-1 ]; return (spB); } /* * =================================================== * matCopyIndirectDouble() : Matrix Copy [B] = [A]. * * Input : MATRIX *spA -- Pointer to matrix A. * Output : MATRIX *spB -- Pointer to matrix B. * =================================================== */ #ifdef __STDC__ MATRIX *matCopyIndirectDouble( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *matCopyIndirectDouble( spA ) MATRIX *spA; #endif /* End case not STDC */ { MATRIX *spB; int ii, ij; /* [a] : Make Copy of Matrix */ spB = matAllocIndirect( spA->cpName, DoubleArray, spA->iNoRows, spA->iNoColumns); for(ii = 1; ii <= spA->iNoRows; ii++) for(ij = 1; ij <= spA->iNoColumns; ij++) spB->uMatrix.dpp[ii-1][ij-1] = spA->uMatrix.dpp[ii-1][ij-1]; return ( spB ); } /* * ======================================================================== * matLUDecompositionDouble() : Use method of Crout Reduction with pivoting * and scaled equtions to decompose a (nxn) matrix [A] into [L][U]. * * Input : MATRIX *spA -- Pointer to (nxn) matrix A. * : VECTOR *spPivot -- Pointer to (nx1) pivot vector. * Output : MATRIX *spA -- Pointer to (nxn) [L][U] product. * ======================================================================== */ #ifdef __STDC__ MATRIX *matLUDecompositionIndirect( MATRIX *spA, VECTOR *spPivot) #else /* Start case not STDC */ MATRIX *matLUDecompositionIndirect( spA, spPivot ) MATRIX *spA; VECTOR *spPivot; #endif /* End case not STDC */ { VECTOR *spScale; double dSum; double dNumerator; double dDenominator; int iOrder1, iOrder2; int ii, ij, ik; /* [a] : Make sure that Matrix spA is square */ if((spA->iNoRows != spA->iNoColumns)) { printf("Problem : spA->iNoRows = %4d spA->iNoColumns = %4d\n", spA->iNoRows, spA->iNoColumns); fatalError("Execution halted in matLUDecomposition()", "Matrix spA must be square", (char *) NULL); } /* [b] : Setup Matrix of Scale Factors */ spScale = SetupScaleFactors( spA ); /* [c] : Loop over Columns and Compute Crout Reduction */ ij = 1; Pivot( spA, spScale, spPivot, ij ); iOrder1 = spPivot->uVector.ip[ 0 ]; for (ij = 2; ij <= spA->iNoColumns; ij = ij + 1) { spA->uMatrix.dpp[iOrder1-1][ij - 1] /= spA->uMatrix.dpp[iOrder1-1][ 0 ]; } for (ij = 2; ij <= spA->iNoColumns - 1; ij = ij + 1) { for (ii = ij; ii <= spA->iNoColumns; ii = ii + 1) { iOrder1 = spPivot->uVector.ip[ ii - 1 ]; dSum = 0.0; for (ik = 1; ik <= ij - 1; ik = ik + 1) { iOrder2 = spPivot->uVector.ip[ ik - 1 ]; dSum += spA->uMatrix.dpp[ iOrder1-1 ][ ik - 1 ] * spA->uMatrix.dpp[ iOrder2-1 ][ ij - 1 ]; } spA->uMatrix.dpp[ iOrder1-1 ][ ij - 1 ] -= dSum; } Pivot( spA, spScale, spPivot, ij ); iOrder1 = spPivot->uVector.ip[ ij - 1 ]; for (ik = ij + 1; ik <= spA->iNoColumns; ik = ik + 1) { dSum = 0.0; for (ii = 1; ii <= ij - 1; ii = ii + 1) { iOrder2 = spPivot->uVector.ip[ ii - 1 ]; dSum += spA->uMatrix.dpp[ iOrder1-1 ][ ii - 1 ] * spA->uMatrix.dpp[ iOrder2-1 ][ ik - 1 ]; } dNumerator = spA->uMatrix.dpp[ iOrder1-1 ][ ik-1 ] - dSum; dDenominator = spA->uMatrix.dpp[ iOrder1-1 ][ ij-1 ]; spA->uMatrix.dpp[ iOrder1-1 ][ ik - 1 ] = dNumerator/dDenominator; } } iOrder1 = spPivot->uVector.ip[ spA->iNoColumns - 1 ]; dSum = 0.0; for (ik = 1; ik <= spA->iNoColumns - 1; ik = ik + 1) { iOrder2 = spPivot->uVector.ip[ ik - 1 ]; dSum += spA->uMatrix.dpp[ iOrder1-1 ][ ik - 1 ] * spA->uMatrix.dpp[ iOrder2-1 ][ spA->iNoColumns - 1 ]; } spA->uMatrix.dpp[ iOrder1 - 1 ][ spA->iNoColumns - 1 ] -= dSum; return( spA ); } /* * ====================================================================== * matLUSubstitutionDouble() : Compute Forward and Backward Substitution. * * Input : VECTOR *spPivot -- Pointer to (nx1) pivot vector. * : MATRIX *spLU -- Pointer to (nxn) matrix [L][U]. * : MATRIX *spB -- Pointer to (nxp) matrix [B]. * Output : MATRIX spX -- Pointer to (nxp) solution matrix. * ====================================================================== */ #ifdef __STDC__ MATRIX *matLUSubstitutionIndirect( VECTOR *spPivot, MATRIX *spLU, MATRIX *spB ) #else /* Start case not STDC */ MATRIX *matLUSubstitutionIndirect( spPivot, spLU, spB) VECTOR *spPivot; MATRIX *spLU; MATRIX *spB; #endif /* End case not STDC */ { MATRIX *spDisp; double dDenominator; double dNumerator; double dSum; int iOrder1; int ii, ij, ik; /* [a] : Allocate Memory for Solution Vector */ spDisp = matAllocIndirect( (char *) NULL, DoubleArray, spLU->iNoRows, spB->iNoColumns ); /* [b] : Forward Substitution with Unscrambling of Permutations */ iOrder1 = spPivot->uVector.ip[ 0 ]; for( ik = 1; ik <= spB->iNoColumns; ik = ik + 1) { spDisp->uMatrix.dpp[0][ik-1] = spB->uMatrix.dpp[ iOrder1 - 1 ][ ik-1 ] / spLU->uMatrix.dpp[ iOrder1 - 1 ][ 0 ]; } for( ik = 1; ik <= spB->iNoColumns; ik = ik + 1) { for( ii = 2; ii <= spLU->iNoRows; ii = ii + 1) { iOrder1 = spPivot->uVector.ip[ ii - 1 ]; dSum = 0.0; for( ij = 1; ij <= ii - 1; ij = ij + 1) { dSum += spLU->uMatrix.dpp[ iOrder1 - 1 ][ ij - 1] * spDisp->uMatrix.dpp[ ij - 1][ik - 1]; } dNumerator = spB->uMatrix.dpp[ iOrder1-1 ][ ik-1 ] - dSum; dDenominator = spLU->uMatrix.dpp[ iOrder1-1 ][ ii-1 ]; spDisp->uMatrix.dpp[ ii - 1 ][ik-1] = dNumerator/dDenominator; } } /* [c] : Backsubstitution with Unscrambling of Permutations */ for( ik = 1; ik <= spB->iNoColumns; ik = ik + 1) { for( ii = spLU->iNoRows - 1; ii >= 1; ii = ii - 1) { iOrder1 = spPivot->uVector.ip[ ii - 1 ]; dSum = 0.0; for( ij = ii + 1; ij <= spLU->iNoRows; ij = ij + 1) { dSum += spLU->uMatrix.dpp[ iOrder1 - 1 ][ ij - 1] * spDisp->uMatrix.dpp[ ij - 1][ik-1]; } spDisp->uMatrix.dpp[ii-1][ik-1] -= dSum; } } return(spDisp); } /* * ================================================================ * CholeskyDecompositionIndirect() : [L][L]^T Cholesky Factorization * of Symmetric Indirect Matrix. * * Input : MATRIX *spA -- Pointer to symmetric matrix [A]. * Output : MATRIX *spA -- Pointer to [L][L]^T. * ================================================================ */ #ifdef __STDC__ MATRIX *CholeskyDecompositionIndirect( MATRIX *spA ) #else /* Start case not STDC */ MATRIX *CholeskyDecompositionIndirect( spA ) MATRIX *spA; #endif /* End case not STDC */ { double dSum; int ii, ij, ik; /* [a] : Check Input matrix [A] */ if( spA == NULL ) fatalError("In CholeskyDecompositionIndirect() : Pointer to matrix is NULL", (char *) NULL); if( spA->Type != DoubleArray ) fatalError("In CholeskyDecompositionIndirect() : spA->Type != DoubleArray", (char *) NULL); if( spA->Rep != Indirect ) fatalError("In CholeskyDecompositionIndirect() : spA->Rep != Indirect", (char *) NULL); /* [c] : Compute [A] = [L][L]^T */ for (ij = 1; ij <= spA->iNoColumns; ij = ij + 1) { for (ii = ij; ii <= spA->iNoRows; ii = ii + 1) { if( ii == ij) { dSum = 0.0; for (ik = 1; ik < ii; ik = ik + 1) dSum += spA->uMatrix.dpp[ ii-1 ][ ik-1 ] * spA->uMatrix.dpp[ ii-1 ][ ik-1 ]; if(spA->uMatrix.dpp[ ii-1 ][ ii-1 ] > dSum) spA->uMatrix.dpp[ii-1][ii-1] = sqrt(spA->uMatrix.dpp[ ii-1 ][ ii-1 ] - dSum); else fatalError("In CholeskyDecompositionIndirect()", "Matrix is not positive definite", (char *) NULL); } if( ii > ij) { dSum = 0.0; for (ik = 1; ik < ij; ik = ik + 1) dSum += spA->uMatrix.dpp[ ii-1 ][ ik-1 ] * spA->uMatrix.dpp[ ij-1 ][ ik-1 ]; spA->uMatrix.dpp[ii-1][ij-1] = (spA->uMatrix.dpp[ ii-1 ][ ij-1 ] - dSum) / spA->uMatrix.dpp[ij-1][ij-1]; spA->uMatrix.dpp[ij-1][ii-1] = spA->uMatrix.dpp[ii-1][ij-1]; } } } return ( spA ); } /* * ================================================================= * SetupScaleFactors() : Initialize vector for scale factors. If [A] * is a (nxn) matrix, then spScale will be a (nx1) vector containing * the maxium absolute value in each row of [A]. * * Input : MATRIX spA -- Pointer to (nxn) matrix [A]. * Output : VECTOR spScale -- Pointer to (nx1) vector spScale. * ================================================================= */ #ifdef __STDC__ VECTOR *SetupScaleFactors( MATRIX *spA ) #else /* Start case not STDC */ VECTOR *SetupScaleFactors( spA ) MATRIX *spA; #endif /* End case not STDC */ { VECTOR *spScale; double dTemp; int ii, ij; spScale = vecAlloc("Scale Factor", DoubleArray, spA->iNoRows ); for (ii = 1; ii <= spA->iNoRows; ii = ii + 1) { spScale->uVector.dp[ii-1] = ABS(spA->uMatrix.dpp[ii-1][0]); for (ij = 2; ij <= spA->iNoColumns; ij = ij + 1) { dTemp = ABS(spA->uMatrix.dpp[ii-1][ij-1]); spScale->uVector.dp[ii-1] = MAX( spScale->uVector.dp[ii-1], dTemp ); } if( spScale->uVector.dp[ii-1] == 0) fatalError("Execution halted in SetupScaleFactors()", "Matrix A is Singular !!", (char *) NULL); } return ( spScale ); } /* * ================================================================== * SetupPivotVector() : Initialize (nx1) vector for row permutations. * * Input : MATRIX spA -- Pointer to (nxn) matrix [A]. * Output : VECTOR spPivot -- Pointer to (nx1) vector spPivot. * ================================================================== */ #ifdef __STDC__ VECTOR *SetupPivotVector( MATRIX *spA ) #else /* Start case not STDC */ VECTOR *SetupPivotVector( spA ) MATRIX *spA; #endif /* End case not STDC */ { VECTOR *spPivot; int ii; spPivot = vecAlloc( "Row Pivots", IntegerArray, spA->iNoRows ); for (ii = 1; ii <= spA->iNoRows; ii = ii + 1) spPivot->uVector.ip[ii-1] = ii; return ( spPivot ); } /* * ================================================================== * Pivot() : Find largest element below diagonal and pivot rows * * Input : MATRIX spA -- Pointer to (nxn) matrix [A]. * : VECTOR spScale -- Pointer to (nx1) vector spScale. * : VECTOR spPivot -- Pointer to (nx1) vector spPivot. * : int iColumnNo -- Column No for Pivoting * Output : MATRIX spA -- [A] and [spPivot] are rearranged by * the pivoting procedure. * ================================================================== */ #ifdef __STDC__ Pivot( MATRIX *spA, VECTOR *spScale, VECTOR *spPivot, int iColumnNo ) #else /* Start case not STDC */ Pivot( spA, spScale, spPivot, iColumnNo ) VECTOR *spScale, *spPivot; MATRIX *spA; int iColumnNo; #endif /* End case not STDC */ { double dLargest, dTemp; int ii , iTemp, iOrder; int iPivot; /* [a] : Find Row No of Largest Scaled Equation */ iPivot = iColumnNo; iOrder = spPivot->uVector.ip[ iColumnNo - 1 ]; dLargest = spA->uMatrix.dpp[ iOrder-1 ][ iColumnNo - 1 ] / spScale->uVector.dp[ iOrder-1 ]; dLargest = ABS( dLargest ); for (ii = iColumnNo + 1; ii <= spA->iNoColumns; ii = ii + 1) { iOrder = spPivot->uVector.ip[ ii - 1 ]; dTemp = ABS((spA->uMatrix.dpp[ iOrder-1 ][ iColumnNo - 1 ] / spScale->uVector.dp[ iOrder-1 ])); if( dTemp > dLargest ) { dLargest = dTemp; iPivot = ii; } } if( dLargest == 0 ) fatalError("Execution halted in Pivot()", "Matrix A is Singular !!", (char *) NULL); /* [b] : Swap Row Nos in Order Vector */ iTemp = spPivot->uVector.ip[ iPivot - 1 ]; spPivot->uVector.ip[ iPivot - 1 ] = spPivot->uVector.ip[ iColumnNo - 1 ]; spPivot->uVector.ip[ iColumnNo - 1 ] = iTemp; } /* * =================================================================== * matGeneralisedEigen() : Solve Eigen Problem [A][x] = [B][x][Lambda]. * where [A] and [B] are small fully populated * symmetric matrices. * * We convert [A][x] = [B][x][Lambda] into [A*][y] = [y][Lambda] in * a four-step procedure: * * (a) : Compute [B] = [L].[L]^T * (b) : Solve [L][C] = [A]. * (c) : Compute [C]^T. * (d) : Solve [L][A*]^T = [C]^T; matrix [A*] is symmetric. * * Use Householder method to solve [A*][y] = [y][Lamda], followed by * backward substitutution for [L]^T.[x] = [y]. * * Written By : M. Austin November 1993 * ================================================================= */ matGeneralisedEigen( spA, spB, spEigenvalue, spEigenvector ) MATRIX *spA, *spB; MATRIX *spEigenvalue; MATRIX *spEigenvector; { MATRIX *spEigenvectorY; MATRIX *spCwork2; MATRIX *spAwork; MATRIX *spCwork; MATRIX *spLLT; int ii, ij, ik, iSize; double sum; /* [a] : Check Input and Allocate Working Arrays */ iSize = spA->iNoRows; spCwork = matAllocIndirect("Working [C]", DoubleArray, iSize, iSize ); spAwork = matAllocIndirect("Working [A]", DoubleArray, iSize, iSize ); /* [b] : Compute Cholesky Decomposition [B] = [L].[L]^T */ spLLT = matCopy( spB ); spLLT = CholeskyDecompositionIndirect( spLLT ); /* [c] : Solve [L][C] = [A] via Forward Substitution */ for (ij = 1; ij <= spA->iNoRows; ij = ij + 1) { for (ii = 1; ii <= spA->iNoColumns; ii = ii + 1) { sum = 0.0; for (ik = 1; ik < ii; ik = ik + 1) sum += spLLT->uMatrix.dpp[ ii-1 ][ ik-1 ] * spCwork->uMatrix.dpp[ ik-1 ][ij-1]; spCwork->uMatrix.dpp[ii-1][ij-1] = (spA->uMatrix.dpp[ ii-1 ][ ij-1 ] - sum) / spLLT->uMatrix.dpp[ ii-1 ][ ii-1 ]; } } /* [d] : Replace [C] by [C]^T. */ spCwork2 = matTranspose ( spCwork ); /* [e] : Solve [L][A*] = [C]^T via Forward Substitution */ for (ij = 1; ij <= spA->iNoRows; ij = ij + 1) { for (ii = 1; ii <= spA->iNoColumns; ii = ii + 1) { sum = 0.0; for (ik = 1; ik < ii; ik = ik + 1) sum += spLLT->uMatrix.dpp[ ii-1 ][ ik-1 ] * spAwork->uMatrix.dpp[ ik-1 ][ij-1]; spAwork->uMatrix.dpp[ii-1][ij-1] = (spCwork2->uMatrix.dpp[ ii-1 ][ ij-1 ] - sum) / spLLT->uMatrix.dpp[ ii-1 ][ ii-1 ]; } } /* [f] : Use Householder Method to Solve Standard Eigenvalue Problem */ iSize = spAwork->iNoRows; spEigenvectorY = matAllocIndirect("Eigenvector [Y]", DoubleArray, iSize, iSize); /* [g] : Use Householder Transformation to convert [spStiff] to Triangular form */ matHouseholder( spAwork, spEigenvalue, spEigenvectorY ); /* [h] : Compute [L]^T.[X] = [Y] via Back Substitution */ for (ij = spAwork->iNoRows; ij >= 1 ; ij = ij - 1) { for (ii = 1; ii <= spAwork->iNoColumns; ii = ii + 1) { sum = 0.0; for (ik = ij + 1; ik <= spAwork->iNoRows; ik = ik + 1) sum += spLLT->uMatrix.dpp[ ij-1 ][ ik-1 ] * spEigenvector->uMatrix.dpp[ ik-1 ][ii-1]; spEigenvector->uMatrix.dpp[ij-1][ii-1] = (spEigenvectorY->uMatrix.dpp[ ij-1 ][ ii-1 ] - sum) / spLLT->uMatrix.dpp[ ij-1 ][ ij-1 ]; } } /* [i] : Cleanup before leaving */ matFree( spAwork ); matFree( spCwork ); matFree( spCwork2 ); matFree( spEigenvectorY ); matFree( spLLT ); } /* * =========================================================================== * matHouseholder() : Compute Eigenvalues/Eigenvectors of Symmetric Matrix Via * Householder's Transformation and QL-Algorithm. * * Input : MATRIX *spK -- Pointer to large matrix [K]. * : MATRIX *spEigenvalue -- Pointer to (nx1) eigenvalue matrix. * : MATRIX *spEigenvector -- Pointer to (nxn) eigenvector matrix. * Output : MATRIX *spEigenvalue -- Pointer to (nx1) eigenvalue matrix. * : MATRIX *spEigenvector -- Pointer to (nxn) eigenvector matrix. * * Written By: Mark Austin November 1993 * =========================================================================== */ matHouseholder( spK, spEigenvalue, spEigenvector ) MATRIX *spK; MATRIX *spEigenvalue; MATRIX *spEigenvector; { double *t1, *t2; int stat; int ii, ij; int iNoEquations; /* [a] : Setup Working Matrices for Eigenvalues and Eigenvectors */ iNoEquations = spK->iNoRows; t1 = (double *) safeCalloc( iNoEquations, sizeof(double), __FILE__, __LINE__ ); t2 = (double *) safeCalloc( iNoEquations, sizeof(double), __FILE__, __LINE__ ); for( ii = 1; ii <= iNoEquations; ii = ii + 1) for( ij = 1; ij <= iNoEquations; ij = ij + 1) spEigenvector->uMatrix.dpp[ ii-1 ][ ij-1 ] = spK->uMatrix.dpp[ ii-1 ][ ij-1 ]; /* [b] : Convert [K] to Tridiagonal Form with Householder Transformation */ matHouseholderTransformation( spEigenvector, t1, t2, iNoEquations); /* [c] : Use QL algorithm to compute Eigenvalues and Eigenvectors */ stat = matQLalgorithm( spEigenvector, t1, t2, iNoEquations); /* [d] : Transfer Working Matrix to Eigenvalues */ for( ii = 1; ii <= iNoEquations; ii = ii + 1) spEigenvalue->uMatrix.dpp[ ii-1 ][0] = t1[ ii-1 ]; free((char *) t1); free((char *) t2); return stat; } /* * ================================================================= * matHouseholderTrans() : Use Householder Transformation to convert * Symmetric Matrix to tridiagonal form. * * Input : * Output : * ================================================================= */ static double TOLERANCE = 1.3e-16; static double EPS = 1.3e-16; matHouseholderTransformation( spK , d, e, iNoEquations) MATRIX *spK; double *d, *e; int iNoEquations; { int ii, ij, ik, im; double dTempf, dTempg; double dTemph, dTempk; ii = iNoEquations; while (--ii) { im = ii - 2; dTempf = spK->uMatrix.dpp[ii][ii-1]; dTempg = 0.0; for (ij = 0; ij <= im; ij++) dTempg += spK->uMatrix.dpp[ii][ij] * spK->uMatrix.dpp[ii][ij]; dTemph = dTempg + dTempf * dTempf; if (dTempg < TOLERANCE) { e[ii] = dTempf; dTemph = 0.0; } else { im++; if (dTempf >= 0) dTempg = -sqrt( dTemph ); else dTempg = sqrt( dTemph ); e[ii] = dTempg; dTemph -= dTempf * dTempg; spK->uMatrix.dpp[ii][ii-1] = dTempf - dTempg; dTempf = 0; for (ij = 0; ij <= im; ij++) { spK->uMatrix.dpp[ij][ii] = spK->uMatrix.dpp[ii][ij] / dTemph; dTempg = 0; for (ik = 0; ik <= ij; ik++) dTempg += spK->uMatrix.dpp[ij][ik] * spK->uMatrix.dpp[ii][ik]; for (ik = ij+1; ik <= im; ik++) dTempg += spK->uMatrix.dpp[ik][ij] * spK->uMatrix.dpp[ii][ik]; e[ij] = dTempg / dTemph; dTempf += dTempg * spK->uMatrix.dpp[ij][ii]; } dTempk = dTempf / (dTemph + dTemph); for (ij = 0; ij <= im; ij++) { dTempf = spK->uMatrix.dpp[ii][ij]; dTempg = e[ij] - dTempk * dTempf; e[ij] = dTempg; for (ik=0; ik <= ij; ik++) spK->uMatrix.dpp[ij][ik] = spK->uMatrix.dpp[ij][ik] - dTempf * e[ik] - dTempg * spK->uMatrix.dpp[ii][ik]; } } d[ii] = dTemph; } d[0] = 0.0; e[0] = 0.0; for (ii=0; ii < iNoEquations; ii++) { im = ii - 1; if (d[ii] != 0.0) for (ij=0; ij <= im; ij++) { dTempg = 0.0; for (ik=0; ik <= im; ik++) dTempg += spK->uMatrix.dpp[ii][ik] * spK->uMatrix.dpp[ik][ij]; for (ik=0; ik <= im; ik++) spK->uMatrix.dpp[ik][ij] -= dTempg * spK->uMatrix.dpp[ik][ii]; } d[ii] = spK->uMatrix.dpp[ii][ii]; spK->uMatrix.dpp[ii][ii] = 1.0; for (ij=0; ij <= im; ij++) spK->uMatrix.dpp[ii][ij] = spK->uMatrix.dpp[ij][ii] = 0.0; } } /* * ============================================================== * matQLalgorithm() : QL Algorithm ... adapted from EISPACK...... * * Input : * Output : * ============================================================== */ matQLalgorithm( spK, d, e, iNoEquations ) MATRIX *spK; double *d, *e; int iNoEquations; { int ii, ij, ik, im, in; int imm, iTemp; int flag; double b,c,f,g,h,p,r,s; for (ii = 2; ii <= iNoEquations; ii = ii + 1) e[ ii-2 ] = e[ ii-1 ]; e[ iNoEquations - 1 ] = 0.0; b = 0.0; f = 0.0; for (im = 1; im <= iNoEquations; im = im + 1) { ij = 0; h = EPS * (fabs(d[im-1]) + fabs(e[im-1])); if (b < h) b = h; for(in = im; in <= iNoEquations; in++) if(fabs(e[ in-1 ]) <= b) break; if (in != im) { flag = 1; while (flag) { if (ij == 30) return 1; ij = ij + 1; imm = im + 1; g = d[ im-1 ]; p = (d[ imm - 1 ] - g)/(e[ im-1 ] + e[ im-1 ]); r = sqrt((double) (p*p + 1.0)); if (p < 0.0) d[ im-1 ] = e[ im-1 ] / (p - r); else d[ im-1 ] = e[ im-1 ] / (p + r); h = g - d[ im-1 ]; for(ii = imm; ii <= iNoEquations; ii = ii + 1) d[ ii-1 ] = d[ ii-1 ] - h; f = f + h; p = d[ in-1 ]; c = 1.0; s = 0.0; for (ii = 1; ii <= (in - im); ii = ii + 1) { iTemp = in - ii; g = c * e[( iTemp - 1 )]; h = c * p; if (fabs(p) >= fabs(e[ iTemp-1 ])) { c = e[ iTemp - 1 ] / p; r = sqrt((double) (c * c + 1.0)); e[ iTemp ] = s * p * r; s = c / r; c = 1.0 / r; } else { c = p / e[ iTemp - 1 ]; r = sqrt((double )(c * c + 1.0)); e[ iTemp ] = s * e[ iTemp - 1 ] * r; s = 1.0 / r; c = c * s; } p = c * d[ iTemp - 1 ] - s * g; d[ iTemp ] = h + s * (c * g + s * d[ iTemp - 1 ]); for (ik = 1; ik <= iNoEquations; ik = ik + 1) { h = spK->uMatrix.dpp[ ik-1 ][ iTemp ]; spK->uMatrix.dpp[ ik-1 ][ iTemp ] = s * spK->uMatrix.dpp[ ik-1 ][ iTemp-1 ] + c * h; spK->uMatrix.dpp[ ik-1 ][ iTemp - 1] = c * spK->uMatrix.dpp[ ik-1 ][ iTemp-1 ] - s * h; } } e[ im-1 ] = s * p; d[ im-1 ] = c * p; if(fabs(e[ im-1 ]) <= b) flag = 0; } } d[ im-1 ] = d[ im-1 ] + f; } for (ii = 2; ii <= iNoEquations; ii = ii + 1) { iTemp = ii - 1; ik = iTemp; p = d[ iTemp-1 ]; for (ij = ii; ij <= iNoEquations; ij = ij + 1) { if (fabs(d[ ij-1 ]) > fabs(p)) { ik = ij; p = d[ ij - 1 ]; } } if (ik != iTemp) { d[ ik-1 ] = d[ iTemp-1 ]; d[ iTemp-1 ] = p; for (ij = 1; ij <= iNoEquations; ij = ij + 1) { p = spK->uMatrix.dpp[ ij-1 ][ iTemp-1 ]; spK->uMatrix.dpp[ ij-1 ][ iTemp-1 ] = spK->uMatrix.dpp[ ij-1 ][ ik-1 ]; spK->uMatrix.dpp[ ij-1 ][ ik-1 ] = p; } } } return 0; }