matrix.c

利用语言编写的有限元分析软件

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

#include <stdio.h>

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

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

/*
#define DEBUG 
*/


/*
 *  ================================================================
 *  MatrixPrint(matrix, ...)                            
 *                                                                 
 *  Usage :  In input file                                      
 *           MatrixPrint(m)                                     
 *        or MatrixPrint(m1, m2, units_m)                      
 *        or MatrixPrint(m1, m2,...,m4, units_m)               
 *  Usage :  In c-code                                       
 *           MatrixPrint(m, (MATRIX *) NULL)                  
 *        or MatrixPrint(m1, m2, units_m, (MATRIX *) NULL)    
 *  ================================================================
 */

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

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

    /* [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

    /* [b] : Find out units_matrix */

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

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

    if( i == 1) {
      switch( p->eRep ) {
        case INDIRECT :
           switch((int) p->eType) {
              case DOUBLE_ARRAY:
                if( UNITS_SWITCH == ON )
                   MatrixUnitsSimplify(p);
                MatrixPrintIndirectDouble(p);
                break;
              case INTEGER_ARRAY:
                if( UNITS_SWITCH == ON )
                   MatrixUnitsSimplify(p);
                MatrixPrintIndirectInteger(p);
                break;
              default:
                   FatalError("In MatrixPrint() : Undefined m->eType",(char *)NULL);
                break;
           }  /* end of case INDIRECT */
           break;
        case SKYLINE :
           switch((int) p->eType) {
              case DOUBLE_ARRAY:
                if( UNITS_SWITCH == ON )
                   MatrixUnitsSimplify(p);
                MatrixPrintSkylineDouble(p);
                break;
              default:
                   FatalError("In MatrixPrint() : Undefined m->eType",(char *)NULL);
                break;
           }  /* end of case SKYLINE */
           break;
        default :
           FatalError("In MatrixPrint() : Undefined m->eRep",(char *)NULL);
           break;
      }  /* end of switch of p->eRep */ 
    }  /* end of if( i==1 ) */
    va_end(arg_ptr);

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

      if(i > 1 && UNITS_SWITCH == ON ) { 
         for(q = first_matrix; q != p;  q = va_arg(arg_ptr, MATRIX *)) {
             ROW_COUNT = (int *) MyCalloc(q->iNoRows, sizeof(int));
             COL_COUNT = (int *) MyCalloc(q->iNoColumns, sizeof(int));
             Max_Row_Counter = 0;
             Max_Col_Counter = 0;

      /*********************************************************************/
      /*   Count the number of columns with units in matrix q              */ 
      /*********************************************************************/
          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) {
                 col_counter = col_counter + 1;
                 COL_COUNT[col_counter-1] = k;  /* record column with units */
                 YESNO_COL = YES;
              }
              else 
                 YESNO_COL = NO;
          }
          row_counter = 0;
          for(j = 1; j<= q->iNoRows; j++) {
              /* Row units */
              d2 = &(q->spRowUnits[j-1]);
              if(d2->length_expnt != 0 || d2->time_expnt   != 0 ||
                 d2->mass_expnt   != 0 || d2->temp_expnt   != 0) {
                 row_counter = row_counter + 1;
                 ROW_COUNT[row_counter-1] = j;  /* record row  with units */
                 YESNO_ROW = YES;
              }
              else
                 YESNO_ROW = NO;
          }
          if(row_counter > Max_Row_Counter) {
             Max_Row_Counter = row_counter;
             MAX_ROW_COUNT   = ROW_COUNT;
          }

          if(col_counter > Max_Col_Counter) {
             Max_Col_Counter = col_counter;
             MAX_COL_COUNT   = COL_COUNT;
          }

          NO_ROW_WITH_UNITS = Max_Row_Counter;
          NO_ROW_WITH_NO_UNITS = q->iNoRows - NO_ROW_WITH_UNITS;

          NO_COL_WITH_UNITS = Max_Col_Counter;
          NO_COL_WITH_NO_UNITS = q->iNoColumns - NO_COL_WITH_UNITS;

      /*********************************************************************/
         if(NO_COL_WITH_NO_UNITS == q->iNoColumns &&
            NO_ROW_WITH_NO_UNITS == q->iNoRows) {     /*  matrix q have no units  */
             printf(" Inconsistent Units \n");
             FatalError("in MatrixPrint(): trying to print out a non-dimensional matrix with units",(char *)NULL);
         }
         else {
            if(p->iNoRows == 1) {
               UNITS_NO_p = p->iNoColumns;
               UNITS_NO_q = NO_COL_WITH_UNITS;
               COL_ROW    = COLUMN;
            }
            if(p->iNoColumns == 1) {
               UNITS_NO_p = p->iNoRows;
               UNITS_NO_q = NO_ROW_WITH_UNITS;
               COL_ROW    = ROW;
            }
            if(UNITS_NO_q < UNITS_NO_p) {
               printf(" Number of units does not match:\n"); 
               printf("   : Too many units in units matrix \n"); 
               printf("   : No of columns with units in matrix of []%dx%d is %d \n",
                            q->iNoRows, q->iNoColumns, NO_COL_WITH_UNITS); 
               printf("   : No of columns of units matrix is %d \n", p->iNoColumns); 
               FatalError(" Fatal Error in MatrixPrint(): ",(char *)NULL); 
            }
            else { /* No of columns with units in matrix q are larger   */
                   /* than or equal to no of columns in units matrix p  */
               switch(COL_ROW)  {
                   case COLUMN :
                    for(ii = 1; ii <= UNITS_NO_q; ii++) {
                        j = MAX_COL_COUNT[ii-1];
                        d1 = &(q->spColUnits[j-1]);
                        for(k = 1; k <= p->iNoColumns; k++){
                            d2 = &(p->spColUnits[k-1]);
			    if(SameUnits(d1,d2) == TRUE) {
                               UnitsCopy( &(q->spColUnits[j-1]), &(p->spColUnits[k-1]) );
                               iFLAG1 = TRUE;
                            } else 
                               iFLAG1 = MAX(iFLAG1, FALSE);
                            if(iFLAG1 == TRUE) break;
                        }
                        if(iFLAG1 == TRUE) iFLAG = TRUE;
                        iFLAG1 = FALSE;
                    }
                    if(iFLAG == FALSE)
                       FatalError(" Fatal Error in MatrixPrint(): inconsistent units ",(char *)NULL);
                   break;
                   case ROW:
                    for(j = 1; j <= UNITS_NO_q; j++) {
                        k = MAX_ROW_COUNT[j-1];
                        d1 = &(q->spRowUnits[k-1]);
                        for(jj = 1; jj <= p->iNoRows; jj++){
                            d2 = &(p->spRowUnits[jj-1]);
			    if(SameUnits(d1,d2) == TRUE) {
                               UnitsCopy( &(q->spRowUnits[k-1]), &(p->spRowUnits[jj-1]) );
                               iFLAG1 = TRUE;
                            } else
                               iFLAG1 = MAX(iFLAG1, FALSE);
                            if(iFLAG1 == TRUE) break;
                        }
                        if(iFLAG1 == TRUE) iFLAG = TRUE;
                        iFLAG1 = FALSE;
                    }
                    if(iFLAG == FALSE)
                       FatalError(" MatrixPrint(): inconsistent units ",(char *)NULL);
                   break;
                   default :
                   break;
                 }
            }
          }
          switch( q->eRep ) {
            case INDIRECT :
              switch((int) q->eType) {
                 case DOUBLE_ARRAY:
                   MatrixUnitsSimplify(q);
                   MatrixPrintIndirectDouble(q);
                   break;
                 case INTEGER_ARRAY:
                   MatrixUnitsSimplify(q);
                   MatrixPrintIndirectInteger(q);
                   break;
                 default:
                   FatalError("In MatrixPrint() : Undefined m->eType",(char *)NULL);
                   break;
              }  /* end of case INDIRECT */
              break;
            case SKYLINE :
              switch((int) q->eType) {
                 case DOUBLE_ARRAY:
                   MatrixUnitsSimplify(q);
                   MatrixPrintSkylineDouble(q);
                   break;
                 default:
                   FatalError("In MatrixPrint() : Undefined m->eType",(char *)NULL);
                   break;
             }  /* end of case SKYLINE */
             break;
           default:
             FatalError("In MatrixPrint() : Undefined m->eRep",(char *)NULL);
             break;
        }  /* end of switch( q->eRep ) */
        free( (char *) ROW_COUNT );
        free( (char *) COL_COUNT );
      }  /* end of for loop */
      va_end(arg_ptr);
    }  /* end of if( i>1 && UNITS_SWITCH == ON ) */
    if(i > 1 &&  UNITS_SWITCH == OFF) {
      printf("*** You need to set units on to use this function \n");
      FatalError("*** Units are off",(char *)NULL);
    }

}

/*
 *  ===================================================================
 *  MatrixPrintVar() : Print a variable number of matrices
 * 
 *  Input  :
 *  Output :
 *  ===================================================================
 */

#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;