matrix.c

有限元程序

/*
 *  ============================================================================= 
 *  ALADDIN Version 2.1.
 *                                                                     
 *  matrix.c : High-Level Functions for Matrix Operations and Printing
 *                                                                     
 *  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, Xiaoguang Chen, and Wane-Jang Lin         March 2000
 *  ============================================================================= 
 */

#include <stdio.h>

#ifdef  __STDC__
#include <stdarg.h>
#else
#include <varargs.h>
#endif

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

/*
#define DEBUG 
*/


/*
 *  ===============================================================================
 *  MatrixPrintCast(matrix, ...) : Print a matrix casting appropriate row or column
 *                                 units to those listed in [ units_m ].
 *                                                                 
 *  Usage :  In Aladdin Input File                                      
 *           MatrixPrintCast(m)                                     
 *        or MatrixPrintCast(m1, m2, units_m)                      
 *        or MatrixPrintCast(m1, m2,...,m4, units_m)               
 *                                                                 
 *        Column units are scaled with [ units_m ] = [ a , b , c ]
 *        Row units are scaled with    [ units_m ] = [ a ; b ; c ]
 *                                                                 
 *  Usage :  In c-code                                       
 *           MatrixPrint(m, (MATRIX *) NULL)                  
 *        or MatrixPrint(m1, m2, units_m, (MATRIX *) NULL)    
 *  ===============================================================================
 */

/* #define MATRIX_DEBUG */

#ifdef  __STDC__
MATRIX *MatrixPrintCast( MATRIX *first_matrix, ... ) {
#else
MATRIX *MatrixPrintCast(va_alist)
va_dcl
{
MATRIX *first_matrix;
#endif

va_list arg_ptr;
MATRIX      *p, *q;
DIMENSIONS *d1,*d2;
int col_counter, *COL_COUNT; 
int row_counter, *ROW_COUNT; 
int NO_COL_WITH_UNITS, NO_ROW_WITH_UNITS; 
int NO_COL_WITH_NO_UNITS, NO_ROW_WITH_NO_UNITS; 
int UNITS_NO_q, UNITS_NO_p, UNITS_SWITCH;
int length, length1, length2; 
int i, j, k, ii, jj, COL_ROW; 
int iFLAG  = FALSE;

#ifdef MATRIX_DEBUG
       printf("Enter MatrixPrintCast() function\n");
#endif

    /* [a] : Get units_matrix */

#ifdef  __STDC__
    va_start(arg_ptr, first_matrix);
#else
    va_start(arg_ptr);
    first_matrix = va_arg(arg_ptr, MATRIX *);
#endif

    UNITS_SWITCH = CheckUnits();
    i = 0;
    q = (MATRIX *) 1;  /* a NON NULL pointer */

    /* [b] : Compute number of matrices in argument list */

    for(p = first_matrix; q != NULL; q = va_arg(arg_ptr, MATRIX *)) {
        if( i != 0) p = q;
        i = i + 1;
    }
    va_end(arg_ptr);

    /* [c] : One matrix in argument list (call standard matrixprint function) */

    if( i == 1) {
        MatrixPrintVar( (MATRIX *) first_matrix, (MATRIX *) NULL );
        return;
    }

    /* [d] : Check that units are on for rest of function */

    if( UNITS_SWITCH == OFF) {
        FatalError("In MatrixPrintCast() : Units must be on to use this function",
                  (char *)NULL);
    }

    /* [e] : Walk along Matrix Argument List and print contents */

#ifdef __STDC__
    va_start(arg_ptr, first_matrix);
#else
    va_start(arg_ptr);
    first_matrix = va_arg(arg_ptr, MATRIX *);
#endif

    for(q = first_matrix; q != p;  q = va_arg(arg_ptr, MATRIX *)) {

        /* [e.1] : Allocate memory for row and columns ... */ 

        ROW_COUNT = (int *) MyCalloc(q->iNoRows,    sizeof(int));
        COL_COUNT = (int *) MyCalloc(q->iNoColumns, sizeof(int));

        /* [e.2] : Count number of columns in matrix q with units          */ 
        /*                                                                 */
        /*         Array COL_COUNT[i] stores the no of the i-th column in  */
        /*         matrix q having units.                                  */
        /*                                                                 */

        col_counter = 0;
        for(k = 1; k <= q->iNoColumns; k++) {
            /* Column units */
            d2 = &(q->spColUnits[k-1]);
            if(d2->length_expnt != 0 || d2->time_expnt   != 0 || d2->mass_expnt != 0 ||
               d2->temp_expnt   != 0 || d2->radian_expnt != 0 ) {
               col_counter = col_counter + 1;
               COL_COUNT[col_counter-1] = k;  /* record column with units */
            } 
        }

        /* [e.3] : Count number of rows in matrix q having units        */ 
        /*                                                              */
        /*         Array ROW_COUNT[i] stores the no of the i-th row in  */
        /*         matrix q having units.                               */
        /*                                                              */

        row_counter = 0;
        for(j = 1; j<= q->iNoRows; j++) {
            d2 = &(q->spRowUnits[j-1]);
            if(d2->length_expnt != 0 || d2->time_expnt   != 0 || d2->mass_expnt != 0 ||
               d2->temp_expnt   != 0 || d2->radian_expnt != 0 ) {
               row_counter = row_counter + 1;
               ROW_COUNT[row_counter-1] = j;  /* record row  with units */
            }
        }

        /* Save cols and rows with (and without) units.            */

        NO_ROW_WITH_UNITS = row_counter;
        NO_ROW_WITH_NO_UNITS = q->iNoRows - NO_ROW_WITH_UNITS;
        NO_COL_WITH_UNITS = col_counter;
        NO_COL_WITH_NO_UNITS = q->iNoColumns - NO_COL_WITH_UNITS;

        /* If matrix "q" is dimensionless, print error message and */
        /* terminate program execution                             */

        if(NO_COL_WITH_NO_UNITS == q->iNoColumns &&
           NO_ROW_WITH_NO_UNITS == q->iNoRows) { 
           FatalError("In MatrixPrintCase()"
                      "Trying to print out a non-dimensional matrix with units",
                       (char *) NULL);
        } 

        /* Units matrix must be larger than (1x1) matrix */

        if(p->iNoRows == 1 && p->iNoColumns == 1) {
           FatalError("In MatrixPrintCase()"
                      "Units matrix must have more than one column or row",
                       (char *) NULL);
        }

        /* p->iNoRows == 1 implies scaling of column units */

        if(p->iNoRows == 1) {
           UNITS_NO_p = p->iNoColumns;
           UNITS_NO_q = NO_COL_WITH_UNITS;
           COL_ROW    = COLUMN;
        }

        /* p->iNoColumns == 1 implies scaling of column units */

        if(p->iNoColumns == 1) {
           UNITS_NO_p = p->iNoRows;
           UNITS_NO_q = NO_ROW_WITH_UNITS;
           COL_ROW    = ROW;
        }

        /* If matrix "p" contains more rows/columns of units than matrix "q" */
        /* print error message and terminate program execution               */

        if(UNITS_NO_q < UNITS_NO_p) {
           FatalError("In MatrixPrintCase()"
                      "The units matrix contains too many rows/columns",
                      (char *) NULL); 
        }

        /* [e.5] : No of columns with units in matrix q are larger   */
        /*         than or equal to no of columns in units matrix p  */

        /* Case 1 : Scale Column Units */

        if( COL_ROW == COLUMN ) {
           for(ii = 1; ii <= UNITS_NO_q; ii++) {
               j  = COL_COUNT[ii-1];
               d1 = &(q->spColUnits[j-1]);
               for(k = 1; k <= p->iNoColumns; k++) {
                   d2 = &(p->spColUnits[k-1]);

                   if((d1->length_expnt  == d2->length_expnt) &&
                      (d1->mass_expnt    == d2->mass_expnt)   &&
                      (d1->time_expnt    == d2->time_expnt)   &&
                      (d1->temp_expnt    == d2->temp_expnt)   &&
                      (d1->radian_expnt  == d2->radian_expnt))  {
                      UnitsCopy( &(q->spColUnits[j-1]), &(p->spColUnits[k-1]) );
                      iFLAG = TRUE;
                   }
               }
            }
            if(iFLAG == FALSE)
               FatalError("In MatrixPrintCast(): Inconsistent Units",
                         (char *) NULL);
        }

        /* Case 2 : Scale Row Units */

        if( COL_ROW == ROW ) {
            for(j = 1; j <= UNITS_NO_q; j++) {
               k  = ROW_COUNT[j-1];
               d1 = &(q->spRowUnits[k-1]);
               for(jj = 1; jj <= p->iNoRows; jj++) {
                   d2 = &(p->spRowUnits[jj-1]);
                   if((d1->length_expnt  == d2->length_expnt) &&
                      (d1->mass_expnt    == d2->mass_expnt)   &&
                      (d1->time_expnt    == d2->time_expnt)   &&
                      (d1->temp_expnt    == d2->temp_expnt)   &&
                      (d1->radian_expnt  == d2->radian_expnt))  {
                         UnitsCopy( &(q->spRowUnits[k-1]), &(p->spRowUnits[jj-1]) );
                         iFLAG = TRUE;
                   }
               }
           }
           if(iFLAG == FALSE)
              FatalError("In MatrixPrintCast(): Inconsistent Units",
                        (char *) NULL );
        }

        /* Now print matrix (without simplification of units) */

        switch((int)p->eRep) {
          case  SEQUENTIAL:
          case  INDIRECT:
                switch((int) p->eType) {
                    case DOUBLE_ARRAY:
                         MatrixPrintIndirectDouble(q);
                         break;
                    case INTEGER_ARRAY:
                         MatrixPrintIndirectInteger(q);
                         break;
                    default:
                         FatalError("In MatrixPrintVar() : Undefined m->eType",
                                   (char *) NULL);
                         break;
                }
                break;
          case  SKYLINE:
                MatrixPrintSkylineDouble(q);
                break;
          default:
                FatalError("In MatrixPrintVar() : Undefined m->eRep",
                          (char *) NULL);
                break;
       } 

       free( (char *) ROW_COUNT );
       free( (char *) COL_COUNT );
   }

   va_end(arg_ptr);

#ifdef MATRIX_DEBUG
       printf("Leave MatrixPrintCast() function\n");
#endif

}


/*
 *  =====================================================================
 *  MatrixPrintVar() : Print a variable number of matrices
 *
 *  In Aladdin input file this function is called with PrintMatrix( .. );
 *  =====================================================================
 */

#ifdef  __STDC__
MATRIX *MatrixPrintVar(MATRIX *first_matrix, ...) {
#else
MATRIX *MatrixPrintVar(va_alist)
va_dcl
{
MATRIX  *first_matrix;
#endif

va_list      arg_ptr;
MATRIX            *p;
int     UNITS_SWITCH;

    /* [a] : Get first function argument */

#ifdef __STDC__
    va_start(arg_ptr, first_matrix);
#else
    va_start(arg_ptr);
    first_matrix = va_arg(arg_ptr, MATRIX *);
#endif

    /* [b] : Print matrices : units switch == ON */

    UNITS_SWITCH = CheckUnits();

    if( UNITS_SWITCH == ON ) {
    for(p = first_matrix; p != NULL;  p = va_arg(arg_ptr, MATRIX *)) {
         
        switch((int)p->eRep) {
            case  SEQUENTIAL:
            case  INDIRECT:
                  switch((int) p->eType) {
                      case DOUBLE_ARRAY:
                           MatrixUnitsSimplify(p);
                           MatrixPrintIndirectDouble(p);
                           break;
                      case INTEGER_ARRAY:
                           MatrixUnitsSimplify(p);
                           MatrixPrintIndirectInteger(p);
                           break;
                      default:
                           FatalError("In MatrixPrintVar() : Undefined m->eType",
                                     (char *) NULL);
                           break;
                  }
                  break;
            case  SKYLINE:
                  MatrixUnitsSimplify(p);
                  MatrixPrintSkylineDouble(p);
                  break;
            default:
                  FatalError("In MatrixPrintVar() : Undefined m->eRep",
                            (char *) NULL);
                  break;
        } 
    }
    }

    /* [c] : Print matrices : units switch == OFF */

    if( UNITS_SWITCH == OFF ) {
    for(p = first_matrix; p != NULL;  p = va_arg(arg_ptr, MATRIX *)) {
         
        switch((int)p->eRep) {
            case  SEQUENTIAL:
            case  INDIRECT:
                  switch((int) p->eType) {
                      case DOUBLE_ARRAY:
                           MatrixPrintIndirectDouble(p);
                           break;
                      default:
                           FatalError("In MatrixPrintVar() : Undefined matrix->eType",
                                     (char *) NULL);
                           break;
                  }
                  break;
            case  SKYLINE:
                  MatrixPrintSkylineDouble(p);
                  break;
            default:
                  FatalError("In MatrixPrintVar() : Undefined matrix->eRep",
                            (char *) NULL);
                  break;
       }
    }
    }

    va_end(arg_ptr);
}



/*
 *  ==============================================================
 *  Allocate Matrix data structure, and array[iNoRows][iNoColumns]
 *  ==============================================================