matrix_indirect.c

有限元程序

MATRIX *spB;
#endif /* End case not STDC */
{
int ii, ij;
int  length1, length2;
DIMENSIONS *d1, *d2;
int UNITS_SWITCH;

     UNITS_SWITCH = CheckUnits();
    /* [a] : Compare Dimensions of Matrices */

       if(MatrixCompareSize( spA, spB ) != TRUE)
          FatalError("Execution halted in MatrixAddIndirectDouble()",
                     "Inconsistent Dimensions",
                     (char *) NULL);

    /* [b] : Check Units and Sub Matrices : [A] = [A] - [B] */

       switch( UNITS_SWITCH) {
         case ON:
            for(ii = 1; ii <= spA->iNoRows; ii++) {
                for(ij = 1; ij <= spA->iNoColumns; ij++) {

                    d1 = UnitsMult(&(spA->spRowUnits[ii-1]),&(spA->spColUnits[ij-1]));
                    d2 = UnitsMult(&(spB->spRowUnits[ii-1]),&(spB->spColUnits[ij-1]));

                    if(SameUnits(d1, d2) == TRUE) 
                      spA->uMatrix.daa[ii-1][ij-1] -= spB->uMatrix.daa[ii-1][ij-1];
                    else {
                     printf("For row No %d, column No %d \n", ii, ij);
                     FatalError("In MatrixAddIndirectDouble(): Inconsistent Units",(char *)NULL);
                    }
                    free((char *) d1->units_name);
                    free((char *) d1);
                    free((char *) d2->units_name);
                    free((char *) d2);
                }
            }
            break;
         case OFF:
            for(ii = 1; ii <= spA->iNoRows; ii++) {
                for(ij = 1; ij <= spA->iNoColumns; ij++) {
                      spA->uMatrix.daa[ii-1][ij-1] -= spB->uMatrix.daa[ii-1][ij-1];
                }
            }
            break;
         default:
            break;
      }

       free((char *) spA->cpMatrixName);
       spA->cpMatrixName = (char *)NULL;
       return (spA);
}

/*
 *  =======================================================================
 *  MatrixMultIndirectDouble() : Matrix Multiplication [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 *MatrixMultIndirectDouble( MATRIX *spA, MATRIX *spB )
#else  /* Start case not STDC */
MATRIX *MatrixMultIndirectDouble( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* End case not STDC */
{
MATRIX    *spC;
int ii, ij, ik;
int  length1, length2, length;
DIMENSIONS *d1, *d2, *d3, *d4, *d5;
int UNITS_SWITCH;

   /* [a] : Compare Inside Dimensions of Matrices */

   if(MatrixCompareInsideDimensions( spA, spB ) != TRUE)
      FatalError("Execution halted in MatrixMultIndirectDouble()",
                 "Inconsistent Inside Dimensions",
                 (char *) NULL);

   /* [b] : Allocate memory for result of matrix multiply */

   spC = MatrixAllocIndirect((char *) NULL, DOUBLE_ARRAY, spA->iNoRows, spB->iNoColumns);
       
   /* [c] : Multiply matrices (with units) */

   UNITS_SWITCH = CheckUnits();

   if (UNITS_SWITCH == ON ) {

       for(ii = 1; ii <= spA->iNoRows; ii++) {
       for(ij = 1; ij <= spB->iNoColumns; ij++) {

           d4 = (DIMENSIONS *) MyCalloc( 1, sizeof(DIMENSIONS) );
           ZeroUnits(d4);

           for(ik = 1; ik <= spA->iNoColumns; ik++) {

               d1 = UnitsMult(&(spA->spRowUnits[ii-1]),&(spA->spColUnits[ik-1]));
               d2 = UnitsMult(&(spB->spRowUnits[ik-1]),&(spB->spColUnits[ij-1]));
               d3 = UnitsMult(d1,d2);

               if( ik==1 && d3->units_name!=(char *)NULL )
                   UnitsCopy(d4, d3);

               if( SameUnits(d3, d4) == TRUE) {
                   spC->uMatrix.daa[ii-1][ij-1] 
                        += spA->uMatrix.daa[ii-1][ik-1]*spB->uMatrix.daa[ik-1][ij-1];
               } else
                   FatalError("In MatrixMultIndirectDouble(): Inconsistent Units",
                             (char *)NULL);

               free((char *) d1->units_name);
               free((char *) d1);
               free((char *) d2->units_name);
               free((char *) d2);
               free((char *) d3->units_name);
               free((char *) d3);
          }
          free((char *) d4->units_name);
          free((char *) d4);
      }
      }

      for(ii = 1; ii <=  spC->iNoRows; ii++)
          UnitsMultRep(&(spC->spRowUnits[ii-1]),&(spB->spRowUnits[0]),&(spA->spRowUnits[ii-1]));
      for(ij = 1; ij <= spC->iNoColumns; ij++)
          UnitsMultRep(&(spC->spColUnits[ij-1]),&(spA->spColUnits[0]),&(spB->spColUnits[ij-1]));

   }

   /* [d] : Multiply matrices (without units) */

   if (UNITS_SWITCH == OFF ) {
       for(ii = 1; ii <= spA->iNoRows; ii++)
       for(ij = 1; ij <= spB->iNoColumns; ij++)
           for(ik = 1; ik <= spA->iNoColumns; ik++)
               spC->uMatrix.daa[ii-1][ij-1]
                    += spA->uMatrix.daa[ii-1][ik-1]*spB->uMatrix.daa[ik-1][ij-1];
   }

   return (spC);
}

/*
 *  =======================================================================
 *  MatrixScaleIndirectDouble() : Multiply Matrix by Scalar
 * 
 *  Input  :  MATRIX *spA        -- Pointer to matrix data structure [A].
 *         :  double scale       -- double : scale factor c.
 *  Output :  MATRIX *spB        -- Pointer to scaled matrix [B] = c.[A].
 *  =======================================================================
 */

#ifdef __STDC__
MATRIX *MatrixScaleIndirectDouble( MATRIX *spA, double dScale)
#else  /* Start case not STDC */
MATRIX *MatrixScaleIndirectDouble( spA, dScale)
MATRIX   *spA;
double dScale;
#endif /* End case not STDC */
{
MATRIX     *spB;
int      ii, ij;

    /* [a] : Check Input matrix [A] */

       if(spA == (MATRIX *) NULL) 
          FatalError("In MatrixScaleIndirectDouble() : [A] == NULL",
                    (char *) NULL);

    /* [b] : Scale Indirect Double Matrix */

       spB = MatrixCopyIndirectDouble( spA );

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

       return ( spB);
}

#ifdef __STDC__
double MatrixContentScaleIndirectDouble(MATRIX *m, int row_no, int col_no )
#else
double MatrixContentScaleIndirectDouble(m, row_no, col_no )
MATRIX  *m;
int row_no;        /* row number    */
int col_no;        /* column number */
#endif
{
int   i, j, UnitsType;
double             da;
DIMENSIONS     *dimen;

#ifdef DEBUG
       printf("*** Enter MatrixContentScaleIndirectDouble() : m->iNoRows    = %4d\n", m->iNoRows);
       printf("                        : m->iNoColumns = %4d\n", m->iNoColumns);
#endif
      if(CheckUnits()==OFF)
         FatalError("You have to set units ON to use this function",
                    "In MatrixContentScaleIndirectDouble",(char *)NULL );
 
      i = row_no;
      j = col_no;

      da = m->uMatrix.daa[i-1][j-1];

#ifdef DEBUG
     printf("\n i = %d, j = %d da = %le \n", i, j, da);
#endif
      if(m->spRowUnits[i-1].scale_factor == 0) {
         printf("==> for Row %d, scale_factor of %s = 0", i, m->spColUnits[i-1].units_name);
         FatalError("Fatal error in MatrixContentScaleIndirectDouble(): ",(char *)NULL);
      }

      if(m->spColUnits[j-1].scale_factor == 0) {
         printf("==> for column %d, scale_factor of %s = 0 \n",j, m->spColUnits[j-1].units_name);
         FatalError("Fatal error in MatrixContentScaleIndirectDouble(): ",(char *)NULL);
      }

      UnitsType = CheckUnitsType();

      if(m->spColUnits[j-1].units_name != NULL) {
         switch(UnitsType) {
           case SI:
             if(!strcmp(m->spColUnits[j-1].units_name, "deg_F") )
                da = ConvertTempUnits(m->spColUnits[j-1].units_name, da, US);
           break;
           case US:
             if(!strcmp(m->spColUnits[j-1].units_name, "deg_C") )
                da = ConvertTempUnits(m->spColUnits[j-1].units_name, da, SI);
           break;
         }
       }

      if(m->spRowUnits[i-1].units_name != NULL) {
         switch(UnitsType) {
           case SI:
             if(!strcmp(m->spRowUnits[i-1].units_name, "deg_F") )
                da = ConvertTempUnits(m->spRowUnits[i-1].units_name, da, US);
           break;
           case US:
             if(!strcmp(m->spRowUnits[i-1].units_name, "deg_C") )
                da = ConvertTempUnits(m->spRowUnits[i-1].units_name, da, SI);
           break;
         }
       }

     /* [a] Scale for Column of matrix m */

          da = da / m->spColUnits[j-1].scale_factor;

     /* [b] Scale for Row of matrix m */

          da = da / m->spRowUnits[i-1].scale_factor;

#ifdef DEBUG
       printf("*** Leave MatrixContentScaleIndirectDouble()\n");
#endif
   return (da);
}

/*
 *  =======================================================================
 *  MatrixNegateIndirectDouble() : Matrix Negation [B] = -[A].
 *  
 *  Input  :  MATRIX     *spA  -- Pointer to matrix A.
 *  Output :  MATRIX     *spB  -- Pointer to matrix B.
 *  =======================================================================
 */

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

     spB = MatrixAllocIndirect((char *) NULL, DOUBLE_ARRAY, spA->iNoRows, spA->iNoColumns);

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

     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]) );
     }
         
     return (spB);
}

/*
 *  =======================================================================
 *  MatrixNegateReplaceIndirectDouble() : Matrix Negation [A] = -[A].
 *  
 *  Input  :  MATRIX     *spA  -- Pointer to matrix A.
 *  Output :  MATRIX     *spB  -- Pointer to matrix A.
 *  =======================================================================
 */

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

     for(ii = 1; ii <= spA->iNoRows; ii = ii + 1)
     for(ij = 1; ij <= spA->iNoColumns; ij = ij + 1)
         spA->uMatrix.daa[ii-1][ij-1] *= -1.0;

     return (spA);
}


/*
 *  =======================================================================
 *  MatrixTransposeIndirectDouble() : Matrix Transpose [B] = [A]^T.
 *  
 *  Input  :  MATRIX *spA  -- Pointer to matrix A.
 *  Output :  MATRIX *spB  -- Pointer to matrix B.
 *  =======================================================================
 */

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

    /* [a] : Transpose Matrix */

       spB = MatrixAllocIndirect((char *) NULL, DOUBLE_ARRAY, spA->iNoColumns, spA->iNoRows);

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

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

/*
 *  ===================================================
 *  MatrixCopyIndirectDouble() : Matrix Copy [B] = [A].
 *  
 *  Input  :  MATRIX *spA  -- Pointer to matrix A.
 *  Output :  MATRIX *spB  -- Pointer to matrix B.
 *  ===================================================
 */

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

     UNITS_SWITCH = CheckUnits();