matrix.c

有限元程序

 *  
 *  Input :    MATRIX spA -- Pointer to Matrix A
 *             MATRIX spB -- Pointer to Matrix B
 *  Output :   MATRIX spA -- Pointer to Replacement Matrix A 
 *  ==============================================================
 */

#ifdef __STDC__
MATRIX *MatrixSubReplace( MATRIX *spA, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *MatrixSubReplace( 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 : *** 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:
                      spA = MatrixSubReplaceIndirectDouble( spA, spB );
                      break;
                 case INTEGER_ARRAY:
                 case COMPLEX_ARRAY:
                 default:
                      FatalError("In MatrixAddReplace() : Undefined spA->eType",
                                (char *) NULL);
                      break;
             }
             break;
        case SKYLINE:
             spC = MatrixSubSkyline( spA, spB );
             MatrixFreeSkyline( spA );
             return( spC );
             break;
        default:
             FatalError("In MatrixSubReplace() : Undefined spA->eRep",
                       (char *) NULL);
             break;
    }

    return ( spA );
}

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

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

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

    /* [b] : Compute Matrix Multiplication for Indirect Storage */

    if((int) spA->eRep == INDIRECT && (int) spB->eRep == INDIRECT) {
       switch((int) spA->eType) {
           case DOUBLE_ARRAY:
                spC = MatrixMultIndirectDouble( spA, spB );
                break;
           case INTEGER_ARRAY:
           case COMPLEX_ARRAY:
           default:
                FatalError("In MatrixMult() : Undefined spA->eType",
                          (char *) NULL);
                break;
       }
    }

    /* [c] : Compute Matrix Multiplication for Skyline-Indirect Storage */

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

    /* [d] : Compute Matrix Multiplication for Skyline-Indirect Storage */

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

    /* [e] : Compute Matrix Multiplication for Skyline-Skyline Storage */

    if((int) spA->eRep == SKYLINE && (int) spB->eRep == SKYLINE) {
       switch((int) spA->eType) {
           case DOUBLE_ARRAY:
       printf(" case : SKYLINE : DOUBLE \n");
                spC = MatrixMultSkyline( spA, spB );
                break;
           case INTEGER_ARRAY:
           case COMPLEX_ARRAY:
           default:
                FatalError("In MatrixMult() : Undefined spA->eType",
                          (char *) NULL);
                break;
       }
    }

    return ( spC );
}

/*
 *  ==========================================================
 *  MatrixPower() : Matrix Power function [c] = [a]^n
 *  
 *  Input :    MATRIX   spA     -- Pointer to Matrix   A
 *             QUANTITY spQ     -- Pointer to Qauntity Q
 *  Output :   MATRIX   spC     -- Pointer to Matrix   C 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *MatrixPower( MATRIX *spA , QUANTITY *spQ )
#else  /* start case not STDC */
MATRIX *MatrixPower( spA, spQ )
MATRIX   *spA;
QUANTITY *spQ;
#endif /* end case not STDC */
{
MATRIX       *spC;
double       dexp; 
int       i, j, n; 
static double EPS = 1.0E-10;

    /* [a] : Check that matrix are square matrix */

    if((int) spA->iNoRows != (int) spA->iNoColumns) {
       printf(" *** In MatrixPower(): \n");
       printf(" *** No of rows of %s    = %d \n", spA->iNoRows, spA->cpMatrixName);
       printf(" *** No of columns of %s = %d \n", spA->iNoColumns, spA->cpMatrixName);
       FatalError(" Matrix must be a quare matrix to use MatrixPower() function",
                 (char *) NULL);
    }

    /* [b] : Check that exponent is an integer */

    dexp = spQ->value;
    n    = (int) dexp;
    if(abs(dexp/((double) n) - 1.0) > EPS) {
       printf(" *** In MatrixPower(): \n");
       printf(" *** ratio  = %lf\n", abs(dexp/(double) n));
       printf(" *** Exponent of = %lf\n", dexp);
       FatalError(" Exponent must be an integer to use MatrixPower() function",
                  (char *) NULL);
    }

    /* [c] : Check that exponent larger than -1 */

    if(dexp < 0.0 && n != (int) -1) {
       printf(" *** In MatrixPower(): \n");
       printf(" *** Exponent of %lf = \n", dexp);
       FatalError(" Exponent must be larger -1 to use MatrixPower() function",
                  (char *) NULL);
    }
     
    /* [d] : Exponent is zero */

    if(n == (int) 0 ) {
       spC = MatrixAllocIndirect((char *) NULL,DOUBLE_ARRAY,spA->iNoRows,spA->iNoColumns);
      if(CheckUnits()==ON) {
         for(i = 1; i <= spA->iNoRows; i++)
             ZeroUnits( &(spC->spRowUnits[i-1]) );
         for(i = 1; i <= spA->iNoColumns; i++)
             ZeroUnits( &(spC->spColUnits[i-1]) );
      }
      for(i = 1; i <= spA->iNoRows; i++)
         for(j = 1; j <= spA->iNoColumns; j++)
             spC->uMatrix.daa[i-1][j-1] = 1.0;
      return(spC);
     }

    /* [e] : Calculate the inverse for Indirect Storge */

     if((int) spA->eRep == INDIRECT && n == (int) -1) {
         switch((int) spA->eType) {
           case DOUBLE_ARRAY:
                spC = MatrixInverseIndirectDouble(spA);
                break;
           case INTEGER_ARRAY:
           case COMPLEX_ARRAY:
           default:
                FatalError("In MatrixPower() : Undefined spA->eType",
                          (char *) NULL);
                break;
         }
     }

    /* [f] : Calculate the inverse for Skyline Storge */

     if((int) spA->eRep == SKYLINE && n == (int) -1) {
         switch((int) spA->eType) {
           case DOUBLE_ARRAY:
                spC = MatrixInverseSkyline(spA);
                break;
           case INTEGER_ARRAY:
           case COMPLEX_ARRAY:
           default:
                FatalError("In MatrixPower() : Undefined spA->eType",
                          (char *) NULL);
                break;
         }
     }


    /* [g] : Compute Matrix Power for Indirect Storage */

    if((int) spA->eRep == INDIRECT && n > 0) {
       switch((int) spA->eType) {
           case DOUBLE_ARRAY:
                if(n == (int) 1)
                   spC = MatrixCopyIndirectDouble(spA);
                else{
                   spC = MatrixCopyIndirectDouble(spA);
                   for(i = 1; i <= n-1; i++) 
                       spC = MatrixMultIndirectDouble( spC, spA );
                }
                break;
           case INTEGER_ARRAY:
           case COMPLEX_ARRAY:
           default:
                FatalError("In MatrixPower() : Undefined spA->eType",
                          (char *) NULL);
                break;
       }
    }

    /* [c] : Compute Matrix Multiplication for Skyline Storage */

    if((int) spA->eRep == SKYLINE && n > 0) {
       switch((int) spA->eType) {
           case DOUBLE_ARRAY:
                if(n == (int) 1)
                   spC = MatrixCopySkyline(spA);
                else{
                   spC = MatrixCopySkyline(spA);
                   for(i = 1; i <= n-1; i++) 
                       spC = MatrixMultSkyline(spC, spA );
                }
                break;
           case INTEGER_ARRAY:
           case COMPLEX_ARRAY:
           default:
                FatalError("In MatrixPower() : Undefined spA->eType",
                          (char *) NULL);
                break;
       }
    }

    return ( spC );
}

/*
 *  ==========================================================
 *  MatrixNegate() : Compute Negative Matrix [B] = -[A];
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *  Output :   MATRIX spB       -- Pointer to Matrix B 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *MatrixNegate( MATRIX *spA )
#else  /* start case not STDC */
MATRIX *MatrixNegate( spA )
MATRIX *spA;
#endif /* end case not STDC */
{
MATRIX *spB;

    switch((int) spA->eRep) {
        case INDIRECT:
             switch((int) spA->eType) {
                 case DOUBLE_ARRAY:
                      spB = MatrixNegateIndirectDouble( spA );
                      break;
                 case INTEGER_ARRAY:
                 case COMPLEX_ARRAY:
                 default:
                      FatalError("In MatrixNegate() : Undefined spA->eType",
                                (char *) NULL);
                      break;
             }
             break;
        case SKYLINE:
             spB = MatrixNegateSkyline( spA );
             break;
        default:
             FatalError("In MatrixNegate() : Undefined spA->eRep",
                       (char *) NULL);
             break;
    }

    return ( spB );
}

/*
 *  ==========================================================
 *  MatrixNegateReplace() : Matrix Replacement [A] = -[A];
 *  
 *  Input :    MATRIX spA       -- Pointer to Matrix A
 *  Output :   MATRIX spA       -- Pointer to Matrix -A 
 *  ==========================================================
 */

#ifdef __STDC__
MATRIX *MatrixNegateReplace( MATRIX *spA )
#else  /* start case not STDC */
MATRIX *MatrixNegateReplace( spA )
MATRIX *spA;
#endif /* end case not STDC */
{

    switch((int) spA->eRep) {
        case INDIRECT:
             switch((int) spA->eType) {
                 case DOUBLE_ARRAY:
                      spA = MatrixNegateReplaceIndirectDouble(spA);
                      break;
                 default:
                      FatalError("In MatrixNegateReplace() : Undefined m->eType",
                                (char *) NULL);
                      break;
             }
             break;
        case SKYLINE:
             spA = MatrixNegateReplaceSkyline( spA );
             break;
        default:
             FatalError("In MatrixNegateReplace() : Undefined spA->eRep",
                       (char *) NULL);
             break;
    }

    return (spA);
}

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

#ifdef __STDC__
MATRIX *MatrixTranspose( MATRIX *spA )
#else  /* start case not STDC */
MATRIX *MatrixTranspose( spA )
MATRIX *spA;
#endif /* end case not STDC */
{
MATRIX *spB;

    switch((int) spA->eRep) {
        case INDIRECT:
             switch((int) spA->eType) {
                 case DOUBLE_ARRAY:
                      spB = MatrixTransposeIndirectDouble( spA );
                      break;
                 case INTEGER_ARRAY:
                 case COMPLEX_ARRAY:
                 default:
                      FatalError("In MatrixTranspose() : Undefined spA->eType",
                                (char *) NULL);
                      break;