matrix.c

有限元程序

 */

#ifdef __STDC__
MATRIX *MatrixAllocate(MATRIX *m)
#else
MATRIX *MatrixAllocate(m)
MATRIX *m;
#endif
{
MATRIX         *m1;
int        iNoRows;
int     iNoColumns;
int              i;

   /* [a] : Check size and values of argument matrix m */

   if((m->iNoRows < 1) || (m->iNoColumns < 2))
       FatalError("In MatrixAllocate() : Argument Matrix Too Small",(char *)NULL);

   /* [b] : Allocate and return Matrix */

   iNoRows    = (int) m->uMatrix.daa[0][0];
   iNoColumns = (int) m->uMatrix.daa[0][1];
 
   m1 = MatrixAllocIndirect((char *) NULL, DOUBLE_ARRAY, iNoRows, iNoColumns);

   /* [c] : Zero out units in row/column units vectors */

   if( CheckUnits() == ON ) {
      for( i=1 ; i<= iNoRows ; i++ )
        ZeroUnits(&(m1->spRowUnits[i-1]));
      for( i=1 ; i<= iNoColumns ; i++ )
        ZeroUnits(&(m1->spColUnits[i-1]));
   }

   return(m1);
}

#ifdef __STDC__
MATRIX *MatrixDiag(MATRIX *m)
#else
MATRIX *MatrixDiag(m)
MATRIX *m;
#endif
{
MATRIX              *m1;
int    matrix_dimension;
int           length, i;

   /* [a] : Check size and values of argument matrix m */

   if((m->iNoRows < 1) || (m->iNoColumns < 2))
       FatalError("In MatrixDiag() : Argument Matrix Too Small",(char *)NULL);

   /* [b] : Check dimensions of matrix */

   matrix_dimension = (int) m->uMatrix.daa[0][0];
   if(matrix_dimension < 1)
      FatalError("In MatrixDiag() : Matrix dimension < 1",(char *)NULL);

   /* [c] : Instantiate and return Matrix                    */
   /*                                                        */
   /* Note : UNITS OF MATRIX ELEMENTS ARE STORED             */
   /*        IN COLUMN_UNITS_BUF AS AN DEFAULT               */

   m1 = MatrixAllocIndirect((char *) NULL, DOUBLE_ARRAY, matrix_dimension, matrix_dimension);

   for(i = 1; i <= matrix_dimension; i++) {
       m1->uMatrix.daa[i-1][i-1]= (double) m->uMatrix.daa[0][1];

       if( CheckUnits() == ON ) {
           UnitsCopy(&(m1->spColUnits[i-1]), &(m->spColUnits[1]));
           ZeroUnits(&(m1->spRowUnits[i-1]));
       }
   }

   return(m1);
}

/*==================================*/
/* Usage :                          */
/* MatrixZero([n]):                 */
/*   generate a nxn matrix of zeros */
/* MatrixZero([n,m]):               */
/*   generate a nxm matrix of zeros */
/*                                  */
/*==================================*/

#ifdef __STDC__
MATRIX *MatrixZero(MATRIX *m)
#else
MATRIX *MatrixZero(m)
MATRIX *m;
#endif
{
MATRIX         *m1;
int        iNoRows;
int     iNoColumns;
int              i;
int              j;


#ifdef DEBUG
      printf("Enter MatrixZero() : \n");
#endif 
   /* Check size and values of argument matrix m */

      if((m->iNoRows < 1) || (m->iNoColumns < 1))
         FatalError("In MatrixZero() : Argument Matrix Too Small",(char *)NULL);

      iNoRows    = (int) m->uMatrix.daa[0][0];

      switch((int) m->iNoColumns) {
        case 1:
          iNoColumns = (int) m->uMatrix.daa[0][0];
        break;
        case 2:
          iNoColumns = (int) m->uMatrix.daa[0][1];
        break;
        default:
          printf("*** too many columns in the matrix. \n Check Zero() in input file \n");
          FatalError("*** Syntax error: two columns at most in Zero([]), (char *) NULL");
        break;
      }

      m1 = MatrixAllocIndirect((char *) NULL, DOUBLE_ARRAY, iNoRows, iNoColumns);
      if(CheckUnits()==ON) {
         for(i = 1; i <= iNoRows; i++)
             ZeroUnits(&(m1->spRowUnits[i-1]));
         for(i = 1; i <= iNoColumns; i++)
             ZeroUnits(&(m1->spColUnits[i-1]));
      }

#ifdef DEBUG
      printf("In MatrixZero():  iNoRows    = %d \n", iNoRows);
      printf("In MatrixZero():  iNoColumns = %d \n", iNoColumns);
#endif

      for(i = 1; i <= iNoRows; i++)
         for(j = 1; j <= iNoColumns; j++)
             m1->uMatrix.daa[i-1][j-1] = 0.0;

#ifdef DEBUG
      printf("Leaving MatrixZero() :  iNoRows = %d \n", iNoRows);
#endif 
      return(m1);
}

#ifdef __STDC__
MATRIX *MatrixOne(MATRIX *m)
#else
MATRIX *MatrixOne(m)
MATRIX *m;
#endif
{
MATRIX         *m1;
int        iNoRows;
int     iNoColumns;
int            i,j;

   /* Check size and values of argument matrix m */

      if((m->iNoRows < 1) || (m->iNoColumns < 1))
         FatalError("In MatrixOne() : Argument Matrix Too Small",(char *)NULL);

      iNoRows    = (int) m->uMatrix.daa[0][0];

      switch((int) m->iNoColumns) {
        case 1:
          iNoColumns = (int) m->uMatrix.daa[0][0];
        break;
        case 2:
          iNoColumns = (int) m->uMatrix.daa[0][1];
        break;
        default:
          printf("*** too many columns in the matrix. \n Check Zero() in input file \n");
          FatalError("*** Syntax error: two columns at most in Zero([]), (char *) NULL");
        break;
      }

      m1 = MatrixAllocIndirect((char *) NULL, DOUBLE_ARRAY, iNoRows, iNoColumns);
      if(CheckUnits()==ON) {
         for(i = 1; i <= iNoRows; i++)
             ZeroUnits(&(m1->spRowUnits[i-1]));
         for(i = 1; i <= iNoColumns; i++)
             ZeroUnits(&(m1->spColUnits[i-1]));
      }

      for(i = 1; i <= iNoRows; i++)
         for(j = 1; j <= iNoColumns; j++)
             m1->uMatrix.daa[i-1][j-1] = 1.0;

      return(m1);
}


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

#ifdef __STDC__
void MatrixFree( MATRIX *spA )
#else  /* start case not STDC */
void MatrixFree( spA )
MATRIX *spA;
#endif /* end case not STDC */
{
      if ( spA==(MATRIX *)NULL )   return;

      switch((int) spA->eRep) {
          case INDIRECT:
               MatrixFreeIndirect( spA );
               break;
          case SKYLINE:
               MatrixFreeSkyline( spA );
               break;
          default:
               FatalError("In MatrixFree() : Undefined spA->eRep",
                         (char *) NULL);
               break;
      }
}


/*
 *  ===================================================
 *  MatrixAdd() : 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 *MatrixAdd( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *MatrixAdd( 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->eRep != spB->eRep) {
       printf("WARNING : Incompatible Matrix Representations in MatrixAdd()\n");
       printf("WARNING : *** spA->eRep = %4d spB->eRep = %4d\n",
               spA->eRep, spB->eRep);
    }

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

    /* [b] : Compute matrix operation */

    switch((int) spA->eRep) {
        case INDIRECT:
             switch((int) spA->eType) {
                 case DOUBLE_ARRAY:
                      spC = MatrixAddIndirectDouble( spA, spB );
                      break;
                 case INTEGER_ARRAY:
                 case COMPLEX_ARRAY:
                 default:
                      FatalError("In MatrixAdd() : Undefined spA->eType",
                                (char *) NULL);
                      break;
             }
             break;
        case SKYLINE:
             spC = MatrixAddSkyline( spA, spB );
             break;
        case SPARSE:
             break;
        default:
             FatalError("In MatrixAdd() : Undefined spA->eRep",
                       (char *) NULL);
             break;
    }

    return ( spC );
}

/*
 *  ==============================================================
 *  MatrixAddReplace() : 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 *MatrixAddReplace( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *MatrixAddReplace( spA , spB )
MATRIX *spA;
MATRIX *spB;
#endif /* end case not STDC */
{
MATRIX *spC;

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

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

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

    /* [b] : Compute matrix operation */

    switch((int) spA->eRep) {
        case INDIRECT:
             switch((int) spA->eType) {
                 case DOUBLE_ARRAY:
                      spA = MatrixAddReplaceIndirectDouble( spA, spB );
                      break;
                 case INTEGER_ARRAY:
                 case COMPLEX_ARRAY:
                 default:
                      FatalError("In MatrixAddReplace() : Undefined spA->eType",
                                (char *) NULL);
                      break;
             }
             break;
        case SKYLINE:
             spC = MatrixAddSkyline( spA, spB );
             MatrixFreeSkyline( spA );
             return( spC );
             break;
        default:
             FatalError("In MatrixAddReplace() : Undefined spA->eRep",
                       (char *) NULL);
             break;
    }

    return ( spA );
}

/*
 *  ==========================================================
 *  MatrixSub() : 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 *MatrixSub( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *MatrixSub( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* end case not STDC */
{
MATRIX *spC;

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

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

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

    /* [b] : Compute Matrix Operation */

    switch((int) spA->eRep) {
        case INDIRECT:
             switch((int) spA->eType) {
                 case DOUBLE_ARRAY:
                      spC = MatrixSubIndirectDouble( spA, spB );
                      break;
                 case INTEGER_ARRAY:
                 case COMPLEX_ARRAY:
                 default:
                      FatalError("In MatrixSub() : Undefined spA->eType",
                                (char *) NULL);
                      break;
             }
             break;
        case SKYLINE:
             spC = MatrixSubSkyline( spA, spB );
             break;
        case SPARSE:
             break;
        default:
             FatalError("In MatrixSub() : Undefined spA->eRep",
                       (char *) NULL);
             break;
    }

    return ( spC );
}

/*
 *  ==============================================================
 *  MatrixSubReplace() : Compute matrix assignment [A] = [A] - [B]