/* * ====================================================================== * Test program for Naive Implementation of Gauss Elimination * * Note : A.x = b corresponds to displacement equations is a * four story shear structure. * * ====================================================================== * Copyright (C) 1996 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 April 1996 * ====================================================================== */ #include #include "matrix.h" #include "vector.h" #include "miscellaneous.h" /* #define NAIVE_GAUSS */ #define SWAP_ROWS VECTOR *NaiveGaussElimination( MATRIX *, VECTOR *); VECTOR *GaussElimination( MATRIX *, VECTOR *); VECTOR *SetupScaleFactors( MATRIX * ); VECTOR *SetupPivotVector( MATRIX * ); int main( void ) { enum { NoFloors = 4 }; MATRIX *spStiff; VECTOR *spDispl, *spEload; double *pdTemp, dTemp; int iSecondRow, iFirstRow; int ii, ij; /* [a] : Allocate Matrices for Stiffness and External Load */ spStiff = matAllocIndirect("Stiffness", DoubleArray, NoFloors, NoFloors ); spEload = vecAlloc("External Load", DoubleArray, NoFloors ); /* [b] : Set values in external Load and Stiffness Matrices */ for(ii = 1; ii <= NoFloors; ii++) spEload->uVector.dp[ ii - 1 ] = 10.0*((double) (ii) / ((double) NoFloors)); spStiff->uMatrix.dpp[ 0 ][ 0 ] = 2000; spStiff->uMatrix.dpp[ 0 ][ 1 ] = -1000; spStiff->uMatrix.dpp[ 1 ][ 0 ] = -1000; spStiff->uMatrix.dpp[ 1 ][ 1 ] = 2000; spStiff->uMatrix.dpp[ 1 ][ 2 ] = -1000; spStiff->uMatrix.dpp[ 2 ][ 1 ] = -1000; spStiff->uMatrix.dpp[ 2 ][ 2 ] = 2000; spStiff->uMatrix.dpp[ 2 ][ 3 ] = -1000; spStiff->uMatrix.dpp[ 3 ][ 2 ] = -1000; spStiff->uMatrix.dpp[ 3 ][ 3 ] = 1000; /* [c] : Test algorithm by swapping two rows */ #ifdef SWAP_ROWS iFirstRow = 1; iSecondRow = 4; printf("*** Swap rows %3d <==> %3d\n", iFirstRow, iSecondRow); pdTemp = spStiff->uMatrix.dpp[ iFirstRow - 1 ]; spStiff->uMatrix.dpp[ iFirstRow - 1 ] = spStiff->uMatrix.dpp[ iSecondRow-1 ]; spStiff->uMatrix.dpp[ iSecondRow-1 ] = pdTemp; dTemp = spEload->uVector.dp[ iFirstRow - 1 ]; spEload->uVector.dp[ iFirstRow - 1 ] = spEload->uVector.dp[ iSecondRow - 1 ]; spEload->uVector.dp[ iSecondRow- 1 ] = dTemp; #endif /* [d] : print stiffness/external load matrices */ matPrint( spStiff ); vecPrint( spEload ); /* [e] : Solve equations and print displacements */ #ifdef NAIVE_GAUSS spDispl = NaiveGaussElimination( spStiff, spEload ); spDispl->cpName = saveString("Displacement", __FILE__, __LINE__ ); #else spDispl = GaussElimination( spStiff, spEload); spDispl->cpName = saveString("Displacement", __FILE__, __LINE__ ); #endif vecPrint( spDispl ); } /* * ================================================================ * Use Gauss Elimination with Scaling and Partial Pivoting to solve * linear equations [A][X] = [B]. [A] is a (nxn) matrix. * * Input : MATRIX spA -- Pointer to (nxn) matrix A. * VECTOR spB -- Pointer to (nx1) vector B. * Output : VECTOR spX -- Pointer to Solution Vector spX. * ================================================================ */ #ifdef __STDC__ VECTOR *GaussElimination( MATRIX *spA, VECTOR *spB) #else /* Start case not STDC */ VECTOR *GaussElimination( spA, spB ) MATRIX *spA; VECTOR *spB; #endif /* End case not STDC */ { VECTOR *spScale, *spPivot, *spX; double multiplier, sum; double dNumerator, dDenominator; int iOrder1, iOrder2; int ii, ij, ik; /* [a] : Check compatibility of matrix sizes */ if((spA->iNoRows != spA->iNoColumns) || (spA->iNoRows != spB->iLength)) { printf("ERROR >> spA->iNoRows = %4d spA->iNoColumns = %4d\n", spA->iNoRows, spA->iNoColumns); printf("ERROR >> spB->iLength = %4d spB->iNoColumns = %4d\n", spB->iLength, 1); fatalError("ERROR >> Execution halted in GaussElimination()", "ERROR >> Matrix Dimensions are Inconsistent Dimensions", (char *) NULL); } /* [b] : Allocate Memory for Displacement Vector */ spX = vecAlloc( (char *) NULL, DoubleArray, spA->iNoRows ); /* [c] : Setup Scale Factors and Order Vector for Partial Pivoting */ spScale = SetupScaleFactors( spA ); spPivot = SetupPivotVector( spA ); vecPrint( spScale ); vecPrint( spPivot ); /* [d] : Elimination of Row Elements */ for (ii = 1; ii <= spA->iNoColumns - 1; ii = ii + 1) { Pivot( spA, spScale, spPivot, ii ); iOrder1 = spPivot->uVector.ip[ ii - 1 ]; for (ij = ii + 1; ij <= spA->iNoRows; ij = ij + 1) { iOrder2 = spPivot->uVector.ip[ ij - 1 ]; dNumerator = spA->uMatrix.dpp[ iOrder2-1 ][ ii-1 ]; dDenominator = spA->uMatrix.dpp[ iOrder1-1 ][ ii-1 ]; multiplier = dNumerator/dDenominator; for (ik = ii + 1; ik <= spA->iNoRows; ik = ik + 1) { spA->uMatrix.dpp[iOrder2-1][ik - 1] -= multiplier * spA->uMatrix.dpp[iOrder1-1][ik - 1]; } spB->uVector.dp[ iOrder2-1 ] -= multiplier * spB->uVector.dp[ iOrder1-1 ]; } } /* [e] : Compute Backsubstitution */ ii = spA->iNoColumns; iOrder2 = spPivot->uVector.ip[ ii - 1 ]; spX->uVector.dp[ii-1] = spB->uVector.dp[ iOrder2-1 ] / spA->uMatrix.dpp[ iOrder2-1 ][ ii-1 ]; for (ii = spA->iNoColumns - 1; ii >= 1; ii = ii - 1) { iOrder1 = spPivot->uVector.ip[ ii - 1 ]; sum = 0.0; for (ij = ii + 1; ij <= spA->iNoColumns; ij = ij + 1) sum += spA->uMatrix.dpp[ iOrder1-1 ][ij-1]*spX->uVector.dp[ij-1]; spX->uVector.dp[ii-1] = (spB->uVector.dp[iOrder1-1] - sum) / spA->uMatrix.dpp[iOrder1-1][ii-1]; } vecPrint( spScale ); vecPrint( spPivot ); return( spX ); } /* * =================================================================== * NaiveGaussElimination() : Gauss Elimiation without partial pivoting * or scaling of equations. * * Input : MATRIX spA -- Pointer to (nxn) matrix [A]. * : VECTOR spB -- Pointer to (nx1) vector [B]. * Output : MATRIX spX -- Pointer to (nx1) solution vector [X]. * =================================================================== */ #ifdef __STDC__ VECTOR *NaiveGaussElimination( MATRIX *spA, VECTOR *spB ) #else /* Start case not STDC */ VECTOR *NaiveGaussElimination( spA, spB ) MATRIX *spA; VECTOR *spB; #endif /* End case not STDC */ { int ii, ij, ik; double multiplier; double sum; VECTOR *spX; /* [a] : Elimination of row elements */ for (ii = 1; ii <= spA->iNoColumns - 1; ii = ii + 1) { for (ij = ii + 1; ij <= spA->iNoRows; ij = ij + 1) { multiplier = spA->uMatrix.dpp[ij-1][ii-1]/spA->uMatrix.dpp[ii-1][ii-1]; for (ik = ii; ik <= spA->iNoRows; ik = ik + 1) { spA->uMatrix.dpp[ij -1][ik - 1] -= multiplier*spA->uMatrix.dpp[ii -1][ik - 1]; } spB->uVector.dp[ij-1] -= multiplier*spB->uVector.dp[ii-1]; } } /* [b] : Allocate Memory for Displacement Vector */ spX = vecAlloc( (char *) NULL, DoubleArray, spA->iNoRows ); /* [c] : Compute Backsubstitution */ ii = spA->iNoColumns; spX->uVector.dp[ii-1] = spB->uVector.dp[ii-1] / spA->uMatrix.dpp[ii-1][ii-1]; for (ii = spA->iNoColumns - 1; ii >= 1; ii = ii - 1) { sum = 0.0; for (ij = ii + 1; ij <= spA->iNoColumns; ij = ij + 1) sum += spA->uMatrix.dpp[ii-1][ij-1]*spX->uVector.dp[ij-1]; spX->uVector.dp[ii-1] = (spB->uVector.dp[ii-1] - sum)/spA->uMatrix.dpp[ii-1][ii-1]; } return ( spX ); }