/*
 *  ======================================================================
 *  Use LU decomposition with row permutations to solve [A][X] = [B].
 *
 *  Note : [A].[x] = [b] corresponds to displacement equations in a four
 *         story shear structure.
 *  ======================================================================
 *  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                                     July 1996
 *  ======================================================================
 */

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

#define SWAP_ROWS
enum { NoFloors = 4 };

int main( void ) {
MATRIX  *spStiff, *spEload;
VECTOR  *spPivot;
double  *pdTemp, dTemp;
int iFirstRow, iSecondRow;
int ii, ij;

   /* [a] : Allocate Matrices for Stiffness and External Load */

   spStiff = matAllocIndirect(    "Stiffness", DoubleArray, NoFloors, NoFloors );
   spEload = matAllocIndirect("External Load", DoubleArray, NoFloors, 1 );

   /* [b] : Set values in external Load and Stiffness Matrices */

   for(ii = 1; ii <= NoFloors; ii++) 
       spEload->uMatrix.dpp[ ii - 1 ][0] = 10.0*((double) (ii) / ((double) NoFloors));

   spStiff->uMatrix.dpp[ 0 ][ 0 ] =  2000;
   spStiff->uMatrix.dpp[ 0 ][ 1 ] = -1000;

   spStiff->uMatrix.dpp[ 1 ][ 0 ] = -1000;
   spStiff->uMatrix.dpp[ 1 ][ 1 ] =  2000;
   spStiff->uMatrix.dpp[ 1 ][ 2 ] = -1000;

   spStiff->uMatrix.dpp[ 2 ][ 1 ] = -1000;
   spStiff->uMatrix.dpp[ 2 ][ 2 ] =  2000;
   spStiff->uMatrix.dpp[ 2 ][ 3 ] = -1000;

   spStiff->uMatrix.dpp[ 3 ][ 2 ] = -1000;
   spStiff->uMatrix.dpp[ 3 ][ 3 ] =  1000;

   /* [c] : Test algorithm by swapping two rows */

#ifdef SWAP_ROWS
   iFirstRow  = 1;
   iSecondRow = 4;

   printf("*** Swap rows %3d <==> %3d\n", iFirstRow, iSecondRow);

   pdTemp = spStiff->uMatrix.dpp[ iFirstRow - 1 ];
   spStiff->uMatrix.dpp[ iFirstRow - 1 ] = spStiff->uMatrix.dpp[ iSecondRow-1 ];
   spStiff->uMatrix.dpp[ iSecondRow-1 ]  = pdTemp;

   dTemp = spEload->uMatrix.dpp[ iFirstRow - 1 ][0];
   spEload->uMatrix.dpp[ iFirstRow - 1 ][0] = spEload->uMatrix.dpp[ iSecondRow - 1 ][0];
   spEload->uMatrix.dpp[ iSecondRow- 1 ][0] = dTemp;
#endif

   /* [d] : print stiffness/external load matrices */

   matPrint( spStiff );
   matPrint( spEload );

   /* [d] : Compute LU decomposition and forward/backward substitution */

   spPivot     = SetupPivotVector( spStiff );
   spStiff     = matLUDecompositionIndirect( spStiff, spPivot);

   spEload         = matLUSubstitutionIndirect( spPivot , spStiff, spEload);
   spEload->cpName = saveString("Solution Vector", __FILE__, __LINE__ );

   /* [e] : Print Displacements */

   vecPrint( spPivot );
   matPrint( spStiff );
   matPrint( spEload );
}