matrix_skyline.c

有限元程序

MATRIX *spA;
#endif /* End case not STDC */
{
MATRIX *spB;
int  ii, ij;
int     *ld;
int  length;

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

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

    /* [b] : Negate Skyline Matrix */
    
       ld = iVectorAlloc ( spA->iNoRows); 
       for(ii = 1; ii <= spA->iNoRows; ii++)
           ld[ii-1] = spA->uMatrix.daa[ii-1][0];
    
       spB = MatrixAllocSkyline( (char *)NULL, DOUBLE_ARRAY,
                                    spA->iNoRows, spA->iNoColumns, ld);

       for(ii = 1; ii <= spA->iNoRows; ii++)
           for(ij = 1; ij <= ld[ii-1]; ij++)
               spB->uMatrix.daa[ii-1][ij] = -spA->uMatrix.daa[ii-1][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]));
       }
           
     free((char *) ld);
     return (spB);
}

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

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

     return (spA);
}

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

    /* [a] : Transpose Matrix */

       ld = iVectorAlloc ( spA->iNoColumns); 
       for(ii = 1; ii <= spA->iNoColumns; ii++)
           ld[ii-1] = spA->uMatrix.daa[ii-1][0];
    
       spB = MatrixAllocSkyline( (char *)NULL, DOUBLE_ARRAY,
                                    spA->iNoRows, spA->iNoColumns, ld);

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

       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]));
       }
      
       free((char *) ld);
       return (spB);
}


/*
 *  ===============================================================================
 *  MatrixMultSkylineIndirectDouble() : Compute Product of Skyline-Indirect Matrix.
 * 
 *  Input  :  MATRIX *spA -- Pointer to matrix data structure [A].
 *         :  MATRIX *spB -- Pointer to matrix data structure [B].
 *  Output :  MATRIX *spC -- Pointer to matrix product [C] = [A]*[B]
 * 
 *  Note : If [A] and [B] are symmetric, then [A]*[B] is symmetric only iff
 *         [A]*[B] = [B]*[A]. In general, this condition will not hold -- so
 *         we store matrix product [ma]*[mb] in INDIRECT form.
 *  ===============================================================================
 */

#ifdef __STDC__
MATRIX *MatrixMultSkylineIndirectDouble( MATRIX *spA, MATRIX *spB )
#else  /* Start case not STDC */
MATRIX *MatrixMultSkylineIndirectDouble( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* End case not STDC */
{
MATRIX *spC;
double       dAvalue;
int ii, ij, ik, init;
int   im,   in;
int  length1, length2, length;
DIMENSIONS *d1, *d2, *d3, *d4, *d5;
int UNITS_SWITCH;

     UNITS_SWITCH = CheckUnits();

    /* [a] : Check for compatible matrix dimensions and data types etc... */

       if((spA == NULL) || (spB == NULL))
          FatalError("In MatrixMultSkyline() : Either spA = NULL or spB = NULL",
                    (char *) NULL);

       if(spA->iNoColumns != spB->iNoRows) 
          FatalError("In MatrixMultSkyline() : spA->iNoColumns != spB->iNoRows",
                    (char *) NULL);

    /* [b] : Allocate Memory for Matrix Product */

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

    /* [c] : Compute elements of matrix product */

       switch( UNITS_SWITCH == ON ) {
         case ON:
            for(ii = 1; ii <= spA->iNoRows; ii++) {
               for(ij = 1; ij <= spB->iNoColumns; ij++) {
                  d4 = (DIMENSIONS *) MyCalloc( 1, sizeof(DIMENSIONS) );
                  ZeroUnits(d4);

                  init = ii - spA->uMatrix.daa[ii-1][0] + 1;

                  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( ik >= init ) {
                         if( SameUnits(d3, d4) == TRUE) {

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

                           if((in-im+1) <= spA->uMatrix.daa[in-1][0])
                               dAvalue = spA->uMatrix.daa[ in-1 ][ in-im+1 ];
                           else
                               dAvalue = 0.0;

                           spC->uMatrix.daa[ ii-1 ][ ij-1 ] +=
                                dAvalue * spB->uMatrix.daa[ ik - 1 ][ ij - 1 ];
                         } 
                         else
                           FatalError(" In MatrixMultSkylineIndirectDouble(): 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]) );

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

              init = ii - spA->uMatrix.daa[ii-1][0] + 1;
              for(ik = init; ik <= spA->iNoColumns; ik++) {
        
                  im = MIN( ii, ik );
                  in = MAX( ii, ik );

                  if((in-im+1) <= spA->uMatrix.daa[in-1][0])
                      dAvalue = spA->uMatrix.daa[ in-1 ][ in-im+1 ];
                  else
                      dAvalue = 0.0;

                  spC->uMatrix.daa[ ii-1 ][ ij-1 ] +=
                       dAvalue * spB->uMatrix.daa[ ik - 1 ][ ij - 1 ];
            
             }
          }
          }
           break;
        default:
           break;
      }

#ifdef DEBUG
       printf("*** Leaving MatrixMultSkylineIndirectDouble()\n");
#endif

      return (spC);
}

#ifdef __STDC__
MATRIX *MatrixMultIndirectSkylineDouble( MATRIX *spA, MATRIX *spB )
#else  /* Start case not STDC */
MATRIX *MatrixMultIndirectSkylineDouble( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* End case not STDC */
{
MATRIX *spC;
double       dBvalue;
int ii, ij, ik, init;
int   im,   in;
int  length1, length2, length;
DIMENSIONS *d1, *d2, *d3, *d4, *d5;
int UNITS_SWITCH;

     UNITS_SWITCH = CheckUnits();

    /* [a] : Check for compatible matrix dimensions and data types etc... */

       if((spA == NULL) || (spB == NULL))
          FatalError("In MatrixMultSkylineIndirect() : Either spA = NULL or spB = NULL",
                    (char *) NULL);

       if(spA->iNoColumns != spB->iNoRows) 
          FatalError("In MatrixMultSkylineIndirect() : spA->iNoColumns != spB->iNoRows",
                    (char *) NULL);

    /* [b] : Allocate Memory for Matrix Product */

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

    /* [c] : Compute elements of matrix product */

       switch( UNITS_SWITCH == ON ) {
         case ON:
            for(ii = 1; ii <= spA->iNoRows; ii++) {
               for(ij = 1; ij <= spB->iNoColumns; ij++) {
                  d4 = (DIMENSIONS *) MyCalloc( 1, sizeof(DIMENSIONS) );
                  d4 = ZeroUnits(d4);

                  init = ii - spB->uMatrix.daa[ii-1][0] + 1;

                  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( ik >= init ) {
                         if( SameUnits(d3, d4) == TRUE) {

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

                           if((in-im+1) <= spB->uMatrix.daa[in-1][0])
                               dBvalue = spB->uMatrix.daa[ in-1 ][ in-im+1 ];
                           else
                               dBvalue = 0.0;

                           spC->uMatrix.daa[ ii-1 ][ ij-1 ] +=
                                spA->uMatrix.daa[ ii - 1 ][ ik - 1 ]* dBvalue;
                         } 
                         else
                           FatalError(" In MatrixMultIndirectSkylineDouble(): 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]));

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

               init = ii - spB->uMatrix.daa[ii-1][0] + 1;
               for(ik = init; ik <= spB->iNoRows; ik++) {
        
                   im = MIN( ik, ij );
                   in = MAX( ik, ij );

                   if((in-im+1) <= spB->uMatrix.daa[in-1][0])
                       dBvalue = spB->uMatrix.daa[ in-1 ][ in-im+1 ];
                   else
                       dBvalue = 0.0;

                   spC->uMatrix.daa[ ii-1 ][ ij-1 ] +=
                        spA->uMatrix.daa[ ii - 1 ][ ik - 1 ]* dBvalue;
            
              }
           }
           }
           break;
        default:
           break;
      }

#ifdef DEBUG
       printf("*** Leaving MatrixMultIndirectSkylineDouble()\n");
#endif

      return (spC);
} 

/*
 *  =======================================================================
 *  MatrixCopySkyline() : Copy Skyline Matrices
 * 
 *  Input  :  MATRIX *spA        -- Pointer to matrix data structure.
 *  Output :  MATRIX *spCopy     -- Pointer to matrix copy
 *  =======================================================================
 */

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