/*
 *  ======================================================================
 *  Program to allocate, print and free a skyline matrix.
 *  ======================================================================
 *  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                                   October 1993
 *  ======================================================================
 */

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

enum { MatrixSize = 6 };

int main( void ) {
MATRIX       *spA;
MATRIX   *spATemp;
MATRIX       *spB;
MATRIX       *spC;
MATRIX      *spAB;
MATRIX      *spLU;
int           *ld;
int        iStart;
int        iRowNo;
int     iColumnNo;
int    iMin, iMax;
int iColumnHeight;
int    iCount = 0;

    /* [a] : Define Profile of Skyline Matrix */
 
    ld = (int *) safeCalloc( MatrixSize, sizeof(int), __FILE__, __LINE__ );

    ld[0] = 1; ld[1] = 2; ld[2] = 3;
    ld[3] = 1; ld[4] = 1; ld[5] = 3;

    /* [b] : Allocate Memory for Skyline Matrix */

    spA = matAllocSkyline( "Skyline", DoubleArray, MatrixSize, MatrixSize, ld );
    spC = matAllocIndirect(    "[C]", DoubleArray, MatrixSize, 2 );

    for( iRowNo = 1; iRowNo <= (int) MatrixSize; iRowNo++)
         for( iColumnNo = 1; iColumnNo <= 2; iColumnNo++ ) 
              spC->uMatrix.dpp[ iRowNo-1 ][ iColumnNo-1 ] = 1.0;

    /* [c] : Initialize Contents of Skyline Matrix */

    for( iColumnNo = 1; iColumnNo <= (int) MatrixSize; iColumnNo++) {

         iColumnHeight = (int) spA->uMatrix.dpp[ iColumnNo-1 ][0];
         iStart = iColumnNo - iColumnHeight + 1;

         for( iRowNo = iStart; iRowNo <= iColumnNo; iRowNo++ ) {
              iCount = iCount + 1;

              iMin = MIN( iRowNo, iColumnNo);
              iMax = MAX( iRowNo, iColumnNo);

              spA->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] = (double) iCount;
         }
   }

   /* [d] : Setup Right-hand Side of Equation */

   spB = matAllocIndirect( "Load Vector [B]", DoubleArray, MatrixSize, 1 );

   spB->uMatrix.dpp[ 0 ][ 0 ] = (double)  6.5;
   spB->uMatrix.dpp[ 1 ][ 0 ] = (double)  9.0;
   spB->uMatrix.dpp[ 2 ][ 0 ] = (double) 13.0;
   spB->uMatrix.dpp[ 3 ][ 0 ] = (double) 16.0;
   spB->uMatrix.dpp[ 4 ][ 0 ] = (double) 18.0;
   spB->uMatrix.dpp[ 5 ][ 0 ] = (double) 30.0;

   /* [e] : Print Skyline Matrix */

   matPrint( spA );
   matPrint( spB );

   /* [f] : Compute [L][D][L]^T Decomposition and Solve Equations */

   spATemp = matCopy( spA );
   spATemp = matLUDecompositionSkyline( spATemp );
   spB     = matLUBacksubstitutionSkyline( spATemp , spB );

   matPrint( spB );

   /* [g] : Test Skyline-Indirect Matrix Multiplication */

   printf("*** Test Functions for Mixed Skyline-Indirect Matrix Multiply \n");

   matPrint( spA );
   matPrint( spC );

   spAB = matMult( spA, spC );

   matPrint( spAB );

   /* [h] : Test Indirect-Skyline Matrix Multiplication */

   /* [i] : Clean up before leaving */

   matFree( spA );
   matFree( spATemp );
   matFree( spB );

}