matrix.c

有限元程序

             }
             break;
        case SKYLINE:
             spB = MatrixTransposeSkyline( spA );
             break;
        case SPARSE:
             break;
        default:
             FatalError("In MatrixTranspose() : Undefined spA->eRep",
                       (char *) NULL);
             break;
    }

    return ( spB );
}

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

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

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

    return ( spB );
}


/*
 *  ==========================================================================
 *  MatrixSolve() : Solve Systems of Linear Equations [A][X] = [B].
 * 
 *  Input :    MATRIX spA -- Pointer to Matrix A
 *             MATRIX spB -- Pointer to r-h-s Matrix B
 *  Output :   MATRIX spX -- Pointer to solution matrix X
 * 
 *  Also calls functions : MatrixLU() - LU decomposition.
 *                         MatrixFB() - Forward/Back substitution.
 *
 *  WARNING : This high-level equation solver works for only one decomposition
 *            followed by multiple forward/backward substitutions.
 *  ==========================================================================
 */

static VECTOR *spPivot = (VECTOR *)NULL;

#ifdef __STDC__
MATRIX *MatrixSolve( MATRIX *spA , MATRIX *spB )
#else  /* start case not STDC */
MATRIX *MatrixSolve( spA, spB )
MATRIX      *spA, *spB;
#endif /* end case not STDC */
{
MATRIX    *spLU;
MATRIX     *spX;

       /* [a] LUP Decomposition followed by Forward/Backward Substitution */

       switch((int) spA->eRep ) {
           case INDIRECT:
                switch((int) spA->eType) {
                    case DOUBLE_ARRAY:
                         if(spPivot != (VECTOR *)NULL) {
                            VectorFreeDouble(spPivot);
                            spPivot = (VECTOR *)NULL;
                         }
                         spPivot = SetupPivotVector( spA );
                         spLU    = LUDecompositionIndirect( spA , spPivot);
                         spX     = LUSubstitutionIndirect( (char *)NULL, spPivot, spLU, spB);
                         break;
                    default:
                         FatalError("In MatrixSolve() : Undefined A->eType", (char *) NULL);
                         break;
                }
                break;
           case SKYLINE:
                spLU = LUDecompositionSkyline( spA );
                spX  = LUBacksubstitutionSkyline( spLU, spB );
                break;
           case SPARSE:
                break;
          default:
                FatalError("In MatrixSolve() : Undefined spA->eRep",
                          (char *) NULL);
                break;
       }

       MatrixFree( spLU );
       return ( spX );
}

/*
 *  ========================================================================
 *  MatrixLU() : Compute LU decompostion with Pivoting of Rows and U_ii = 1.
 * 
 *  Input :    MATRIX spA  -- Pointer to Matrix A
 *  Output :   MATRIX spLU -- Pointer to decomposed matrix LU
 *  ========================================================================
 */


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

       /* [a] : LUP Decomposition */

       switch((int) spA->eRep) {
           case INDIRECT:
                switch((int) spA->eType) {
                    case DOUBLE_ARRAY:
                         if(spPivot != (VECTOR *)NULL) {
                            VectorFreeDouble(spPivot);
                            spPivot = (VECTOR *)NULL;
                         }
                         spPivot = SetupPivotVector( spA );
                         spLU    = LUDecompositionIndirect( spA , spPivot);
                         break;
                    default:
                         FatalError("In MatrixLU() : Undefined A->eType", (char *) NULL);
                         break;
                }
                break;
           case SKYLINE:
                spLU = LUDecompositionSkyline( spA );
                break;
           case SPARSE:
                break;
          default:
                FatalError("In MatrixLU() : Undefined spA->eRep",
                          (char *) NULL);
                break;
       }

       return ( spLU );
}

/*
 *  =======================================================
 *  MatrixFB() : Compute forward and backward substitution.
 * 
 *  Input :    MATRIX spLU -- Pointer to Matrix A
 *             MATRIX spB  -- Pointer to r.h.s. matrix B
 *  Output :   MATRIX spX  -- Pointer to solution matrix X
 *  =======================================================
 */

#ifdef __STDC__
MATRIX *MatrixFB( MATRIX *spLU, MATRIX *spB )
#else  /* start case not STDC */
MATRIX *MatrixFB( spLU, spB )
MATRIX *spLU;
MATRIX  *spB;
#endif /* end case not STDC */
{
MATRIX *spX;

       /* [a] : Forward/Backward Substitution */

       switch((int) spLU->eRep) {
           case INDIRECT:
                switch((int) spLU->eType) {
                    case DOUBLE_ARRAY:
                         spX = LUSubstitutionIndirect((char *) NULL, spPivot, spLU, spB);
                         break;
                    default:
                         FatalError("In MatrixFB() : Undefined lu->eType",
                                   (char *) NULL);
                         break;
                }
                break;
           case SKYLINE:
                spX = LUBacksubstitutionSkyline( spLU, spB );
                break;
          default:
                FatalError("In MatrixFB() : Undefined spLU->eRep",
                          (char *) NULL);
                break;
       }

       return ( spX );
}

/*
 *  =============================================================
 *  MatrixInverse() : Compute matrix inverse
 * 
 *  Input :    MATRIX spA    -- Pointer to Matrix A
 *  Output :   MATRIX spAinv -- Pointer to inversed matrix spAinv
 *  =============================================================
 */

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

       switch((int) spA->eRep) {
           case INDIRECT:
                switch((int) spA->eType) {
                    case DOUBLE_ARRAY:
                         spAinv = MatrixInverseIndirectDouble( spA );
                         break;
                    default:
                         FatalError("In MatrixInverse() : Undefined spA->eType",
                                   (char *) NULL);
                         break;
                }
                break;
           case SKYLINE:
                spAinv = MatrixInverseSkyline( spA );
                break;
          default:
                FatalError("In MatrixInverse() : Undefined spA->eRep",
                          (char *) NULL);
                break;
       }

       return ( spAinv );
}

/*
 *  =======================================================
 *  MatrixDet() : Compute determinant of Matrix
 * 
 *  Input :    MATRIX     m  -- Pointer to Matrix m
 *  Output :   QUANTITY   q  -- Matrix Determinant.
 *  =======================================================
 */

#ifdef __STDC__
QUANTITY  *MatrixDet(MATRIX *m)
#else
QUANTITY  *MatrixDet(m)
MATRIX    *m;
#endif
{
MATRIX   *m1;
QUANTITY  *q;

#ifdef DEBUG
       printf("*** Enter MatrixDet() : m->iNoRows    = %4d\n", m->iNoRows);
       printf("                      : m->iNoColumns = %4d\n", m->iNoColumns);
#endif
      q         = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));
      if(CheckUnits() == ON) {
         q->dimen  = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
         ZeroUnits(q->dimen);
      }
      else
         q->dimen = (DIMENSIONS *)NULL;

      m1 = MatrixCopy(m);
      if( m1->eRep==SKYLINE )   m1 = MatrixSkylineToIndirect(m1);
          
      q->value = dMatrixDet( m1->uMatrix.daa, m1->iNoRows, m1->iNoColumns );

      MatrixFree(m1);

#ifdef DEBUG
       printf("*** Leave MatrixDet()\n");
#endif

   return (q);
}

/*
 *  ===========================================================
 *  MatrixMax() : Get the maximum number in the matrix elements
 *  MatrixMin() : Get the minimum number in the matrix elements
 * 
 *  Input :    MATRIX     m  -- Pointer to Matrix m
 *  Output :   QUANTITY   q  -- Max or Min Number
 *  ===========================================================
 */

#ifdef __STDC__
QUANTITY  *MatrixMax(MATRIX *m)
#else
QUANTITY  *MatrixMax(m)
MATRIX    *m;
#endif
{
int     i, j;
QUANTITY  *q;

      q  = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));
      if(CheckUnits() == ON) {
         q->dimen  = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
         ZeroUnits(q->dimen);
      }
      else
         q->dimen = (DIMENSIONS *)NULL;

      q->value = m->uMatrix.daa[0][0];
      for( i=1 ; i <= m->iNoRows ; i++ ) {
          for( j=1 ; j <= m->iNoColumns ; j++ ) {
              q->value = MAX( q->value, m->uMatrix.daa[i-1][j-1] );
          }
      }

   return (q);
}

#ifdef __STDC__
QUANTITY  *MatrixMin(MATRIX *m)
#else
QUANTITY  *MatrixMin(m)
MATRIX    *m;
#endif
{
int     i, j;
QUANTITY  *q;

      q = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));
      if(CheckUnits() == ON) {
         q->dimen  = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
         ZeroUnits(q->dimen);
      }
      else
         q->dimen = (DIMENSIONS *)NULL;

      q->value = m->uMatrix.daa[0][0];
      for( i=1 ; i <= m->iNoRows ; i++ ) {
          for( j=1 ; j <= m->iNoColumns ; j++ ) {
              q->value = MIN( q->value, m->uMatrix.daa[i-1][j-1] );
          }
      }

   return (q);
}

QUANTITY *MatrixL2Norm(m)
MATRIX   *m;
{
QUANTITY *q;
double   sum = 0.0;
int i;

#ifdef DEBUG
       printf("*** Enter MatrixL2Norm() : m->iNoRows    = %4d\n", m->iNoRows);
       printf("                         : m->iNoColumns = %4d\n", m->iNoColumns);
#endif

   /* Check that matrix is either a column (or row) vector */

      if((m->iNoRows != 1) && (m->iNoColumns != 1)) {
          FatalError("In MatrixL2Norm() : Matrix is not a row (or column) vector",(char *)NULL);
      }

   /* Compute L2 Norm for column/row vector */

   sum       =  dVmatrixL2Norm(m->uMatrix.daa, m->iNoRows, m->iNoColumns);
   q         = (QUANTITY *) MyCalloc(1,sizeof(QUANTITY));
   q->value  = sum;
   if(CheckUnits() == ON) {
      q->dimen  = (DIMENSIONS *) MyCalloc(1,sizeof(DIMENSIONS));
      ZeroUnits(q->dimen);
   } else
      q->dimen  = (DIMENSIONS *)NULL;

   return (q);
}

/* 
 *  =================
 *  Matrix Dimensions 
 *  ================= 
 */ 

#ifdef __STDC__
MATRIX *MatrixDimension(MATRIX *m1)
#else
MATRIX *MatrixDimension(m1)
MATRIX *m1;
#endif
{
MATRIX *m2;

      m2 = MatrixAllocIndirect("Dimensions", DOUBLE_ARRAY, 1, 2);