/*
 *  ======================================================================
 *  Program to Test Computation of Matrix Triple Product
 *  ======================================================================
 *  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 1993
 *  ======================================================================
 */

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

int main( void ) {
enum { ProblemSize = 3 };
MATRIX *matTripleProduct1( MATRIX * , MATRIX * , MATRIX * );
MATRIX *matTripleProduct2( MATRIX * , MATRIX * , MATRIX * );
MATRIX *spA, *spB, *spC, *spTp;
int ii, ij;

   /* [a] : Allocate Memory for Small 5 x 5 Matrices */

   spA = matAllocIndirect( "A", DoubleArray, ProblemSize, ProblemSize );
   spB = matAllocIndirect( "B", DoubleArray, ProblemSize, ProblemSize );
   spC = matAllocIndirect( "C", DoubleArray, ProblemSize, 1 );

   /* [b] : Initialize Values in Matrix */

   for( ii = 1; ii <= ProblemSize; ii++ ) {
        spC->uMatrix.dpp[ ii - 1 ][ 0 ] = ii;
        for( ij = 1; ij <= ProblemSize; ij++ )  {
             spA->uMatrix.dpp[ ii - 1 ][ ij - 1 ] = ii - ij;
             spB->uMatrix.dpp[ ii - 1 ][ ij - 1 ] = ii - ij;
        }
   }

   /* [c] : Print Matrices "A", "B" and "C"  */

   matPrint( spA );
   matPrint( spB );
   matPrint( spC );

   /* [d] : Compute and Print Matrix Triple Products */

   spTp         = matTripleProduct1( spA, spB, spC );
   spTp->cpName = saveString("Triple Product 1", __FILE__, __LINE__ );
   matPrint(spTp);

   spTp         = matTripleProduct2( spA, spB, spC );
   spTp->cpName = saveString("Triple Product 2", __FILE__, __LINE__ );
   matPrint(spTp);

   /* [e] : Free Matrix Memory */

   matFree( spA );
   matFree( spB );
   matFree( spC );
}

/*
 *  =====================================================================
 *  matTripleProduct1() : Use 3-level looping construct to compute matrix
 *                        triple product.
 * 
 *  Input  : MATRIX * spA = pointer to (nxn) matrix.
 *         : MATRIX * spB = pointer to (nxn) matrix.
 *         : MATRIX * spC = pointer to (nx1) matrix.
 *  Output : MATRIX *     = pointer to (nx1) triple product matrix.
 *  =====================================================================
 */

MATRIX *matTripleProduct1( MATRIX *spA, MATRIX *spB, MATRIX *spC ) {
MATRIX   *spTp;
int ii, ij, ik;

   spTp = matAllocIndirect( "Triple Product", DoubleArray, spC->iNoRows, 1 );

   for(ii = 0; ii < spA->iNoRows; ii++)
   for(ik = 0; ik < spA->iNoColumns; ik++)
   for(ij = 0; ij < spC->iNoRows; ij++)
       spTp->uMatrix.dpp[ ii ][0] += ( spA->uMatrix.dpp[ii][ik] *
                                       spB->uMatrix.dpp[ik][ij] *
                                       spC->uMatrix.dpp[ij][0] );

   return ( spTp );
}

/*
 *  =====================================================================
 *  matTripleProduct2() : Use pair of matrix multiplies to compute matrix
 *                        triple product.
 * 
 *  Input  : MATRIX * spA = pointer to (nxn) matrix.
 *         : MATRIX * spB = pointer to (nxn) matrix.
 *         : MATRIX * spC = pointer to (nx1) matrix.
 *  Output : MATRIX *     = pointer to (nx1) triple product matrix.
 *  =====================================================================
 */

MATRIX *matTripleProduct2( MATRIX *spA, MATRIX *spB, MATRIX *spC ) {
MATRIX   *spTp;
MATRIX *spTemp;

   spTemp = matMult( spB, spC    );
   spTp   = matMult( spA, spTemp );
   matFree( spTemp );

   return ( spTp );
}