/*
 *  ======================================================================
 *  Test Program for Inverse Iteration to Solve Symmetric General
 *  Eigenvalue Problem
 *  
 *  [A] [X] = [L] [B] [X] : where [L] = Dominant Eigenvalue
 *                                [X] = Dominant Eigenvector
 *  
 *  ======================================================================
 *  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                                      June 1995
 *  ======================================================================
 */

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

int main( void ) {
MATRIX  *spMass, *spStiff;
MATRIX      *spEigenvalue;
MATRIX     *spEigenvector;
enum { NoDof = 2, One = 1 };

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

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

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

   spStiff->uMatrix.dpp[ 0 ][ 0 ] =  200;
   spStiff->uMatrix.dpp[ 0 ][ 1 ] = -100;
   spStiff->uMatrix.dpp[ 1 ][ 0 ] = -100;
   spStiff->uMatrix.dpp[ 1 ][ 1 ] =  200;

   spMass->uMatrix.dpp[ 0 ][ 0 ] = 1;
   spMass->uMatrix.dpp[ 1 ][ 1 ] = 2;

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

   matPrint( spStiff );
   matPrint( spMass  );

   /* [d] : Allocate Memory for Eigenvalue and Eigenvector Matrix */

   spEigenvector = matAllocIndirect("Eigenvector", DoubleArray, NoDof, One );
   spEigenvalue  = matAllocIndirect( "Eigenvalue", DoubleArray,   One, One );

   /* [e] : Solve [K][X] = [L][M][X] : [L] = Dominant Eigenvalue    */
   /*                                  [X] = Dominant Eigenvector   */

   matInverseIteration( spStiff, spMass , spEigenvalue, spEigenvector );

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

   matPrint( spEigenvalue  );
   matPrint( spEigenvector );
}


/* 
 *  ========================================================================
 *  matInverseIteration() : Use Inverse Iteration to Solve Symmetric General
 *  Eigenvalue Problem
 *  
 *  [K] [X] = [L] [M] [X] : where [L] = Dominant Eigenvalue
 *                                [X] = Dominant Eigenvector
 *  
 *  Input :   MATRIX  *spK           -- pointer to [K] Matrix 
 *            MATRIX  *spM           -- pointer to [M] Matrix 
 *            MATRIX  *spEigenvalue  -- pointer to [L] Matrix 
 *            MATRIX  *spEigenvector -- pointer to [X] Matrix 
 *
 *  Output :  MATRIX  *spEigenvalue  -- pointer to [L] Matrix 
 *            MATRIX  *spEigenvector -- pointer to [X] Matrix 
 *  ========================================================================
 */

static double dTolerance = 0.000001;

matInverseIteration( MATRIX          *spK, MATRIX *spM,
                     MATRIX *spEigenvalue, MATRIX *spEigenvector ) {
MATRIX        *spLU;
MATRIX *spX1, *spXK;
MATRIX *spY1, *spYK;
MATRIX       *spNum;
MATRIX     *spDenom;
VECTOR     *spPivot;
int          ii, ij;
double dConvergence;
double dEigenvalueNew = 0;
double dEigenvalueOld = 0;
enum { One = 1, MaxIterations = 20 };

  /* [a] : Initialize [X1] vector */

  spX1 = matAllocIndirect("Eigenvector", DoubleArray, spM->iNoRows, One );
  for(ii = 1; ii <= spM->iNoRows; ii = ii + 1)
      spX1->uMatrix.dpp[ ii - 1 ][0] = 1.0;

  /* [b] : Compute [Y1] = [M].[X1] */

  spY1 = matMultIndirectDouble( spM, spX1 );

  /* [c] : Iteration Procedure for First Mode Eigenvalue and Eigenvector */

  spPivot = SetupPivotVector( spK );
  spLU    = matLUDecompositionIndirect( spK , spPivot);

  for(ii = 1; ii <= MaxIterations; ii = ii + 1) {

     /* [c.1] : Solve Equations [K][XK] = [Y1] */
    
     spXK = matLUSubstitutionIndirect( spPivot , spLU, spY1 );

     /* [c.2] : Computer [YK] = [M][XK] */

     spYK = matMultIndirectDouble( spM , spXK );

     /* [c.3] : Compute Rayleigh Quotient */

     spNum   = matMultIndirectDouble( matTransposeIndirectDouble( spXK ), spY1 );
     spDenom = matMultIndirectDouble( matTransposeIndirectDouble( spXK ), spYK );

     dEigenvalueNew = spNum->uMatrix.dpp[0][0] / spDenom->uMatrix.dpp[0][0];

     /* [c.4] : Normalize [YK] */

     for(ij = 1; ij <= spM->iNoRows; ij = ij + 1)
         spY1->uMatrix.dpp[ ij-1 ][0] = spYK->uMatrix.dpp[ ij-1 ][0] /
                                        sqrt( spDenom->uMatrix.dpp[0][0] );

     /* [c.5] : Test Convergence */

     dConvergence = fabs(dEigenvalueOld - dEigenvalueNew)/dEigenvalueNew;
     if(dConvergence < dTolerance) {

        printf("\n*** INVERSE ITERATION CONVERGED IN %2d ITERATIONS \n", ii);

        spEigenvalue->uMatrix.dpp[0][0] = dEigenvalueNew;
        for(ij = 1; ij <= spM->iNoRows; ij = ij + 1)
            spEigenvector->uMatrix.dpp[ ij-1 ][0] = spXK->uMatrix.dpp[ ij - 1 ][0] /
                                                    sqrt( spDenom->uMatrix.dpp[0][0] );

        break;
     } else {
        dEigenvalueOld = dEigenvalueNew;
     }
  }
}