/*
 *  ==============================================================================
 *  Use Householder Transformation and QL-Algorithm to compute Eigenvalues/vectors
 *  of (5x5) symmetric matrix.  
 *  ==============================================================================
 *  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                                              June 1996 
 *  ==============================================================================
 */

#include <stdio.h>
#include <math.h>
#include "matrix.h"
#include "vector.h"
#include "miscellaneous.h"

enum { NoNodes = 5 };

int main( void ) {
MATRIX *spMass, *spStiff;
MATRIX     *spEigenvalue;
MATRIX    *spEigenvector;
int       ii, ij,  iSize;

   /* [a] : Allocate Memory for Matrices */

   spStiff = matAllocIndirect( "Stiffness",  DoubleArray, NoNodes, NoNodes );
   spMass  = matAllocIndirect(      "Mass",  DoubleArray, NoNodes, NoNodes );

   /* [b] : Fill-In Contents of Stiffness and Mass Matrices */

   for(ii = 1; ii <= NoNodes; ii++) {
       spMass->uMatrix.dpp[ ii-1 ][ ii-1 ] = 1.0;
       for(ij = 1; ij <= NoNodes; ij++) {
           if((ii - ij) == 1 || (ii - ij) == -1)
               spStiff->uMatrix.dpp[ ii-1 ][ ij-1 ] = -1000;
           else if(ii == ij && (ii == 1))
               spStiff->uMatrix.dpp[ ii-1 ][ ij-1 ] =  1000;
           else if(ii == ij && (ii != 1))
               spStiff->uMatrix.dpp[ ii-1 ][ ij-1 ] =  2000;
           else
               spStiff->uMatrix.dpp[ ii-1 ][ ij-1 ] =  0;
       }
   }

   /* [c] : Print Stiffness and Mass Matrices */

   matPrint( spStiff );
   matPrint( spMass  );

   /* [d] : Allocate Memory for Eigenvalues and Eigenvectors */

   iSize = spStiff->iNoRows;
   spEigenvalue  = matAllocIndirect( "Eigenvalues", DoubleArray, iSize,     1);
   spEigenvector = matAllocIndirect("Eigenvectors", DoubleArray, iSize, iSize);

   /* [e] : Use Householder Transformation to convert [spStiff] to Triangular form */

   matHouseholder( spStiff, spEigenvalue, spEigenvector );

   /* [f] : Print Eigenvalues and Eigenvectors */

   matPrint( spEigenvalue  );
   matPrint( spEigenvector );
}