/*
 *  =================================================================== 
 *  High-Level Functions for Matrix Operations and Printing
 *  =================================================================== 
 *  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 1992 - October 1993
 *  ===================================================================
 */

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

/*
 *  ==========================================
 *  matPrint() : Output matrix.
 *  
 *  Input  : MATRIX spA -- pointer to matrix A
 *  Output : void.
 *  ==========================================
 */

#ifdef __STDC__
void matPrint( MATRIX *spA )
#else  /* start case not STDC */
void matPrint( spA )
MATRIX *spA;
#endif /* end case not STDC */
{
       switch((int) spA->Rep) {
           case Sequential:
           case Indirect:
                switch((int) spA->Type) {
                    case DoubleArray:
                         matPrintIndirectDouble( spA );
                         break;
                    case IntegerArray:
                         matPrintIndirectInteger( spA );
                         break;
                    case ComplexArray:
                         fatalError("In matPrint() : spA->Type not implemented",
                                   (char *) NULL);
                         break;
                    default:
                         fatalError("In matPrint() : Undefined spA->Type",
                                   (char *) NULL);
                         break;
                }
                break;
           case Skyline:
                matPrintSkylineDouble( spA );
                break;
           case Sparse:
           default:
                fatalError("In matPrint() : Undefined spA->Rep", (char *) NULL);
                break;
       }
}

/*
 *  ==========================================
 *  matDump() : Dump Matrix.
 *  
 *  Input  : MATRIX spA -- pointer to matrix A
 *  Output : void.
 *  ==========================================
 */

#ifdef __STDC__
void matDump( MATRIX *spA )
#else  /* start case not STDC */
void matDump( spA )
MATRIX *spA;
#endif /* end case not STDC */
{

       switch((int) spA->Rep) {
           case Skyline:
                matDumpSkyline( spA );
                break;
           case Indirect:
           case Sparse:
           default:
                fatalError("In matDump() : Undefined spA->Rep",
                          (char *) NULL);
                break;
       }
}

/*
 *  ===========================================
 *  matFree() : Free Matrix Storage
 *  
 *  Input  : Matrix spA -- pointer to matrix A.  
 *  Output : void.
 *  ===========================================
 */

#ifdef __STDC__
void matFree( MATRIX *spA )
#else  /* start case not STDC */
void matFree( spA )
MATRIX *spA;
#endif /* end case not STDC */
{
      switch((int) spA->Rep) {
          case Indirect:
               matFreeIndirect( spA );
               break;
          case Skyline:
               matFreeSkyline( spA );
               break;
          default:
               fatalError("In matFree() : Undefined spA->Rep",
                         (char *) NULL);
               break;
      }
}


/*
 *  ============================================================
 *  matScale() : Scale a Matrix by a constant [c] = d.[a]
 *  
 *  Input :    MATRIX spA          -- Pointer to matrix A
 *             double dScaleFactor -- scale factor.
 *  Output :   MATRIX spC          -- Pointer to scaled matrix C 
 *  ============================================================
 */

#ifdef __STDC__
MATRIX *matScale(MATRIX *spA, double factor )
#else
MATRIX *matScale(spA, factor )
MATRIX  *spA;
double  factor;
#endif
{
MATRIX *spB;

   switch((int) spA->Rep) {
       case Indirect:
            switch((int) spA->Type) {
                case DoubleArray:
                     spB = matScaleIndirectDouble( spA, factor );
                     break;
                default:
                     fatalError("In matScale() : Undefined spA->eType", (char *) NULL);
                     break;
            }
            break;
       case Skyline:
            spB = matScaleSkyline( spA, factor );
            break;
       default:
            fatalError("In matScale() : Undefined spA->eRep", (char *) NULL);
            break;
   }

   return ( spB );
}


/*
 *  ===================================================
 *  matAdd() : Matrix Add Operation [c] = [a] + [b].
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *             MATRIX spB       -- Pointer to Matrix B
 *  Output :   MATRIX spC       -- Pointer to Matrix C 
 *  ===================================================
 */

#ifdef __STDC__
MATRIX *matAdd( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *matAdd( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* end case not STDC */
{
MATRIX *spC;

    /* [a] : Check compatibility of matrix storage schemes, data types, and sizes */

    if(spA->Rep != spB->Rep) {
       printf("WARNING : Incompatible Matrix Representations in matAdd()\n");
       printf("WARNING : *** spA->Rep = %4d spB->Rep = %4d\n",
              spA->Rep, spB->Rep);
    }

    if(spA->Type != spB->Type) {
       printf("WARNING : Incompatible Matrix Types in matAdd()\n");
       printf("WARNING : *** spA->Type = %4d spB->Type = %4d\n",
              spA->Type, spB->Type);
    }

    /* [b] : Compute matrix operation */

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spC = matAddIndirectDouble( spA, spB );
                      break;
                 case IntegerArray:
                 case ComplexArray:
                 default:
                      fatalError("In matPrint() : Undefined spA->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
             spC = matAddSkyline( spA, spB );
             break;
        case Sparse:
        default:
             fatalError("In matAdd() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return ( spC );
}

/*
 *  ==============================================================
 *  matAddReplace() : Compute matrix assignment [A] = [A] + [B]
 *  
 *  Input :    MATRIX spA -- Pointer to Matrix A
 *             MATRIX spB -- Pointer to Matrix B
 *  Output :   MATRIX spA -- Pointer to Replacement Matrix A 
 *  ==============================================================
 */

#ifdef __STDC__
MATRIX *matAddReplace( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *matAddReplace( spA , spB )
MATRIX *spA;
MATRIX *spB;
#endif /* end case not STDC */
{

    /* [a] : Check that matrix types are compatible */

    if(spA->Rep != spB->Rep) {
       printf("WARNING : Incompatible Matrix Representations in matAddReplace()\n");
       printf("WARNING : *** m1->Rep = %4d m2->Rep = %4d\n", spA->Rep, spB->Rep);
    }

    if(spA->Type != spB->Type) {
       printf("WARNING : Incompatible Matrix Types in matAdd()\n");
       printf("WARNING : *** m1->Type = %4d m2->Type = %4d\n", spA->Type, spB->Type);
    }

    /* [b] : Compute matrix operation */

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spA = matAddReplaceIndirectDouble( spA, spB );
                      break;
                 case IntegerArray:
                 case ComplexArray:
                 default:
                      fatalError("In matAddReplace() : Undefined spA->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
        default:
             fatalError("In matAddReplace() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return ( spA );
}

/*
 *  ==========================================================
 *  matSub() : Matrix Subtraction Operation [c] = [a] - [b]
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *             MATRIX spB       -- Pointer to Matrix B
 *  Output :   MATRIX spC       -- Pointer to Matrix C 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *matSub( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *matSub( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* end case not STDC */
{
MATRIX *spC;

    /* [a] : Check that matrix types are compatible */

    if( spA->Rep != spB->Rep ) {
        printf("WARNING : Incompatible Matrix Representations in matSub()\n");
        printf("WARNING : *** m1->Rep = %4d m2->Rep = %4d\n",
                spA->Rep, spB->Rep);
    }

    if( spA->Type != spB->Type ) {
        printf("WARNING : Incompatible Matrix Types in matSub()\n");
        printf("WARNING : *** m1->Type = %4d m2->Type = %4d\n",
                      spA->Type, spB->Type);
    }

    /* [b] : Compute Matrix Operation */

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spC = matSubIndirectDouble( spA, spB );
                      break;
                 case IntegerArray:
                 case ComplexArray:
                 default:
                      fatalError("In matSub() : Undefined spA->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
             spC = matSubSkyline( spA, spB );
             break;
        case Sparse:
        default:
             fatalError("In matSub() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return ( spC );
}

/*
 *  ==============================================================
 *  matSubReplace() : Compute matrix assignment [A] = [A] - [B]
 *  
 *  Input :    MATRIX spA -- Pointer to Matrix A
 *             MATRIX spB -- Pointer to Matrix B
 *  Output :   MATRIX spA -- Pointer to Replacement Matrix A 
 *  ==============================================================
 */

#ifdef __STDC__
MATRIX *matSubReplace( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *matSubReplace( spA , spB )
MATRIX *spA;
MATRIX *spB;
#endif /* end case not STDC */
{

    /* [a] : Check that matrix types are compatible */

    if(spA->Rep != spB->Rep) {
       printf("WARNING : Incompatible Matrix Representations in matAddReplace()\n");
       printf("WARNING : *** spA->Rep = %4d spB->Rep = %4d\n",
              spA->Rep, spB->Rep);
    }

    if(spA->Type != spB->Type) {
       printf("WARNING : Incompatible Matrix Types in matAdd()\n");
       printf("WARNING : *** spA->Type = %4d spB->Type = %4d\n",
              spA->Type, spB->Type);
    }

    /* [b] : Compute matrix operation */

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spA = matSubReplaceIndirectDouble( spA, spB );
                      break;
                 case IntegerArray:
                 case ComplexArray:
                 default:
                      fatalError("In matAddReplace() : Undefined spA->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
        default:
             fatalError("In matSubReplace() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return ( spA );
}

/*
 *  ==========================================================
 *  matMult() : Matrix Multiplication [c] = [a].[b]
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *             MATRIX spB       -- Pointer to Matrix B
 *  Output :   MATRIX spC       -- Pointer to Matrix C 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *matMult( MATRIX *spA , MATRIX *spB )
#else  /* start case not STDC */
MATRIX *matMult( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* end case not STDC */
{
MATRIX *spC;

#ifdef MYDEBUG
       printf("*** Enter matMult()\n");
#endif

    /* [a] : Check that matrix data types are compatible */

    if(spA->Type != spB->Type) {
       printf("WARNING : Incompatible Matrix Types in matAdd()\n");
       printf("WARNING : *** spA->Type = %4d spB->Type = %4d\n",
              spA->Type, spB->Type);
    }

    /* [b] : Compute Matrix Multiplication for Indirect Storage */

    if((int) spA->Rep == Indirect && (int) spB->Rep == Indirect) {
       switch((int) spA->Type) {
           case DoubleArray:
                spC = matMultIndirectDouble( spA, spB );
                break;
           case IntegerArray:
           case ComplexArray:
           default:
                fatalError("In matMult() : Undefined spA->Type",
                          (char *) NULL);
                break;
       }
    }

    /* [c] : Compute Matrix Multiplication for Skyline-Indirect Storage */

    if((int) spA->Rep == Skyline && (int) spB->Rep == Indirect) {
       switch((int) spA->Type) {
           case DoubleArray:
                spC = matMultSkylineIndirectDouble( spA, spB );
                break;
           case IntegerArray:
           case ComplexArray:
           default:
                fatalError("In matMult() : Undefined spA->Type",
                          (char *) NULL);
                break;
       }
    }

    /* [d] : Compute Matrix Multiplication for Skyline-Indirect Storage */

    if((int) spA->Rep == Indirect && (int) spB->Rep == Skyline) {
       switch((int) spA->Type) {
           case DoubleArray:
                spC = matMultIndirectSkylineDouble( spA, spB );
                break;
           case IntegerArray:
           case ComplexArray:
           default:
                fatalError("In matMult() : Undefined spA->Type",
                          (char *) NULL);
                break;
       }
    }

    return ( spC );
}

/*
 *  ==============================================================
 *  matMultReplace() : Compute matrix assignment [A] = [A] * [B]
 *  
 *  Input :    MATRIX spA -- Pointer to Matrix A
 *             MATRIX spB -- Pointer to Matrix B
 *  Output :   MATRIX spA -- Pointer to Replacement Matrix A 
 *  ==============================================================
 */

#ifdef __STDC__
MATRIX *matMultReplace( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *matMultReplace( spA , spB )
MATRIX *spA;
MATRIX *spB;
#endif /* end case not STDC */
{

   /* [a] : Check that matrix representations and dimensions are compatible */

   if(spA->Rep != spB->Rep) {
      printf("WARNING : Incompatible Matrix Representations in matMultReplace()\n");
      printf("WARNING : *** spA->Rep = %4d spB->Rep = %4d\n", spA->Rep, spB->Rep);
   }

   if(spA->Type != spB->Type) {
      printf("WARNING : Incompatible Matrix Types in matMultReplace()\n");
      printf("WARNING : *** spA->Type = %4d spB->Type = %4d\n", spA->Type, spB->Type);
   }

    /* [b] : Compute [A] = [A]*[B] matrix operation */

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spA = matMultReplaceIndirectDouble( spA, spB );
                      break;
                 case IntegerArray:
                 case ComplexArray:
                 default:
                      fatalError("In matAddReplace() : Undefined spA->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
        default:
             fatalError("In matMultReplace() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return ( spA );
}


/*
 *  ==========================================================
 *  matNegate() : Compute Negative Matrix [B] = -[A];
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *  Output :   MATRIX spB       -- Pointer to Matrix B 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *matNegate( MATRIX *spA )
#else  /* start case not STDC */
MATRIX *matNegate( spA )
MATRIX *spA;
#endif /* end case not STDC */
{
MATRIX *spB;

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spB = matNegateIndirectDouble( spA );
                      break;
                 case IntegerArray:
                 case ComplexArray:
                 default:
                      fatalError("In matNegate() : Undefined spA->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
        default:
             fatalError("In matNegate() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return ( spB );
}

/*
 *  ==========================================================
 *  MatrixNegatReplace() : Matrix Replacement [A] = -[A];
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *  Output :   MATRIX spA       -- Pointer to Matrix -A 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *MatrixNegatReplace( MATRIX *spA )
#else  /* start case not STDC */
MATRIX *MatrixNegatReplace( spA )
MATRIX *spA;
#endif /* end case not STDC */
{

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spA = matNegateIndirectDouble(spA);
                      break;
                 default:
                      fatalError("In matNegate() : Undefined m->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
        default:
             fatalError("In MatrixNegatReplace() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return (spA);
}

/*
 *  ==========================================================
 *  matTranspose() : Compute Matrix Transpose [B] = [A]^T.
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *  Output :   MATRIX spB       -- Pointer to Matrix B 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *matTranspose( MATRIX *spA )
#else  /* start case not STDC */
MATRIX *matTranspose( spA )
MATRIX *spA;
#endif /* end case not STDC */
{
MATRIX *spB;

#ifdef MYDEBUG
       printf("*** In matTranspose() : Type = %3d\n", (int) spA->Type);
#endif

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spB = matTransposeIndirectDouble( spA );
                      break;
                 case IntegerArray:
                 case ComplexArray:
                 default:
                      fatalError("In matTranspose() : Undefined spA->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
        case Sparse:
        default:
             fatalError("In matTranspose() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return ( spB );
}

/*
 *  ==========================================================
 *  matCopy() : Make Copy of Matrix [B] = [A].
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *  Output :   MATRIX spB       -- Pointer to Matrix B 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *matCopy( MATRIX *spA )
#else  /* start case not STDC */
MATRIX *matCopy( spA )
MATRIX *spA;
#endif /* end case not STDC */
{
MATRIX *spB;

    switch((int) spA->Rep) {
        case Indirect:
             switch((int) spA->Type) {
                 case DoubleArray:
                      spB = matCopyIndirectDouble( spA );
                      break;
                 case IntegerArray:
                 case ComplexArray:
                 default:
                      fatalError("In matCopy() : Undefined spA->Type",
                                (char *) NULL);
                      break;
             }
             break;
        case Skyline:
             spB = matCopySkyline( spA );
             break;
        default:
             fatalError("In matPrint() : Undefined spA->Rep",
                       (char *) NULL);
             break;
    }

    return ( spB );
}


/*
 *  ==========================================================================
 *  matSolve() : Solve Systems of Linear Equations [A][X] = [B].
 * 
 *  Input :    MATRIX spA -- Pointer to (nxn) matrix A
 *             MATRIX spB -- Pointer to (nxp) r-h-s Matrix B
 *  Output :   MATRIX spX -- Pointer to (nxp) solution matrix.
 * 
 *  Also calls functions : MatrixLU() - LU decomposition.
 *                         MatrixFB() - Forward/Back substitution.
 *
 *  WARNING : This high-level equation solver works for only one decomposition
 *            followed by multiple forward/backward substitutions.
 *  ==========================================================================
 */

static VECTOR *spPivot;

#ifdef __STDC__
MATRIX *matSolve( MATRIX *spA , MATRIX *spB )
#else  /* start case not STDC */
MATRIX *matSolve( spA, spB )
MATRIX      *spA, *spB;
#endif /* end case not STDC */
{
MATRIX    *spLU;
MATRIX     *spX;

       /* [a] LUP Decomposition followed by Forward/Backward Substitution */

       switch((int) spA->Rep ) {
           case Indirect:
                switch((int) spA->Type) {
                    case DoubleArray:
                         spPivot = SetupPivotVector( spA );
                         spLU    = matLUDecompositionIndirect( spA , spPivot);
                         spX     = matLUSubstitutionIndirect( spPivot, spLU, spB);
                         break;
                    default:
                         fatalError("In MatrixLU() : Undefined A->Type", (char *) NULL);
                         break;
                }
                break;
           case Skyline:
                spLU = matLUDecompositionSkyline( spA );
                spX  = matLUBacksubstitutionSkyline( spLU, spB );
                break;
           case Sparse:
           default:
                fatalError("In matSolve() : Undefined spA->Rep",
                          (char *) NULL);
                break;
       }

       return ( spX );
}