matrix_indirect.c

有限元程序

    /* [a] : Make Copy of Matrix */

       spB = MatrixAllocIndirect( spA->cpMatrixName, DOUBLE_ARRAY, spA->iNoRows, spA->iNoColumns);

       switch( UNITS_SWITCH) {
         case ON:
             for(ii = 1; ii <= spA->iNoRows; ii++) {
                 UnitsCopy(&(spB->spRowUnits[ii-1]), &(spA->spRowUnits[ii-1]));
                 for(ij = 1; ij <= spA->iNoColumns; ij++) {
                    spB->uMatrix.daa[ii-1][ij-1] = spA->uMatrix.daa[ii-1][ij-1];
                 }
             }
             for(ij = 1; ij <= spA->iNoColumns; ij++)
                 UnitsCopy(&(spB->spColUnits[ij-1]), &(spA->spColUnits[ij-1]));
             break;
         case OFF:
             for(ii = 1; ii <= spA->iNoRows; ii++)
                 for(ij = 1; ij <= spA->iNoColumns; ij++)
                    spB->uMatrix.daa[ii-1][ij-1] = spA->uMatrix.daa[ii-1][ij-1];
             break;
         default:
             break;
       }

       return ( spB );
}


#ifdef __STDC__
MATRIX *MatrixInverseIndirectDouble ( MATRIX *spA )
#else  /* Start case not STDC */
MATRIX *MatrixInverseIndirectDouble ( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
MATRIX     *spAinv;
MATRIX *lu, *F, *X;
double          *d;
VECTOR       *indx;
int           i, j;
int         length;
int   UNITS_SWITCH;
     
     UNITS_SWITCH = CheckUnits();
     /* [0] : Check whether matrix spA is a square matrix    */
        
        if(spA->iNoRows != spA->iNoColumns) {
           printf("******matrix %s is not a square matrix :\n", spA->cpMatrixName);
           printf("******       %s->iNoRows    = %d\n", spA->cpMatrixName, spA->iNoRows); 
           printf("******       %s->iNoColumns = %d\n", spA->cpMatrixName, spA->iNoColumns); 
           printf("******check input file \n");
           FatalError("****** Only the square matrix has inverse. In MatrixInverse()",(char *)NULL);
        }

     /* [a] : Set up the column vector of ones */

     F    = MatrixAllocIndirect("F", DOUBLE_ARRAY, spA->iNoRows, 1);
     for(j = 1; j<= spA->iNoRows; j++)
            F->uMatrix.daa[j-1][0] = 0.0;

     if( UNITS_SWITCH == ON ) {
        ZeroUnits(&(F->spColUnits[0]));
        for(i = 1; i <= spA->iNoRows; i++)
            ZeroUnits(&(F->spRowUnits[i-1]));
     }

     /* [b] : Compute [L][U] decomposition of spA (copy) */

     indx = SetupPivotVector(spA);
     lu   = LUDecompositionIndirect(spA,indx);

     /* [c] : Solve series of equations for inverse */

     spAinv = MatrixAllocIndirect("Inverse(spA)", DOUBLE_ARRAY, spA->iNoRows, spA->iNoColumns);
     for(i = 1; i <= spA->iNoColumns; i++) {
         F->uMatrix.daa[i-1][0] = 1.0;
         X = LUSubstitutionIndirect((char *)NULL,indx,lu,F);
         for(j = 1; j<= spA->iNoRows; j++)
             spAinv->uMatrix.daa[j-1][i-1] = X->uMatrix.daa[j-1][0];
         MatrixFreeIndirect(X);
         F->uMatrix.daa[i-1][0] = 0.0;
     }

     /* [d] : Assign units to spAinv */
     if( UNITS_SWITCH == ON ) {
        for(i = 1; i <= spA->iNoRows; i++){
           UnitsPowerRep( &(spAinv->spColUnits[i-1]), &(spA->spRowUnits[i-1]), -1.0, YES );
           UnitsPowerRep( &(spAinv->spRowUnits[i-1]), &(spA->spColUnits[i-1]), -1.0, YES );
        }
     }

     /* [e] : Free memory, return result */

     VectorFree(indx);
     MatrixFree( lu );
     MatrixFree( F );

     return (spAinv);
}

/* 
 *  =====================================================================
 *  LUDecompositionDouble() : Use method of Crout Reduction with pivoting
 *  and scaled equtions to decompose a (nxn) matrix [A] into [L][U].
 *  
 *  Input  :  MATRIX     *spA  -- Pointer to (nxn) matrix A.
 *         :  VECTOR *spPivot  -- Pointer to (nx1) pivot vector.
 *  Output :  MATRIX     *spA  -- Pointer to (nxn) [L][U] product.
 *  =====================================================================
 */

#ifdef __STDC__
MATRIX *LUDecompositionIndirect( MATRIX *spA, VECTOR *spPivot)
#else  /* Start case not STDC */
MATRIX *LUDecompositionIndirect( spA, spPivot )
MATRIX      *spA;
VECTOR *spPivot;
#endif /* End case not STDC */
{
MATRIX         *spLU;
VECTOR      *spScale;
double          dSum;
double    dNumerator;
double  dDenominator;
int iOrder1, iOrder2;
int       ii, ij, ik;

    /* [a] : Make sure that Matrix spA is square    */

       if((spA->iNoRows != spA->iNoColumns)) {
           (void) printf("Problem : spA->iNoRows = %4d spA->iNoColumns = %4d\n",
                                    spA->iNoRows,      spA->iNoColumns);
           FatalError("Execution halted in LUDecomposition()",
                      "Matrix spA must be square",
                      (char *) NULL);
       }

    /* [b] : Setup Matrix of Scale Factors */

       spLU    = MatrixCopyIndirectDouble( spA );
       spScale = SetupScaleFactors( spLU );
    
    /* [c] : Loop over Columns and Compute Crout Reduction */

       ij = 1;
       Pivot( spLU, spScale, spPivot, ij );
       for (ij = 2; ij <= spLU->iNoColumns; ij = ij + 1) {
            iOrder1 = spPivot->uVector.ia[ 0 ];

            spLU->uMatrix.daa[iOrder1-1][ij - 1] =
                 spLU->uMatrix.daa[iOrder1-1][ij - 1]/spLU->uMatrix.daa[iOrder1-1][ 0 ];
       }

       for (ij = 2; ij <= spLU->iNoColumns - 1; ij = ij + 1) {
            for (ii = ij; ii <= spLU->iNoColumns; ii = ii + 1) {

                 iOrder1 = spPivot->uVector.ia[ ii - 1 ];

                 dSum = 0.0;
                 for (ik = 1; ik <= ij - 1; ik = ik + 1) {
                      iOrder2 = spPivot->uVector.ia[ ik - 1 ];
                      dSum += spLU->uMatrix.daa[ iOrder1-1 ][ ik - 1 ] *
                              spLU->uMatrix.daa[ iOrder2-1 ][ ij - 1 ];
                 }

                 spLU->uMatrix.daa[ iOrder1-1 ][ ij - 1 ] -= dSum;
            }

            Pivot( spLU, spScale, spPivot, ij );

            iOrder1 = spPivot->uVector.ia[ ij - 1 ];

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

                 dSum = 0.0;
                 for (ii = 1; ii <= ij - 1; ii = ii + 1) {
                      iOrder2 = spPivot->uVector.ia[ ii - 1 ];
                      dSum += spLU->uMatrix.daa[ iOrder1-1 ][ ii - 1 ] *
                              spLU->uMatrix.daa[ iOrder2-1 ][ ik - 1 ];
                 }

                 dNumerator   = spLU->uMatrix.daa[ iOrder1-1 ][ ik-1 ] - dSum;
                 dDenominator = spLU->uMatrix.daa[ iOrder1-1 ][ ij-1 ];

                 spLU->uMatrix.daa[ iOrder1-1 ][ ik - 1 ] = dNumerator/dDenominator;
            }
       }

       iOrder1 = spPivot->uVector.ia[ spLU->iNoColumns - 1 ];

       dSum = 0.0;
       for (ik = 1; ik <= spLU->iNoColumns - 1; ik = ik + 1) {
            iOrder2 = spPivot->uVector.ia[ ik - 1 ];
            dSum += spLU->uMatrix.daa[ iOrder1-1 ][ ik - 1 ] *
                    spLU->uMatrix.daa[ iOrder2-1 ][ spLU->iNoColumns - 1 ];
       }

       spLU->uMatrix.daa[ iOrder1 - 1 ][ spLU->iNoColumns - 1 ] -= dSum;

       VectorFree( spScale );
       return( spLU );
}


/* 
 *  =====================================================================
 *  LUSubstitutionDouble() : Compute Forward and Backward Substitution.
 *  
 *  Input  :  cpMatrixNameOfSolution -- Name of Solution Matrix
 *         :  VECTOR *spPivot        -- Pointer to (nx1) pivot vector.
 *         :  MATRIX    *spLU        -- Pointer to (nxn) matrix [L][U].
 *         :  MATRIX     *spB        -- Pointer to (nx1) matrix [B].
 *  Output :  MATRIX      spX        -- Pointer to (nx1) solution matrix.
 *  =====================================================================
 */

#ifdef __STDC__
MATRIX *LUSubstitutionIndirect( char *cpMatrixNameOfSolution, VECTOR *spPivot, MATRIX *spLU, MATRIX *spB )
#else  /* Start case not STDC */
MATRIX *LUSubstitutionIndirect( cpMatrixNameOfSolution, spPivot, spLU, spB)
char *cpMatrixNameOfSolution;
VECTOR              *spPivot;
MATRIX                 *spLU;
MATRIX                  *spB;
#endif /* End case not STDC */
{
MATRIX            *spDisp;
double       dDenominator;
double         dNumerator;
double               dSum;
int       iOrder1, ii, ij;
int                length;
DIMENSIONS   *d, *d1, *d2;
int UNITS_SWITCH;

     UNITS_SWITCH = CheckUnits();
   /* [a] : Allocate Memory for Displacement Vector */

      spDisp = MatrixAllocIndirect( cpMatrixNameOfSolution, DOUBLE_ARRAY, spLU->iNoRows, 1 );

   /* [0] : Calculate the units of displacements */

   if( UNITS_SWITCH == ON ) {
      for (ii = 1; ii <= spLU->iNoRows; ii++) {
         d1 = UnitsMult( &(spB->spRowUnits[ii-1]), &(spB->spColUnits[0]) );
         d2 = UnitsMult( &(spLU->spRowUnits[ii-1]), &(spLU->spColUnits[ii-1]) );
         UnitsDivRep( &(spDisp->spRowUnits[ii-1]), d1, d2, YES); 

         free((char *) d1->units_name);
         free((char *) d1);
         free((char *) d2->units_name);
         free((char *) d2);
      }
      ZeroUnits(&(spDisp->spColUnits[0]));
   }

   /* [b] : Forward Substitution with Unscrambling of Permutations */

      iOrder1 = spPivot->uVector.ia[ 0 ];
      spDisp->uMatrix.daa[0][0] = spB->uMatrix.daa[ iOrder1 - 1 ][ 0 ] /
                                  spLU->uMatrix.daa[ iOrder1 - 1 ][ 0 ];

      for( ii = 2; ii <= spLU->iNoRows; ii = ii + 1) {

           iOrder1 = spPivot->uVector.ia[ ii - 1 ];

           dSum = 0.0;
           for( ij = 1; ij <= ii - 1; ij = ij + 1) {
                dSum += spLU->uMatrix.daa[ iOrder1 - 1 ][ ij - 1] *
                        spDisp->uMatrix.daa[ ij - 1][0];
           }               

           dNumerator   = spB->uMatrix.daa[ iOrder1-1 ][0] - dSum;
           dDenominator = spLU->uMatrix.daa[ iOrder1-1 ][ ii-1 ];

           spDisp->uMatrix.daa[ ii - 1 ][0] = dNumerator/dDenominator;
      }               

   /* [c] : Backsubstitution with Unscrambling of Permutations */

      for( ii = spLU->iNoRows - 1; ii >= 1; ii = ii - 1) {

           iOrder1 = spPivot->uVector.ia[ ii - 1 ];

           dSum = 0.0;
           for( ij = ii + 1; ij <= spLU->iNoRows; ij = ij + 1) {
                dSum += spLU->uMatrix.daa[ iOrder1 - 1 ][ ij - 1] *
                        spDisp->uMatrix.daa[ ij - 1][0];
           }               

           spDisp->uMatrix.daa[ ii - 1 ][0] -= dSum;
      }               

      return(spDisp);
}


/*
 *  ================================================================
 *  CholeskyDecompositionIndirect() : [L][L]^T Cholesky Factorization 
 *                                    of Symmetric Indirect Matrix.
 * 
 *  Input  :  MATRIX *spA   -- Pointer to symmetric matrix [A].
 *  Output :  MATRIX *spA   -- Pointer to [L][L]^T.
 *  ================================================================
 */

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

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

       if( spA == NULL )
           FatalError("In CholeskyDecompositionIndirect() : Pointer to matrix is NULL",
                     (char *) NULL);

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

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

    /* [c] : Compute [A] = [L][L]^T */

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

            if( ii == ij) {
                dSum  = 0.0;
                for (ik = 1; ik < ii; ik = ik + 1)
                     dSum += spA->uMatrix.daa[ ii-1 ][ ik-1 ] *
                             spA->uMatrix.daa[ ii-1 ][ ik-1 ];

                if(spA->uMatrix.daa[ ii-1 ][ ii-1 ] > dSum)
                   spA->uMatrix.daa[ii-1][ii-1] = sqrt(spA->uMatrix.daa[ ii-1 ][ ii-1 ] - dSum);
                else
                   FatalError("In CholeskyDecompositionIndirect()",
                              "Matrix is not positive definite",
                              (char *) NULL);
            }

            if( ii > ij) {
                dSum  = 0.0;
                for (ik = 1; ik < ij; ik = ik + 1)
                     dSum += spA->uMatrix.daa[ ii-1 ][ ik-1 ] *
                             spA->uMatrix.daa[ ij-1 ][ ik-1 ];

                spA->uMatrix.daa[ii-1][ij-1] = (spA->uMatrix.daa[ ii-1 ][ ij-1 ] - dSum) /
                                                spA->uMatrix.daa[ij-1][ij-1];

                spA->uMatrix.daa[ij-1][ii-1] = spA->uMatrix.daa[ii-1][ij-1];
            }

            }
       }

       return ( spA );
}


/* 
 *  =================================================================
 *  SetupScaleFactors() : Initialize vector for scale factors. If [A]
 *  is a (nxn) matrix, then spScale will be a (nx1) vector containing
 *  the maximum absolute value in each row of [A].
 *  
 *  Input :   MATRIX spA       -- Pointer to (nxn) matrix [A].
 *  Output :  VECTOR spScale   -- Pointer to (nx1) vector spScale.
 *  =================================================================
 */ 

#ifdef __STDC__
VECTOR *SetupScaleFactors( MATRIX *spA )
#else  /* Start case not STDC */
VECTOR *SetupScaleFactors( spA )
MATRIX *spA;
#endif /* End case not STDC */
{
VECTOR *spScale;
double    dTemp;
int      ii, ij;