matrix_skyline.c

有限元程序

     UNITS_SWITCH = CheckUnits();
    /* [a] : Check for compatible matrix dimensions and data types */

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

       if(spA->eRep != spA->eRep) 
          FatalError("In MatrixAddSkyline() : spA->eRep != spB->eRep",
                    (char *) NULL );

       if(spA->eType != spB->eType) 
          FatalError("In MatrixAddSkyline() : spA->eType != spB->eType",
                    (char *) NULL);

    /* [b] : Compute Skyline Profile for [A] + [B] */

       ld = iVectorAlloc( spA->iNoColumns ); 
       for(ii = 0; ii < spA->iNoColumns; ii++) 
           ld[ii] = (int) MAX (spA->uMatrix.daa[ii][0], spB->uMatrix.daa[ii][0]);

    /* [c] : Allocate Skyline Matrix for [C] : [C] = [A] + [B] */

       spC = MatrixAllocSkyline((char *)NULL,DOUBLE_ARRAY,spA->iNoRows, spA->iNoColumns, ld);

    /* [d] : Compute [C] = [A] + [B] */

       switch( UNITS_SWITCH) {
         case ON:
            for(ii = 1; ii <= spA->iNoRows; ii++)
                UnitsCopy( &(spC->spRowUnits[ii-1]), &(spA->spRowUnits[ii-1]) );
            for(ij = 1; ij <= spA->iNoColumns; ij++)
                UnitsCopy( &(spC->spColUnits[ij-1]), &(spA->spColUnits[ij-1]) );

         /* [b] : Check Units and Add Matrices : [C] = [A] + [B] */

            for(ii = 0; ii < spA->iNoColumns; ii++) {
                min_ldi = (int) MIN ( spA->uMatrix.daa[ii][0], spB->uMatrix.daa[ii][0]);
                for(ij = 1; ij <= spA->iNoColumns; ij++) {

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

                    if(SameUnits(d1, d2) == TRUE) {
                       if( ij <= min_ldi )
                          spC->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij] 
                                                   + spB->uMatrix.daa[ii][ij];
                       else if( (spA->uMatrix.daa[ii][0] > min_ldi) && (ij <= ld[ii]) ) 
                          spC->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij];
                       else if( (spB->uMatrix.daa[ii][0] > min_ldi) && (ij <= ld[ii]) ) 
                          spC->uMatrix.daa[ii][ij] = spB->uMatrix.daa[ii][ij];
                    }
                    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 = 0; ii < spA->iNoColumns; ii++) {
                min_ldi = (int) MIN ( spA->uMatrix.daa[ii][0], spB->uMatrix.daa[ii][0]);
                for(ij = 1; ij <= min_ldi; ij++)
                    spC->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij] +
                                               spB->uMatrix.daa[ii][ij];

                if( spA->uMatrix.daa[ii][0] > min_ldi) 
                    for(ij = min_ldi+1; ij <= ld[ii]; ij++) 
                        spC->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij];
                 
                if( spB->uMatrix.daa[ii][0] > min_ldi) 
                    for(ij = min_ldi+1; ij <= ld[ii]; ij++) 
                        spC->uMatrix.daa[ii][ij] = spB->uMatrix.daa[ii][ij];
            }
            break;
         default:
            break;
       }
       spC = MatrixReallocSkyline ( spC );

       free((char *) ld);
       return ( spC );
}

/*
 *  =======================================================================
 *  MatrixSubSkyline() : Subtract two Skyline Matrices
 * 
 *  Input  :  MATRIX *spA        -- Pointer to matrix data structure.
 *         :  MATRIX *spB        -- Pointer to matrix data structure.
 *  Output :  MATRIX *spC        -- Pointer to matrix differnce.
 *  =======================================================================
 */

#ifdef __STDC__
MATRIX *MatrixSubSkyline( MATRIX *spA , MATRIX *spB )
#else  /* Start case not STDC */
MATRIX *MatrixSubSkyline( spA, spB )
MATRIX *spA, *spB;
#endif /* End case not STDC */
{
MATRIX *spC;
int *ld, min_ldi;
int  ii, ij;
int  length, length1, length2;
DIMENSIONS *d1, *d2;
int UNITS_SWITCH;

     UNITS_SWITCH = CheckUnits();
    /* [a] : Check for compatible matrix dimensions and data types */

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

       if(spA->eRep != spA->eRep) 
          FatalError("In MatrixSubSkyline() : spA->eRep != spB->eRep",
                    (char *) NULL );

       if(spA->eType != spB->eType) 
          FatalError("In MatrixSubSkyline() : spA->eType != spB->eType",
                    (char *) NULL);

    /* [b] : Compute Skyline Profile for [A] - [B] */

       ld = iVectorAlloc( spA->iNoRows ); 
       for(ii = 0; ii < spA->iNoRows; ii++) 
           ld[ii] = (int) MAX (spA->uMatrix.daa[ii][0], spB->uMatrix.daa[ii][0]);

    /* [c] : Allocate Skyline Matrix for [C] : [C] = [A] - [B] */

       spC = MatrixAllocSkyline((char *)NULL, DOUBLE_ARRAY,
                        spA->iNoRows, spA->iNoColumns, ld);

    /* [d] : Compute [C] = [A] - [B] */

       switch( UNITS_SWITCH) {
         case ON:
            for(ii = 1; ii <= spA->iNoRows; ii++)
                UnitsCopy( &(spC->spRowUnits[ii-1]), &(spA->spRowUnits[ii-1]) );
            for(ij = 1; ij <= spA->iNoColumns; ij++)
                UnitsCopy( &(spC->spColUnits[ij-1]), &(spA->spColUnits[ij-1]) );

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

            for(ii = 0; ii < spA->iNoRows; ii++) {
                min_ldi = (int) MIN ( spA->uMatrix.daa[ii][0], spB->uMatrix.daa[ii][0]);
                for(ij = 1; ij <= spA->iNoColumns; ij++) {

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

                    if(SameUnits(d1, d2) == TRUE) {
                       if( ij <= min_ldi )
                          spC->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij]
                                                   - spB->uMatrix.daa[ii][ij];
                       else if( (spA->uMatrix.daa[ii][0] > min_ldi) && (ij <= ld[ii]) ) 
                          spC->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij];
                       else if( (spB->uMatrix.daa[ii][0] > min_ldi) && (ij <= ld[ii]) ) 
                          spC->uMatrix.daa[ii][ij] = -spB->uMatrix.daa[ii][ij];
                    }
                    else {
                       printf("For row No %d, column No %d \n", ii, ij);
                       FatalError("In MatrixSubIndirectDouble(): 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 = 0; ii < spA->iNoRows; ii++) {
                min_ldi = (int) MIN ( spA->uMatrix.daa[ii][0], spB->uMatrix.daa[ii][0]);
                for(ij = 1; ij <= min_ldi; ij++)
                    spC->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij] -
                                               spB->uMatrix.daa[ii][ij];

                if( spA->uMatrix.daa[ii][0] > min_ldi) 
                    for(ij = min_ldi+1; ij <= ld[ii]; ij++) 
                        spC->uMatrix.daa[ii][ij] = spA->uMatrix.daa[ii][ij];
                 
                if( spB->uMatrix.daa[ii][0] > min_ldi) 
                    for(ij = min_ldi+1; ij <= ld[ii]; ij++) 
                        spC->uMatrix.daa[ii][ij] = -spB->uMatrix.daa[ii][ij];
            }
            break;
         default:
            break;
       }
       spC = MatrixReallocSkyline ( spC );

       free((char *) ld);
       return ( spC );
}

/*
 *  =======================================================================
 *  MatrixMultSkyline() : Multiply two Skyline Matrices.
 * 
 *  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 *MatrixMultSkyline( MATRIX *spA, MATRIX *spB )
#else  /* Start case not STDC */
MATRIX *MatrixMultSkyline( spA, spB )
MATRIX *spA;
MATRIX *spB;
#endif /* End case not STDC */
{
MATRIX *spC;
int ii, ij, ik, init;
int im, in, kka, kkb;
int  length1, length2, length;
DIMENSIONS *d1, *d2, *d3, *d4, *d5;
int UNITS_SWITCH;

#define DEBUG

#ifdef DEBUG
       printf("*** Enter MatrixMultSkyline()\n");
#endif
     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->eRep != spB->eRep) 
          FatalError("In MatrixMultSkyline() : spA->eRep != spB->eRep",
                    (char *) NULL);

       if(spA->eType != spB->eType) 
          FatalError("In MatrixMultSkyline() : spA->eType != spB->eType",
                    (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) {
         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 = (int) MAX((ii - spA->uMatrix.daa[ii-1][0]+1),
                                   (ij - spB->uMatrix.daa[ij-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) MAX(ik,ii);
                           kka = (int) MIN(ik,ii);

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

                           if((im-kka+1) <= spA->uMatrix.daa[ im-1 ][0] &&
                              (in-kkb+1) <= spB->uMatrix.daa[ in-1 ][0] )
                               spC->uMatrix.daa[ ii-1 ][ ij-1 ] +=
                                  spA->uMatrix.daa[im-1][im-kka+1]
                                 *spB->uMatrix.daa[in-1][in-kkb+1];
                        }
                        else
                          FatalError("In MatrixMultSkyline(): 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 = (int) MAX((ii - spA->uMatrix.daa[ii-1][0]+1),
                                   (ij - spB->uMatrix.daa[ij-1][0]+1));
                  for(ik = init; ik <= spA->iNoColumns; ik++) {
                     im  = (int) MAX(ik,ii);
                     kka = (int) MIN(ik,ii);

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

                     if((im-kka+1) <= spA->uMatrix.daa[im-1][0] &&
                        (in-kkb+1) <= spB->uMatrix.daa[in-1][0] )
                         spC->uMatrix.daa[ ii-1 ][ ij-1 ] +=
                         spA->uMatrix.daa[im-1][im-kka+1]*spB->uMatrix.daa[in-1][in-kkb+1];
                  }
               }
            }
           break;
        default:
           break;
      }

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

      return (spC);
}
#undef DEBUG

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

#ifdef __STDC__
MATRIX *MatrixNegateSkyline( MATRIX *spA )
#else  /* Start case not STDC */
MATRIX *MatrixNegateSkyline( spA )