/*
 *  ======================================================================
 *  Generate Upper Triangular Matrix and Solve by Backsubstitution
 *
 *  Copyright (C) 1994-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 1994
 *  ======================================================================
 */

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

VECTOR *BackSubstitution( MATRIX *, VECTOR * );

enum { NoNodes = 4 };

int main( void ) {
VECTOR  *spB, *spX;
MATRIX  *spU;
int   ii, ij;

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

    spU = matAllocIndirect(" [U] ", DoubleArray, NoNodes, NoNodes);
    spB =         vecAlloc(" [B] ", DoubleArray, NoNodes);

    /* [b] : Initialize contents of upper-triangular matrix and r.h.s. vector */

    for(ii = 1; ii <= (int) NoNodes; ii++) {
        spB->uVector.dp[ ii - 1 ] = (int) NoNodes - 2*ii;
        for(ij = ii; ij <= (int) NoNodes; ij++) {
            spU->uMatrix.dpp[ ii-1 ][ ij-1 ] = 1;
        }
    }

    /* [c] : Print Matrices */

    matPrint( spU );
    vecPrint( spB );

    /* [d] : Solve equations and print displacements */

    spX = BackSubstitution( spU, spB );

    vecPrint( spX );
}

/* 
 *  =========================================================================
 *  BackSubstitution() : Solve Upper Triangular System via back substitution.
 *  
 *  Input  : MATRIX spU       -- Pointer to (nxn) matrix [U].
 *         : VECTOR spB       -- Pointer to (nx1) vector [B].
 *  Output : VECTOR spX       -- Solution vector is stored in [X].
 *  =========================================================================
 */

VECTOR *BackSubstitution( MATRIX *spA, VECTOR *spB ) {
VECTOR   *spX;
int    ii, ij;
double   dSum;

   /* [a] : Allocate memory for solution vector */

   spX = vecAlloc(" [X] ", DoubleArray, (int) spB->iLength );

   /* [b] : Compute backsubstitution */

   ii = spA->iNoColumns;
   spX->uVector.dp[ii-1] = spB->uVector.dp[ii-1] / spA->uMatrix.dpp[ii-1][ii-1];
   for (ii = spA->iNoColumns - 1; ii >= 1; ii = ii - 1) {
        dSum = 0.0;
        for (ij = ii + 1; ij <= spA->iNoColumns; ij = ij + 1)
             dSum += spA->uMatrix.dpp[ii-1][ij-1]*spX->uVector.dp[ij-1];

        spX->uVector.dp[ii-1] = (spB->uVector.dp[ii-1] - dSum)/spA->uMatrix.dpp[ii-1][ii-1];
   }

   return ( spX );
}