/*
 *  ======================================================================
 *  Evaluate geometric series of matrices.
 *  ======================================================================
 *  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 <stdio.h>
#include "matrix.h"
#include "miscellaneous.h"

int main( void ) {
MATRIX *spA, *spI,  *spSum1,  *spSum2, *spResidual;
MATRIX   *spTemp1, *spTemp2, *spTemp3,    *spTemp4;
enum { ProblemSize = 2 , SeriesLength = 15 };
int  ii, ij;

   /* [a] : Allocate and initialize "A" and "I" matrices */

   spA = matAllocIndirect("[A]", DoubleArray, ProblemSize, ProblemSize );
   spI = matAllocIndirect("[I]", DoubleArray, ProblemSize, ProblemSize );
   for (ii = 1; ii <= ProblemSize; ii = ii + 1) {
        spI->uMatrix.dpp[ii-1][ii-1] = 1.0;
        for (ij = 1; ij <= ProblemSize; ij = ij + 1)
             spA->uMatrix.dpp[ii-1][ij-1] = ii + ij;
   }

   matPrint( spI );
   matPrint( spA );

   /* [b] : Use Horner's rule to compute matrix series summation */

   spSum1 = matCopy (spA);
   for (ii = 1; ii <= SeriesLength - 1; ii = ii + 1) {
        spSum1 = matAddReplace( spSum1, spI ); 
        spSum1 = matMultReplace( spSum1, spA ); 
   }
   spSum1 = matAddReplace( spSum1, spI ); 
   spSum1->cpName = saveString("Series Summation 1", __FILE__, __LINE__ );

   matPrint( spSum1);

   /* [c] : Use summation rule to compute matrix series directly */

   spTemp1 = matCopy (spA);
   for (ii = 1; ii <= SeriesLength; ii = ii + 1)
        spTemp1 = matMultReplace( spTemp1, spA ); 
   
   spTemp2 = matSub( spI, spTemp1 ); 
   spTemp3 = matSub( spI, spA ); 
   spTemp4 = matSolve (spTemp3, spI);
   spSum2  = matMult( spTemp2 , spTemp4 );
   spSum2->cpName = saveString("Series Summation 2", __FILE__, __LINE__ );

   matPrint( spSum2);

   /* [d] : Evaluate and print difference in series summations */

   spResidual = matSub ( spSum2, spSum1 );
   spResidual->cpName = saveString("Summation 2 - Summation 1", __FILE__, __LINE__ );
   matPrint( spResidual );

   /* [e] : Clean up and leave program */

   matFree( spA );
   matFree( spI );
   matFree( spSum1 );
   matFree( spSum2 );
   matFree( spTemp1 );
   matFree( spTemp2 );
   matFree( spTemp3 );
   matFree( spResidual );
}