matrix_indirect.c

有限元分析源代码

/*
 *  ============================================================================= 
 *  ALADDIN Version 1.0 :
 *    matrix_indirect.c : Functions for Matrices having INDIRECT storage pattern.
 *                                                                     
 *  Copyright (C) 1995 by Mark Austin, Xiaoguang Chen, and Wane-Jang Lin
 *  Institute for Systems Research,                                           
 *  University of Maryland, College Park, MD 20742                                   
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software for any 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                                         December 1995
 *  ============================================================================= 
 */

#include <stdio.h>
#include <math.h>

#ifdef __STDC__
#include <stddef.h>
#include <stdlib.h>
#include <stdarg.h>
#endif

#include "units.h"
#include "matrix.h"
#include "vector.h"
#include "defs.h"
#include "miscellaneous.h"

/* #define DEBUG */
/* #define LU_DEBUG */


/*
 *  =======================================================================
 *  MatrixAllocIndirect() : Allocate memory for Matrix data structure with
 *  INDIRECT storage pattern
 *  
 *  Input  :  char *cpMatrixName -- Pointer to name of matrix.
 *         :  DATA_TYPE eType    -- Data type to be stored in Matrix.
 *         :  int  iNoRows       -- No of Rows in Matrix.
 *         :  int  iNoColumns    -- No of Columns in Matrix.
 *  Output :  MATRIX *Matrix     -- Pointer to matrix data structure.
 *  =======================================================================
 */

#ifdef __STDC__
MATRIX *MatrixAllocIndirect( char *cpMatrixName, DATA_TYPE  eType,
                             int  iNoRows, int   iNoColumns)
#else  /* Start case not STDC */
MATRIX *MatrixAllocIndirect( cpMatrixName, eType, iNoRows, iNoColumns)
char     *cpMatrixName;
DATA_TYPE        eType;
int            iNoRows;
int         iNoColumns;
#endif /* End case not STDC */
{
MATRIX *spMatrix;

   /* [a] : Allocate matrix parent data structure, and cpMatrixName */

      spMatrix = (MATRIX *) MyMalloc( sizeof(MATRIX) );
      if(cpMatrixName != (char *)NULL)
         spMatrix->cpMatrixName = SaveString(cpMatrixName);
      else 
         spMatrix->cpMatrixName = (char *) NULL;

   /* [b] : Set parameters and allocate memory for matrix uMatrix */

      spMatrix->eRep       = INDIRECT;
      spMatrix->eType      = eType;
      spMatrix->iNoRows    = iNoRows;
      spMatrix->iNoColumns = iNoColumns;
      if( CheckUnits() == ON ) {
         spMatrix->spRowUnits 
         = (DIMENSIONS *) MyCalloc(iNoRows, sizeof(DIMENSIONS));
         spMatrix->spColUnits
         = (DIMENSIONS *) MyCalloc(iNoColumns, sizeof(DIMENSIONS));
      }
      else {
         spMatrix->spRowUnits = (DIMENSIONS *)NULL;
         spMatrix->spColUnits = (DIMENSIONS *)NULL;
      }

      switch((int) spMatrix->eType) {
          case DOUBLE_ARRAY:
               spMatrix->uMatrix.daa 
               = MatrixAllocIndirectDouble( iNoRows, iNoColumns);
               break;
          case INTEGER_ARRAY:
               spMatrix->uMatrix.iaa 
               = MatrixAllocIndirectInteger( iNoRows, iNoColumns);
               break;
          case COMPLEX_ARRAY:
               FatalError("In MatrixAllocIndirect() : spMatrix->eType not implemented",
                         (char *) NULL);
               break;
          default:
               FatalError("In MatrixAllocIndirect() : Undefined spMatrix->eType",
                         (char *) NULL);
               break;
      }

      return (spMatrix);
}

/*
 *  =======================================================================
 *  MatrixAllocIndirectDouble() : Allocate INDIRECT storage pattern for
 *  DOUBLE data types
 *  
 *  Input  :  int  iNoRows       -- No of Rows in Matrix.
 *         :  int  iNoColumns    -- No of Columns in Matrix.
 *  Output :  double **Matrix    -- Pointer to matrix data structure.
 *  =======================================================================
 */

#ifdef __STDC__
double **MatrixAllocIndirectDouble( int iNoRows, int iNoColumns)
#else  /* Start case not STDC */
double **MatrixAllocIndirectDouble( iNoRows, iNoColumns)
int    iNoRows;
int iNoColumns;
#endif /* End case not STDC */
{
double **Matrix;
int ii;

      Matrix = (double **) MyCalloc( iNoRows, sizeof(double *));
      for(ii = 1; ii <= iNoRows; ii++) 
          Matrix[ii-1] = (double *) MyCalloc( iNoColumns, sizeof(double));

      return (Matrix);
}

/*
 *  =======================================================================
 *  MatrixAllocIndirectInteger() : Allocate INDIRECT storage pattern for
 *                                 INTEGER data types
 *  
 *  Input  :  int  iNoRows       -- No of Rows in Matrix.
 *         :  int  iNoColumns    -- No of Columns in Matrix.
 *  Output :  double **Matrix    -- Pointer to matrix data structure.
 *  =======================================================================
 */

#ifdef __STDC__
int **MatrixAllocIndirectInteger( int iNoRows, int iNoColumns)
#else  /* Start case not STDC */
int **MatrixAllocIndirectInteger( iNoRows, iNoColumns)
int    iNoRows;
int iNoColumns;
#endif /* End case not STDC */
{
int **Matrix;
int ii;

      Matrix = (int **) MyCalloc( iNoRows, sizeof(int *));
      for(ii = 1; ii <= iNoRows; ii++) 
          Matrix[ii-1] = (int *) MyCalloc( iNoColumns, sizeof(int));

      return (Matrix);
}


/*
 *  =======================================================================
 *  MatrixPrintIndirectDouble() : Print a Matrix [iNoRows][iNoColumns] of
 *                                data type DOUBLE 
 *  
 *  MatrixPrintIndirectInteger() : Print a Matrix [iNoRows][iNoColumns] of
 *                                 data type INTEGER
 *  
 *  Where --  COLUMNS_ACROSS_PAGE  = Number of columns printed across page.
 *            spA->iNoRows    = Number of rows in matrix.
 *            spA->iNoColumns = Number of columns in matrix.
 *            iFirstColumn         = Number of first column in block
 *            iLastColumn          = Number of last  column in block
 *            ib                   = Current No of Matrix Block.
 * 
 *  Input  :  MATRIX *spA         -- Pointer to matrix data structure.
 *  Output :  void
 *  =======================================================================
 */

enum { COLUMNS_ACROSS_PAGE = 6 };                                        /* Item [a] */

#ifdef __STDC__                                                  
void MatrixPrintIndirectDouble( MATRIX *spA )
#else  /* Start case not STDC */
void MatrixPrintIndirectDouble( spA )
MATRIX  *spA;
#endif /* End case not STDC */
{
int ii, ij, ib;    
int iNoBlocks; 
int iFirstColumn, iLastColumn; 
double da;
int UNITS_SWITCH;

     UNITS_SWITCH = CheckUnits();
    /* [a] : Compute no of blocks of rows to be printed */               /* Item [c] */

    if( spA->iNoColumns % ((int) COLUMNS_ACROSS_PAGE) == 0 )                          
        iNoBlocks = (spA->iNoColumns/((int) COLUMNS_ACROSS_PAGE));
    else
        iNoBlocks = (spA->iNoColumns/((int) COLUMNS_ACROSS_PAGE)) + 1;

    /* [b] : Loop over blocks of rows */

    for( ib = 1; ib <= iNoBlocks; ib++ ) {                               /* Item [d] */

         iFirstColumn = (ib-1)*((int) COLUMNS_ACROSS_PAGE) + 1;
         iLastColumn  = (int) MIN( ib*((int) COLUMNS_ACROSS_PAGE) , spA->iNoColumns );

         /* [c] : Print title of matrix at top of each block */          /* Item [e] */
       
         if( spA->cpMatrixName != NULL ) 
             printf("\nMATRIX : \"%s\"\n\n", spA->cpMatrixName);
         else 
             printf("\nMATRIX : \"UNTITLED\"\n\n");

         /* [d] : Label row and column nos */

         printf ("row/col         ");
         for( ii = iFirstColumn; ii <= iLastColumn; ii++ )
             printf("       %3d   ", ii);
         printf("\n");

         switch( UNITS_SWITCH) {
           case ON:
               printf ("        units   " );
               for( ii = iFirstColumn; ii <= iLastColumn; ii++ )
                  if(spA->spColUnits[ii-1].units_name != NULL) {
                     printf("%10s   ",spA->spColUnits[ii-1].units_name);
                  }
                  else
                    printf("             ");
               printf("\n");

             /* [e] : Print Contents of Matrix */   /* Item [f] */

               for( ii = 1; ii <= spA->iNoRows; ii++ ) {
                 printf(" %3d ", ii);
                 if(spA->spRowUnits[ii-1].units_name != NULL) {
                    printf("%8s ",spA->spRowUnits[ii-1].units_name);
                 }
                 else
                    printf("         ");
                 for(ij  = iFirstColumn; ij <= iLastColumn; ij++) {
                     da = MatrixContentScaleIndirectDouble(spA, ii, ij);
                     printf(" %12.5e", da );
                 }
                 printf("\n");
               }
               break;
           case OFF:

               /* [e] : Print Contents of Matrix */                            /* Item [f] */

               for( ii = 1; ii <= spA->iNoRows; ii++ ) {
                 printf(" %3d ", ii);
                 printf("         ");
                 for( ij  = iFirstColumn; ij <= iLastColumn; ij++)
                     printf(" %12.5e", spA->uMatrix.daa[ ii-1 ][ ij-1 ]);
                 printf("\n");
               }
               break;
           default:
               break;
         }
    }
}

#ifdef __STDC__                                                  
void MatrixPrintIndirectInteger( MATRIX *spA )
#else  /* Start case not STDC */
void MatrixPrintIndirectInteger( spA )
MATRIX  *spA;
#endif /* End case not STDC */
{
int ii, ij, ib;    
int iNoBlocks; 
int iFirstColumn, iLastColumn; 
int UNITS_SWITCH;

     UNITS_SWITCH = CheckUnits();
    /* [a] : Compute no of blocks of rows to be printed */               /* Item [c] */

    if( spA->iNoColumns % ((int) COLUMNS_ACROSS_PAGE) == 0 )                          
        iNoBlocks = (spA->iNoColumns/((int) COLUMNS_ACROSS_PAGE));
    else
        iNoBlocks = (spA->iNoColumns/((int) COLUMNS_ACROSS_PAGE)) + 1;

    /* [b] : Loop over blocks of rows */

    for( ib = 1; ib <= iNoBlocks; ib++ ) {                               /* Item [d] */

         iFirstColumn = (ib-1)*((int) COLUMNS_ACROSS_PAGE) + 1;
         iLastColumn  = (int) MIN( ib*((int) COLUMNS_ACROSS_PAGE) , spA->iNoColumns );

         /* [c] : Print title of matrix at top of each block */          /* Item [e] */
       
         if( spA->cpMatrixName != NULL ) 
             (void) printf("\nMATRIX : \"%s\"\n\n", spA->cpMatrixName);
         else 
             (void) printf("\nMATRIX : \"UNTITLED\"\n\n");

         /* [d] : Label row and column nos */

         (void) printf ("row/col        ");
         for( ii = iFirstColumn; ii <= iLastColumn; ii++ )
              (void) printf("%3d          ", ii);
         (void) printf("\n");

         switch( UNITS_SWITCH) {
           case ON:
               printf ("      units   " );
               for( ii = iFirstColumn; ii <= iLastColumn; ii++ )
                  if(spA->spColUnits[ii-1].units_name != NULL) {
                     printf("%10s   ",spA->spColUnits[ii-1].units_name);
                  }
               printf("\n");

               /* [e] : Print Contents of Matrix */                            /* Item [f] */

               for( ii = 1; ii <= spA->iNoRows; ii++ ) {
                    (void) printf(" %3d ", ii);
                    if(spA->spRowUnits[ii-1].units_name != NULL) {
                       printf("%8s ",spA->spRowUnits[ii-1].units_name);
                    }
                    else
                       printf("         ");
                    for( ij  = iFirstColumn; ij <= iLastColumn; ij++)
                         (void) printf(" %12d", spA->uMatrix.iaa[ ii-1 ][ ij-1 ]);
                    (void) printf("\n");
               }
               break;
           case OFF:
               /* [e] : Print Contents of Matrix */                            /* Item [f] */

               for( ii = 1; ii <= spA->iNoRows; ii++ ) {
                    printf(" %3d ", ii);
                    printf("         ");
                    for( ij  = iFirstColumn; ij <= iLastColumn; ij++)
                         (void) printf(" %12d", spA->uMatrix.iaa[ ii-1 ][ ij-1 ]);
                    printf("\n");
               }
               break;
           default:
               break;
         }
    }
}


/*
 *  =======================================================================
 *  MatrixFreeIndirect() : Free memory for INDIRECT storage.
 *  
 *  Input  :  MATRIX *Matrix  -- Pointer to matrix data structure.
 *  Output :  void
 *  =======================================================================
 */

#ifdef __STDC__
void MatrixFreeIndirect( MATRIX *spA )
#else  /* Start case not STDC */
void MatrixFreeIndirect( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
int i;

   if( spA == (MATRIX *)NULL )   return;

   /* [a] : Free body of matrix[iNoRows][iNoColumns] */

      switch(spA->eType) {
          case DOUBLE_ARRAY:
               MatrixFreeIndirectDouble( spA->uMatrix.daa, spA->iNoRows);
               break;