/* * ======================================================================= * Test Program for Subspace Algorithm to solve [K][x] = [Lambda][M][x] * * Assume (a) : [K] and [M] are (nxn) matrices stored with Skyline Pattern * : compute first m eigenvalues/eigenvectors. * ======================================================================= * 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 November 1993 * ======================================================================= */ #include #include #include "matrix.h" #include "vector.h" #include "miscellaneous.h" #include "defs.h" /* #define MYDEBUG */ /* #define MYDEBUG2 */ /* Program defines and error message buffers */ enum { ElementNdof = 2, NoElements = 2, NoEigenValues = 1 }; int main( void ) { MATRIX *spEigenvalue; MATRIX *spEigenvector; MATRIX *spStiff; MATRIX *spMass; MATRIX *spDest; MATRIX *spKe, *spMe; VECTOR *spProfile; int iRowNo; int iColumnNo; int iColumnHeight; int iMin, iMax; int iMinRow; int iSize; int ii, ij, ik; /* [a] : Instantiate and Print Element Level Mass and Stiffness Matrices */ spMe = ElementMassMatrix(); spKe = ElementStiffnessMatrix(); matPrint( spMe ); matPrint( spKe ); /* [b] : Setup and Print Destination Array */ spDest = DestinationMatrix(); matPrint( spDest ); /* [c] : Compute Size and Profile of Skyline Matrix */ iSize = 0; for( ii = 1; ii <= spDest->iNoRows; ii = ii + 1) { for( ij = 1; ij <= spDest->iNoColumns; ij = ij + 1) { iSize = MAX(iSize, spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]); } } spProfile = vecAlloc( "Skyline Profile", IntegerArray, iSize); for( ii = 1; ii <= (int) NoElements; ii = ii + 1) { iMinRow = spDest->uMatrix.ipp[ ii-1 ][ 0 ]; for( ij = 2; ij <= (int) ElementNdof; ij = ij + 1) { if( spDest->uMatrix.ipp[ ii-1 ][ ij-1 ] != 0) iMinRow = MIN(iMinRow, spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]); } for( ij = 1; ij <= (int) ElementNdof; ij = ij + 1) { iColumnNo = spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]; if((1 + iColumnNo - iMinRow) > spProfile->uVector.ip[ iColumnNo-1 ]) spProfile->uVector.ip[ iColumnNo-1 ] = 1 + iColumnNo - iMinRow; } } vecPrint( spProfile ); /* [d] : Assemble Global Mass and Stiffness Matrices */ spStiff = matAllocSkyline("Global Stiffness", DoubleArray, iSize, iSize, spProfile->uVector.ip ); spMass = matAllocSkyline("Global Mass", DoubleArray, iSize, iSize, spProfile->uVector.ip ); for(ii = 1; ii <= (int) NoElements; ii = ii + 1) { for( ij = 1; ij <= (int) ElementNdof; ij = ij + 1) { iRowNo = spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]; if (iRowNo != 0) { for( ik = ij; ik <= (int) ElementNdof; ik = ik + 1) { iColumnNo = spDest->uMatrix.ipp[ ii-1 ][ ik-1 ]; if (iColumnNo != 0) { iMin = MIN( iRowNo, iColumnNo); iMax = MAX( iRowNo, iColumnNo); spStiff->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] += spKe->uMatrix.dpp[ ij-1 ][ ik-1 ]; spMass->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] += spMe->uMatrix.dpp[ ij-1 ][ ik-1 ]; } } } } } /* matPrint( spStiff ); matPrint( spMass ); */ /* [e] : Allocate Memory for Eigenvalues and Eigenvectors */ iSize = spStiff->iNoRows; spEigenvalue = matAllocIndirect( "Eigenvalues", DoubleArray, NoEigenValues, 1 ); spEigenvector = matAllocIndirect("Eigenvectors", DoubleArray, iSize, NoEigenValues); /* [f] : Use Subspace Iteration to Compute Eigenvalues and Eigenvectors */ matSubspace( spStiff, spMass, spEigenvalue, spEigenvector , NoEigenValues); /* [g] : Print EigenValues and EigenVectors */ matPrint( spEigenvalue ); matPrint( spEigenvector ); /* [h] : Clean up before leaving */ vecFree( spProfile ); matFree( spMass ); matFree( spStiff ); matFree( spEigenvector ); matFree( spEigenvalue ); matFree( spDest ); matFree( spMe ); matFree( spKe ); return; } /* * ===================================================== * Setup (2x2) element level Mass and Stiffness Matrices * ===================================================== */ #ifdef __STDC__ MATRIX *ElementMassMatrix() #else /* Start case not STDC */ MATRIX *ElementMassMatrix() #endif /* End case not STDC */ { MATRIX *spMe; /* [a] : Setup (2x2) Element Mass Matrix */ spMe = matAllocIndirect("Element Mass [Me]", DoubleArray, (int) ElementNdof, (int) ElementNdof); switch( (int) ElementNdof ) { case 2: spMe->uMatrix.dpp[ 0 ][ 0 ] = 2.0; spMe->uMatrix.dpp[ 1 ][ 0 ] = 1.0; spMe->uMatrix.dpp[ 0 ][ 1 ] = 1.0; spMe->uMatrix.dpp[ 1 ][ 1 ] = 2.0; break; default: break; } return ( spMe ); } #ifdef __STDC__ MATRIX *ElementStiffnessMatrix() #else /* Start case not STDC */ MATRIX *ElementStiffnessMatrix() #endif /* End case not STDC */ { MATRIX *spKe; /* [a] : Setup (2x2) Element Mass Matrix */ spKe = matAllocIndirect("Element Stiffness [Ke]", DoubleArray, ElementNdof, ElementNdof); switch( ElementNdof ) { case 2: spKe->uMatrix.dpp[ 0 ][ 0 ] = 1000; spKe->uMatrix.dpp[ 1 ][ 0 ] = -1000; spKe->uMatrix.dpp[ 0 ][ 1 ] = -1000; spKe->uMatrix.dpp[ 1 ][ 1 ] = 1000; break; default: break; } return ( spKe ); } /* * ======================== * Setup Destination Matrix * ======================== */ #ifdef __STDC__ MATRIX *DestinationMatrix() #else /* Start case not STDC */ MATRIX *DestinationMatrix() #endif /* End case not STDC */ { MATRIX *spDest; int ii; /* [a] : Setup Destination Array */ spDest = matAllocIndirect("Destination Array", IntegerArray, (int) NoElements, (int) 2); for( ii = 1; ii <= (int) NoElements; ii = ii + 1) { spDest->uMatrix.ipp[ ii-1 ][ 0 ] = (int) (ii % (NoElements + 1)); spDest->uMatrix.ipp[ ii-1 ][ 1 ] = (int) ((ii + 1) % (NoElements + 1)); } return ( spDest ); } /* * ======================================================================== * MatrixAssembleSkyline() : Assemble Skyline Matrix from Indirect Matrices * * Input : MATRIX *spA -- Pointer to large skyline matrix. * : MATRIX *spB -- Pointer to small matrix. * : int *iaRowDest -- integer destination array for rows. * : int *iaColDest -- integer destination array for columns. * Output : MATRIX *spA -- Pointer to large matrix. * ======================================================================== */ #ifdef __STDC__ MATRIX *MatrixAssembleSkyline( MATRIX *spA, MATRIX *spB, int *iaRowDest, int *iaColDest ) #else /* Start case not STDC */ MATRIX *MatrixAssembleSkyline( spA, spB, iaRowDest, iaColDest ) MATRIX *spA; MATRIX *spB; int *iaRowDest; int *iaColDest; #endif /* End case not STDC */ { int ii, ij; int iRowNo; int iColNo; /* [a] : Loop over Row Numbers in spB */ for(ii = 1; ii <= spB->iNoRows; ii = ii + 1) { iRowNo = iaRowDest[ ii-1 ]; if(iRowNo > 0 && iRowNo <= spA->iNoRows ) { /* [b] : Loop over Column Numbers in spB */ for(ij = 1; ij <= spB->iNoColumns; ij = ij + 1) { iColNo = iaColDest[ ij-1 ]; if(iColNo > 0 && iColNo <= spA->iNoColumns ) spA->uMatrix.dpp[ iRowNo-1 ][ iColNo-1 ] += spB->uMatrix.dpp[ ii-1 ][ ij-1 ]; } } } return ( spA ); }