matrix_skyline.c

有限元程序

/*
 *  ============================================================================= 
 *  ALADDIN Version 2.1.
 *                                                                     
 *  matrix_skyline.c : Functions for Operations on Skyline Matrices.
 *                                                                     
 *  Copyright (C) 1995-2000 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 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 software may not be sold or included in commercial software
 *     products without a license. 
 *  5. This notice is to remain intact.
 *                                                                    
 *  Written by: Mark Austin and Lanheng Jin                              May 1994
 *  ============================================================================= 
 */

#include <stdio.h>
#include <stdlib.h>

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

/* #define MYDEBUG */

/*
 *  =======================================================================
 *  MatrixAllocSkyline() : Allocate (and free) memory for SKYLINE matrix
 *                         data structure. 
 *  
 *  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.
 *         :  int  *ld           -- Length of each column in skyline profile
 *  Output :  MATRIX *spA        -- Pointer to matrix data structure.
 *  =======================================================================
 */

#ifdef __STDC__
MATRIX *MatrixAllocSkyline( char *cpMatrixName,  DATA_TYPE eType,
                            int  iNoRows,  int  iNoColumns,  int  *ld)
#else  /* start case not STDC */
MATRIX *MatrixAllocSkyline( cpMatrixName,  eType,  iNoRows, iNoColumns, ld)
char     *cpMatrixName;
DATA_TYPE        eType;
int            iNoRows;
int         iNoColumns;
int                *ld;
#endif /* end case not STDC */
{
MATRIX *spA;
int      ii;

   /* [a] : Check input parameters/arrays */

      if(ld == (int *)NULL) 
         FatalError("In MatrixAllocSkyline() : ld[] = NULL",
                   (char *) NULL);

   /* [b] : Allocate matrix parent data structure, and name */

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

   /* [c] : Set parameters and allocate memory for matrix array */

      spA->eRep       = SKYLINE;
      spA->eType      = eType;
      spA->iNoRows    = iNoRows;
      spA->iNoColumns = iNoColumns;

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

      switch((int) spA->eType) {
          case DOUBLE_ARRAY:
               spA->uMatrix.daa = (double **)
                                  MyCalloc( iNoColumns, sizeof(double*));
               for(ii = 0; ii < iNoColumns; ii = ii + 1) {
                   spA->uMatrix.daa[ ii ] = (double *)
                                            MyCalloc((ld[ii]+1), sizeof(double));
                   spA->uMatrix.daa[ ii ][0] = ld[ ii ];
               }
               break;
          default:
               FatalError("In MatrixAllocSkyline() : Undefined spA->eType",
                         (char *) NULL);
               break;
      }

      return ( spA );
}

/*
 *  =======================================================================
 *  MatrixFreeSkyline() : Free memory in Skyline Matrix
 *  
 *  Input  :  MATRIX *Matrix     -- Pointer to matrix data structure.
 *  Output :  void
 *  =======================================================================
 */

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

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

   /* [a] : Free body of skyline matrix[][] */

      switch( spA->eType ) {
          case DOUBLE_ARRAY:
               for(ii = 1; ii <= spA->iNoColumns; ii++)
                   free ((char *) (spA->uMatrix.daa[ ii-1 ]));
               free ((char *) (spA->uMatrix.daa));
               break;
          default:
               FatalError("In MatrixFreeSkyline() : Undefined spA->eType",
                         (char *) NULL);
               break;
      }

   /* [b] : Free name and parent data structure for matrix */

      if( CheckUnits() == ON ) {
         for( ii=1 ; ii<=spA->iNoRows ; ii++ )
             free ((char *) spA->spRowUnits[ii-1].units_name );
         for( ii=1 ; ii<=spA->iNoColumns ; ii++ )
             free ((char *) spA->spColUnits[ii-1].units_name );
         free ((char *) spA->spRowUnits);
         free ((char *) spA->spColUnits);
      }
      free ((char *) spA->cpMatrixName);
      free ((char *) spA);
      spA = (MATRIX *)NULL;
}

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

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

#ifdef __STDC__                                                  
void MatrixPrintSkylineDouble( MATRIX *spA )
#else  /* Start case not STDC */
void MatrixPrintSkylineDouble( spA )
MATRIX  *spA;
#endif /* End case not STDC */
{
int ii, ij, ik, im, 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("\nSKYLINE MATRIX : \"%s\"\n\n", spA->cpMatrixName);
         else 
             printf("\nSKYLINE MATRIX : \"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 = MatrixContentScaleSkyline(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++) {

                        ik = (int) MIN( ii, ij );
                        im = (int) MAX( ii, ij );

                        if((im-ik+1) <= spA->uMatrix.daa[im-1][0])
                             printf(" %12.5e ", spA->uMatrix.daa[ im-1 ][ im-ik+1 ]);
                        else
                             printf(" %12.5e ", 0.0);

                   }
                   printf("\n");
              }
              break;
           default:
              break;
         }
    }
}


/*
 *  ================================================================
 *  Matrix Operations : Addition, Subtraction, Multiplication, Copy.
 *  ================================================================
 */

/*
 *  ========================================================================
 *  MatrixReallocSkyline() : Reallocate Memory for Skyline Matrix. Eliminate
 *                           matrix elements that are smaller than MINREAL
 * 
 *  Input  :  MATRIX *spA  -- Pointer to matrix data structure.
 *  Output :  MATRIX *spB  -- Pointer to new skyline matrix.
 *  ========================================================================
 */

#ifdef __STDC__
MATRIX *MatrixReallocSkyline ( MATRIX *spA )
#else  /* Start case not STDC */
MATRIX *MatrixReallocSkyline ( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
MATRIX     *spB;
int ii, ij, *ld;
int      length;

    if( spA == NULL )
        FatalError("In MatrixReallocSkyline() : spA = NULL",
                  (char *) NULL);

    ld = iVectorAlloc( spA->iNoColumns ); 
    for(ii=0; ii < spA->iNoColumns; ii = ii + 1) {
        for (ij = spA->uMatrix.daa[ii][0]; (ij > 1) &&
             (ABS(spA->uMatrix.daa[ii][ij]) < MINREAL); ij = ij - 1);
        ld[ii] = ij;
    }

    spB = MatrixAllocSkyline( spA->cpMatrixName, DOUBLE_ARRAY,
                              spA->iNoRows, spA->iNoColumns, ld);

    for(ii = 0; ii < spA->iNoColumns; ii++)
        for(ij = 1; ij <= spB->uMatrix.daa[ii][0]; ij++)
            spB->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij];

    if( CheckUnits() == ON ) {
       for(ii = 1; ii <= spA->iNoRows; ii++)
           UnitsCopy( &(spB->spRowUnits[ii-1]),  &(spA->spRowUnits[ii-1]) );
       for(ij = 1; ij <= spA->iNoColumns; ij++)
           UnitsCopy( &(spB->spColUnits[ij-1]),  &(spA->spColUnits[ij-1]) );
    }

    MatrixFreeSkyline(spA);
    free((char *) ld);
    return (spB); 
}

/*
 *  =======================================================================
 *  MatrixAddSkyline() : Add Skyline Matrices
 * 
 *  Input  :  MATRIX *spA        -- Pointer to matrix data structure.
 *         :  MATRIX *spB        -- Pointer to matrix data structure.
 *  Output :  MATRIX *spC        -- Pointer to matrix sum
 *  =======================================================================
 */

#ifdef __STDC__
MATRIX *MatrixAddSkyline( MATRIX *spA , MATRIX *spB )
#else  /* Start case not STDC */
MATRIX *MatrixAddSkyline( spA, spB )
MATRIX *spA, *spB;
#endif /* End case not STDC */
{
MATRIX *spC;
int *ld, min_ldi;
int  ii, ij;
int  length, length1, length2;
DIMENSIONS *d1, *d2;
int UNITS_SWITCH;